# Linear Algebra Manipulations with Cholesky Factors

**URL:** <https://discourse.julialang.org/t/linear-algebra-manipulations-with-cholesky-factors/133899>\
**Category:** Performance\
**Tags:** linearalgebra, sparse\
**Created:** [November 15, 2025, 3:42pm UTC](https://discourse.julialang.org/t/linear-algebra-manipulations-with-cholesky-factors/133899 "2025-11-15T15:42:14Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![gideonsimpson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gideonsimpson/32/1928_2.png) [@gideonsimpson](https://discourse.julialang.org/u/gideonsimpson)\
**Post date:** [November 15, 2025, 3:42pm UTC](https://discourse.julialang.org/t/linear-algebra-manipulations-with-cholesky-factors/133899/1 "2025-11-15T15:42:14Z")

</div>

I’m working on a problem where I need to do some Gaussian random variable sampling, which would naive be implemented as:

```julia-auto
nx = 10
e = ones(nx)
A = spdiagm(-1 => -e[1:end-1], 0 => 2e, 1 => -e[1:end-1])
chol = cholesky(A, perm=1:size(A));
z = randn(nx);
w = chol.U\z;

```

But since I am going to have to generate the `w` variable many times, I would like to optimize the solve step to minimize memory allocations. My first instinct was to do

```julia-auto
F = lu(chol.U);
w = zeros(nx); # preallocate
ldiv!(w, F, z); # in place solve

```

but `chol.U` is the wrong data type. I’m ok with doing a bit of up front computational work here, as I’m going to do many samples, but I’m really trying to avoid any use of dense matrix representations.

---

<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:** [November 15, 2025, 6:51pm UTC](https://discourse.julialang.org/t/linear-algebra-manipulations-with-cholesky-factors/133899/2 "2025-11-15T18:51:20Z")

</div>

> [@gideonsimpson](#):
>
> `A = spdiagm(-1 => -e[1:end-1], 0 => 2e, 1 => -e[1:end-1])`

If your actual matrix is tridiagonal like this, I would tend to use the `SymTridiagonal` type (i.e. a _structured_ sparse matrix) rather than a generic sparse matrix, and exploit the specialized solvers for tridiagonal matrices.

> [@gideonsimpson](#):
>
> `F = lu(chol.U);`

Once you have the cholesky factor, it is already upper triangular, so there is no need to do an additional factorization.

> [@gideonsimpson](#):
>
> `chol = cholesky(A, perm=1:size(A));`

By specifying the identity permutation you may be greatly increasing fill-in, depending on the sparsity pattern. (Also, this code doesn’t run… I think you meant `size(A,1)`.) Generally, I would recommend allowing it to select the permutation for you (and account for that in your linear algebra).

> [@gideonsimpson](#):
>
> `w = chol.U\z;`

If you are only updating `z` and keeping the matrix fixed, this is a sparse triangular solve and hence should be fast.

(It does allocate a new vector for the lhs. Unfortunately, there doesn’t seem to currently be a `ldiv!` implementation for `chol.U`, though it should be straightforward to add one. On the other hand, for large matrices the cost of the allocation of the solution vector shouldn’t be a significant portion of the time.)

What is the typical size `nx` in your actual application? If a tiny size like `nx = 10` is typical, then the best choice of algorithm changes.