I’m working on a problem where I need to do some Gaussian random variable sampling, which would naive be implemented as:
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
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.
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.
Once you have the cholesky factor, it is already upper triangular, so there is no need to do an additional factorization.
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).
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.