# What is a recommended way to solve a linear equation?

**URL:** <https://discourse.julialang.org/t/what-is-a-recommended-way-to-solve-a-linear-equation/128420>\
**Category:** New to Julia\
**Tags:** question, linearalgebra, linearsolve\
**Created:** [April 26, 2025, 2:28am UTC](https://discourse.julialang.org/t/what-is-a-recommended-way-to-solve-a-linear-equation/128420 "2025-04-26T02:28:18Z")\
**Posts on this page:** 12\
**Page:** 2

<div class="post-metadata">

**Author:** ![ziolai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ziolai/32/23422_2.png) [@ziolai](https://discourse.julialang.org/u/ziolai)\
**Post date:** [May 1, 2025, 10:27am UTC](https://discourse.julialang.org/t/what-is-a-recommended-way-to-solve-a-linear-equation/128420/21 "2025-05-01T10:27:29Z")

</div>

Suppose solving Stokes flow. Suppose storing the Jacobian as a BlockArray as in e.g. [Incompressible Navier-Stokes equations via DifferentialEquations.jl · Ferrite.jl](https://ferrite-fem.github.io/Ferrite.jl/stable/tutorials/ns_vs_diffeq/) (deliberately omitting the non-linear convective term here for simplicity). What is required to instruct LinSolve to perform Schur complement preconditioning (or any other type of segregated solve as preconditioner)? Use case is Navier-Stokes with chemistry for instance.

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [May 7, 2025, 3:56am UTC](https://discourse.julialang.org/t/what-is-a-recommended-way-to-solve-a-linear-equation/128420/22 "2025-05-07T03:56:19Z")

</div>

> [@stevengj](#):
>
> Is your matrix symmetric? Is it symmetric positive-definite (SPD)? If it is SPD, use sparse Cholesky.

Here is a gap. What if `A` is real-symmetric but not posdef. What algorithm would the built-in solver `A\b` resort to? And I know there is an idea to convert it to posdef:

1. left precondition: we solve (A^\top A) x =A^\top b
2. right precondition: we solve (AA^\top)y = b, then x := A^\top y

would these methods be competent? Allegedly the condition number is squared.

**PS** I want to look into the multiple dispatch of LinearAlgebra.jl, thus I need a tool like Debugger.jl, but I don’t know where is wrong currently → [ERROR: Debugger.jl requires a LineEditREPL type of REPL](https://discourse.julialang.org/t/error-debugger-jl-requires-a-lineeditrepl-type-of-repl/128760). ☹

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [May 7, 2025, 3:59am UTC](https://discourse.julialang.org/t/what-is-a-recommended-way-to-solve-a-linear-equation/128420/23 "2025-05-07T03:59:23Z")

</div>

These ideas will both be significantly slower than what we actually do (lu factorization)

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [May 7, 2025, 4:20am UTC](https://discourse.julialang.org/t/what-is-a-recommended-way-to-solve-a-linear-equation/128420/24 "2025-05-07T04:20:50Z")

</div>

Do you mean this function?

> <https://github.com/JuliaSparse/SparseArrays.jl/blob/16bbcbc78831dcdc7328ec00690a32ca12a90f5e/src/solvers/cholmod.jl#L1934-L1946>

I find the dispatch a little bit involved. Since there are 2 properties: Symmetric and Sparse.  
I mean, is there a concise type `SymSparse`, just like `SymTridiagonal`? Otherwise the style looks not very neat.

---

<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 7, 2025, 12:14pm UTC](https://discourse.julialang.org/t/what-is-a-recommended-way-to-solve-a-linear-equation/128420/25 "2025-05-07T12:14:10Z")

</div>

> [@Oscar\_Smith](#):
>
> These ideas will both be significantly slower than what we actually do (lu factorization)

Slower and much less accurate, since squaring the condition number will also double the number of digits you lose to roundoff errors.

You can use `bunchkaufman(A) \ b` to have an LU-like factorization that exploits symmetry but does not require positive-definiteness, but this is only implemented for dense matrices.

---

<div class="post-metadata">

**Author:** ![fdekerme](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fdekerme/32/43574_2.png) [@fdekerme](https://discourse.julialang.org/u/fdekerme)\
**Post date:** [May 7, 2025, 3:32pm UTC](https://discourse.julialang.org/t/what-is-a-recommended-way-to-solve-a-linear-equation/128420/26 "2025-05-07T15:32:16Z")

</div>

Thank you very much, this is exactly what I’m looking for!

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [May 7, 2025, 4:07pm UTC](https://discourse.julialang.org/t/what-is-a-recommended-way-to-solve-a-linear-equation/128420/27 "2025-05-07T16:07:32Z")

</div>

> [@stevengj](#):
>
> You can use `bunchkaufman(A) \ b` to have an LU-like factorization that exploits symmetry but does not require positive-definiteness, but this is only implemented for dense matrices.

With a sparse matrix you can do `ldlt(A) \ b`, which ~~I think is the same concept, just different name.~~ EDIT: No, they’re different: in Bunch-Kaufman, D is not necessarily diagonal, see below.

~~(It’s a bit of a mystery why LinearAlgebra.jl defines these as two separate decompositions. If you read the `ldlt(::SymTridiagonal)` docstring, it sounds like the difference is that `bunchkaufman` uses pivoting while `ldlt` does not, but `ldlt(::SparseMatrixCSC)` _does_ in fact use pivoting…)~~

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [May 7, 2025, 11:54pm UTC](https://discourse.julialang.org/t/what-is-a-recommended-way-to-solve-a-linear-equation/128420/28 "2025-05-07T23:54:10Z")

</div>

Besides, here are some facts that make me confused.

1. `ldlt` is **not** defined for a dense `Symmetric` matrix.
2. But `ldlt` is defined for a sparse `Symmetric` matrix.
3. In both dense and sparse mode, `ldlt` is **not involved** in the default dispatch of `A \ b`. (the default choice is `\(lu(A), b)` instead).
4. For a sparse `A`, `\(ldlt(A), b)` works.

Then the question is: Should I use `\(ldlt(A), b)` rather than `\(lu(A), b)`, if `A` is sparse, symmetric, not posdef?

Additionally, the following doc is obscure to me

```julia
  ldlt(A::SparseMatrixCSC; shift = 0.0, check = true, perm=nothing) -> CHOLMOD.Factor

  Compute the LDL' factorization of a sparse matrix A. A must be a SparseMatrixCSC or a
  Symmetric/Hermitian view of a SparseMatrixCSC. Note that even if A doesn't have the type tag, it     
  must still be symmetric or Hermitian. A fill-reducing permutation is used.

```

Why should it enforce “must still be”?  
If I have an asymmetric `A` and I do `A = Symmetric(A)`, then my intended meaning is _essentially_ `A = Matrix(Symmetric(A))`. Therefore, I don’t think this interface is reasonable.

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [May 8, 2025, 2:35am UTC](https://discourse.julialang.org/t/what-is-a-recommended-way-to-solve-a-linear-equation/128420/29 "2025-05-08T02:35:30Z")

</div>

> [@danielwe](#):
>
> > [@stevengj](#):
> >
> > You can use `bunchkaufman(A) \ b` to have an LU-like factorization that exploits symmetry but does not require positive-definiteness, but this is only implemented for dense matrices.
> 
> With a sparse matrix you can do `ldlt(A) \ b`, which I think is the same concept, just different name.

No, there’s a difference: in Bunch-Kaufman, the D is block diagonal with 1x1 and 2x2 blocks, and is represented as a `Tridiagonal` (seems like this could have been `SymTridiagonal`?). This is not mentioned in the documentation. In LDLt, D is diagonal, however the factorization may fail even if the matrix is nonsingular.

EDIT: documentation PR to clarify `bunchkaufman` submitted: [Clarify the structure of Bunch-Kaufman factors by danielwe · Pull Request #1341 · JuliaLang/LinearAlgebra.jl · GitHub](https://github.com/JuliaLang/LinearAlgebra.jl/pull/1341)

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [May 8, 2025, 4:18am UTC](https://discourse.julialang.org/t/what-is-a-recommended-way-to-solve-a-linear-equation/128420/30 "2025-05-08T04:18:39Z")

</div>

> [@WalterMadelim](#):
>
> - `ldlt` is **not** defined for a dense `Symmetric` matrix.
> - But `ldlt` is defined for a sparse `Symmetric` matrix.

This follows from what the underlying libraries provide: LAPACK for dense and CHOLMOD for sparse. Sparse linear algebra uses different algorithms and implementations compared to dense, so you can’t expect one-to-one feature parity. As explained above, for dense matrices, there’s the similar but not identical `bunchkaufman`.

> [@WalterMadelim](#):
>
> - In both dense and sparse mode, `ldlt` is **not involved** in the default dispatch of `A \ b`. (the default choice is `\(lu(A), b)` instead).
> - For a sparse `A`, `\(ldlt(A), b)` works.

Regarding the last point, it’s not always the case that `ldlt(A)` works, even if `A` is nonsingular. It almost always fails on a randomly sampled matrix such as `ldlt(Symmetric(sprandn(1000, 1000, 0.025)))`. Quoting [Wikipedia](https://en.wikipedia.org/wiki/Cholesky_decomposition#LDL_decomposition):

> Some indefinite matrices for which no Cholesky decomposition exists have an LDL decomposition with negative entries in **D** : it suffices that the first _n_ − 1 leading principal minors of **A** are non-singular.[12]

That’s probably why `ldlt` is not used by default—it’s not considered worthwhile to try and possibly fail with `ldlt` when `lu` provides a robust alternative.

It’s also not worthwhile for `\` to spend too much time choosing the ideal factorization, since it will only be used once and then thrown away. I think that’s why `factorize` tries a little harder—the resulting factorization object can be reused multiple times, so it’s worth the effort to check a few more properties. However, if you really care about exploiting structure, you should use your problem-specific knowledge and choose the appropriate factorization yourself. For example, I’m sure there exist nonrandom symmetric sparsity structures, perhaps from PDE discretizations, where you can be confident that the leading principal minors are well-conditioned but the matrix is indefinite. If you know this is the case, reaching for `ldlt` makes sense.

* * *

All that said, sparse symmetric/Hermitian `\` _does_ attempt a Cholesky factorization, with LU as a fallback if Cholesky fails (that is, if the matrix is not positive definite). It might seem like it would be better to replace `cholesky` with `ldlt` here, since `ldlt` works for a larger class of matrices than `cholesky` and they are otherwise very similar mathematically and implemented by the same library. However, `cholesky` is implemented using a high-performance supernodal algorithm, while `ldlt` uses a non-supernodal (simplicial) algorithm. I can’t tell you exactly what that means, but it’s clear from the CHOLMOD user guide that the supernodal algorithm is substantially faster for large matrices (see [SuiteSparse/CHOLMOD/Doc/CHOLMOD\_UserGuide.pdf at v7.10.3 · DrTimothyAldenDavis/SuiteSparse · GitHub](https://github.com/DrTimothyAldenDavis/SuiteSparse/blob/v7.10.3/CHOLMOD/Doc/CHOLMOD_UserGuide.pdf) at the top of page 107). Thus, trying `cholesky` before falling back to LU makes some sense even if trying `ldlt` is not considered worthwhile.

EDIT: documentation PR to clarify performance vs applicability tradeoff between `cholesky` and `ldlt`: [Clarify pros, cons and limitations of Cholesky and LDLt by danielwe · Pull Request #621 · JuliaSparse/SparseArrays.jl · GitHub](https://github.com/JuliaSparse/SparseArrays.jl/pull/621)

* * *

> [@WalterMadelim](#):
>
> If I have an asymmetric `A` and I do `A = Symmetric(A)`, then my intended meaning is _essentially_ `A = Matrix(Symmetric(A))`.

I think you’re misunderstanding the docstring. If you do `A = Symmetric(A)`, then you have applied the type tag, so all is good—`ldlt` will see your matrix as identical to `SparseMatrixCSC(Symmetric(A))`, which is what you wanted. However, if you supply a plain `B::SparseMatrixCSC`, _not_ wrapped in symmetric, and `B` is not exactly symmetric element-by-element, then `ldlt` will error. That’s what the docstring is explaining.

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [May 8, 2025, 6:12am UTC](https://discourse.julialang.org/t/what-is-a-recommended-way-to-solve-a-linear-equation/128420/31 "2025-05-08T06:12:57Z")

</div>

> [@danielwe](#):
>
> That’s what the docstring is explaining.

The expression of this sentence is not very proper, particularly the word “even”.

> Note that even if A doesn’t have the type tag, it must still be symmetric or Hermitian.

I would say:

- Note that **if** A doesn’t have a type tag, **then** _itself_ must be symmetric or Hermitian.

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [May 8, 2025, 6:48am UTC](https://discourse.julialang.org/t/what-is-a-recommended-way-to-solve-a-linear-equation/128420/32 "2025-05-08T06:48:34Z")

</div>

I agree that the wording could be improved. When you have specific issues and suggestions for improvement, don’t hesitate to submit them on github! (I just submitted my first pull request to SparseArrays.jl myself, see the link in my previous post.) Here’s a link to the specific lines you’re suggesting to change: [SparseArrays.jl/src/solvers/cholmod.jl at main · JuliaSparse/SparseArrays.jl · GitHub](https://github.com/JuliaSparse/SparseArrays.jl/blob/main/src/solvers/cholmod.jl#L1716-L1717). If you want, you can submit a suggested change right there in the browser by clicking the little pencil icon in the upper right corner.

[Previous page](https://discourse.julialang.org/t/what-is-a-recommended-way-to-solve-a-linear-equation/128420.md?page=1)