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.