Suppose solving Stokes flow. Suppose storing the Jacobian as a BlockArray as in e.g. Incompressible Navier-Stokes equations via DifferentialEquations.jl · Ferrite.jl (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.
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:
- left precondition: we solve (A^\top A) x =A^\top b
- 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. 
These ideas will both be significantly slower than what we actually do (lu factorization)
1 Like
Do you mean this function?
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.
1 Like
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.
2 Likes
Thank you very much, this is exactly what I’m looking for!
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…)
2 Likes
Besides, here are some facts that make me confused.
ldltis not defined for a denseSymmetricmatrix.- But
ldltis defined for a sparseSymmetricmatrix. - In both dense and sparse mode,
ldltis not involved in the default dispatch ofA \ b. (the default choice is\(lu(A), b)instead). - 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
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.
1 Like
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
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.
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:
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 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
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.
4 Likes
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.
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. 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.
4 Likes