Help understanding solves with sparse cholesky factors

I recognize there are a few threads about this or something similar already, but I would really appreciate some clarification on something. I am aware that cholesky for sparse matrices permutes, and that you can ostensibly control it with the perm=... keyword. But I don’t get why this doesn’t work:

using LinearAlgebra, SparseArrays

# Make some random sparse pos-def matrix. I'm aware this is not very sparse,
# just a simple example:
tmp = sprandn(256, 256, 0.1)
M   = Symmetric(tmp'*tmp)

# Factorize the matrix in a few ways, making copies and stuff to try and reduce
# the risk of re-used buffers or any other confounding thing.
Mf_dense = cholesky(Symmetric(Matrix(M)))
Mf_1 = cholesky(copy(M), perm=collect(1:256))    # no  permutation
Mf_2 = cholesky(copy(M), perm=collect(256:-1:1)) # yes permutation

# Three solves:
sample_x = rand(256)
s_dense  = Mf_dense.L\sample_x
s_1      = Mf_1.L\sample_x
s_2      = Mf_2.PtL\sample_x

# First test: the un-permuted sparse solve agrees with the dense solve.
@assert s_dense ≈ s_1 "dense solve and un-permuted sparse solve don't match"

# Second test, which fails: the PtL solve of the permuted one agrees with the
# sparse solve (and thus also the dense solve).
@assert s_1 ≈ s_2 "un-permuted sparse solve and permuted sparse solve don't agree"

Unless I’m taking crazy pills, P A P^T = L L^T, then A = (P^T L) (P^T L)^T, which is why the sparse Cholesky struct offers the .PtL field in the first place. So I really don’t get why this second one doesn’t pass.

Can anybody help me here? I can’t even manually play with permutations in such a way to get Mf_1.L\sample_x back. s_1 and s_2 don’t even agree up to a permutation. If you look at hcat(sort(s_1), sort(s_2)), the values are completely different. What’s going on?

If it is a legitimate bug and I had to guess what it is, the ldiv! call internally is treating PtL as a triangular matrix, which it isn’t anymore. Or something like that. But that seems unlikely.

Yes, but P^T L is not lower-triangular, so this is not the same as the Cholesky factor \hat{L} (Mf_1.L) you get from the un-permuted factorization A = \hat{L}\hat{L}^T. P^T L is a factor of some sort, but it’s not a Cholesky factor.

1 Like

Oof, of course. Thank you.