Could this be faster using some decomposition?

General Usage Performance
1

Hi,

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 :

function rand_normalized(M)
   e = abs.(randn(2000))
   return e ./ sqrt(e'M*e)
end

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 ?

2

You could cut down on allocations with

e = randn(2000)
e .= abs.(e)
e .*= 1/sqrt(dot(e, M, e))
return e

which only allocates 1 array rather than the 4 in your code above.

3

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.

4

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.

1 Like
5

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 abs…

Hum… my M does not seem to be sdp… (although it should be!) What could i do for a matrix that is not sdp ?