# Could this be faster using some decomposition?

**URL:** <https://discourse.julialang.org/t/could-this-be-faster-using-some-decomposition/75323>\
**Category:** Performance\
**Tags:** linearalgebra\
**Created:** [January 27, 2022, 9:20pm UTC](https://discourse.julialang.org/t/could-this-be-faster-using-some-decomposition/75323 "2022-01-27T21:20:36Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [January 27, 2022, 9:20pm UTC](https://discourse.julialang.org/t/could-this-be-faster-using-some-decomposition/75323/1 "2022-01-27T21:20:36Z")

</div>

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 :

```julia
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 ?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 27, 2022, 9:56pm UTC](https://discourse.julialang.org/t/could-this-be-faster-using-some-decomposition/75323/2 "2022-01-27T21:56:42Z")

</div>

You could cut down on allocations with

```julia
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.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [January 27, 2022, 9:57pm UTC](https://discourse.julialang.org/t/could-this-be-faster-using-some-decomposition/75323/3 "2022-01-27T21:57:32Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 27, 2022, 10:00pm UTC](https://discourse.julialang.org/t/could-this-be-faster-using-some-decomposition/75323/4 "2022-01-27T22:00:19Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [January 27, 2022, 10:02pm UTC](https://discourse.julialang.org/t/could-this-be-faster-using-some-decomposition/75323/5 "2022-01-27T22:02:00Z")

</div>

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 ?