# Help understanding solves with sparse cholesky factors

**URL:** <https://discourse.julialang.org/t/help-understanding-solves-with-sparse-cholesky-factors/81165>\
**Category:** Numerics\
**Created:** [May 16, 2022, 6:58pm UTC](https://discourse.julialang.org/t/help-understanding-solves-with-sparse-cholesky-factors/81165 "2022-05-16T18:58:22Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [May 16, 2022, 6:58pm UTC](https://discourse.julialang.org/t/help-understanding-solves-with-sparse-cholesky-factors/81165/1 "2022-05-16T18:58:22Z")

</div>

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:

```julia

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.

---

<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:** [May 16, 2022, 7:23pm UTC](https://discourse.julialang.org/t/help-understanding-solves-with-sparse-cholesky-factors/81165/2 "2022-05-16T19:23:44Z")

</div>

> [@cgeoga](#):
>
> Unless I’m taking crazy pills, P A P^T = L L^T, then A = (P^T L) (P^T L)^T,

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.

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [May 16, 2022, 7:28pm UTC](https://discourse.julialang.org/t/help-understanding-solves-with-sparse-cholesky-factors/81165/3 "2022-05-16T19:28:55Z")

</div>

Oof, of course. Thank you.