This is more maths than Julia optimisation, but I do not recall if there is a trick in the linear algebra book to make this more efficiently.
I have a covariance matrix
M, of size
(2000,2000), and i want to normalize vectors w.r.t. this matrix :
e = abs.(randn(2000))
return e ./ sqrt(e'M*e)
I am using this function inside a loop and therefore if there is some precomputations depending on
M that we can do, lets do them.
Any idea ?
You could cut down on allocations with
e = randn(2000)
e .= abs.(e)
e .*= 1/sqrt(dot(e, M, e))
which only allocates 1 array rather than the 4 in your code above.
Yeah i already have done technical stuff, this runs through
LoopVectorization, so no allocations etc… I was more thinking about a math trick. I only presented to you a ‘readable’ version.
If you precompute a Cholesky factorization M = R^*R (since M must be SPD), then you could use
norm(R*x) instead. This saves half of the flops because half of R is zero.
That’s it. Thanks a lot this is perfect.
Can’t we simplify even more using the fact that
e is gaussian ? Hum maybe not due to the
M does not seem to be sdp… (although it should be!) What could i do for a matrix that is not sdp ?