What is a recommended way to solve a linear equation?

[WalterMadelim]

I wonder what is a fair way to solve (maybe approximately) a linear equation system Ax = b, at this time. The A may be very large scale and sparse. I might want to use an iterative solution approach (e.g. CG, MINRES etc.), which is allegedly more appropriate than the direct (factorization) methods.

I noticed LinearSolve.jl, but I don’t think his documentation teaches people well.
I am not aware of how to use that package properly after reading.
e.g. Solving Linear Systems in Julia · LinearSolve.jl

What does p, t mean?
And I think he should stress the usage of matrix-free operator, since this is the novelty of this package. Otherwise I can use A\b directly.
This is to say, I think there should be a proper example, e.g. generate a 100-by-100 sparse A, and showcase how to prepare the input arg A—whether it should be a SparseMatrixCSC type, or a Function that takes x as an argument and returns an image vector. And how to prepare the b vector—should it be a SparseVector or something.

Another example, the introduction here Accelerating your Linear Solves · LinearSolve.jl

A Krylov subspace method with proper preconditioning will be better than direct solvers when the matrices get large enough. You could always precondition a sparse matrix with iLU as an easy choice, though the tolerance would need to be tuned in a problem-specific way.

What does “iLU” mean? This terms comes without other information. I guess there should be a link or something associated.

And there seems to be multiple solver packages, I don’t know which one should I choose. e.g. I want to use MINRES, there is IterativeSolversJL_MINRES and KrylovJL_MINRES. Can someone tell me the status quo?

Thanks, I’d appreciate.

[ufechner7]

@WalterMadelim

What does “iLU” mean?

In the given sentence, iLU refers to incomplete LU factorization, a preconditioning technique used in solving large sparse linear systems. Here’s the breakdown:

1. Preconditioning Purpose:
   iLU approximates the LU factorization of a matrix while preserving sparsity, making it computationally efficient. It reduces fill-in (uncontrolled non-zero entries in LU factors) compared to full LU decomposition

• Role in Krylov Subspace Methods:
  When paired with iterative Krylov subspace methods (e.g., GMRES, conjugate gradient), iLU accelerates convergence by transforming the original system into one with better numerical properties. This is crucial for large matrices where direct solvers become infeasible due to memory constraints.
• Tolerance Tuning:
  The “tolerance” refers to parameters controlling the accuracy of the iLU approximation (e.g., drop thresholds for small entries, level of fill). These must be adjusted per problem to balance preconditioner quality and computational cost.

[ufechner7]

@WalterMadelim

What does p, t mean?

the parameters p and t have the following meanings:

p:

• p stands for parameters of the system or operator.
• These parameters can represent any additional, possibly time-varying or system-dependent values that influence the behavior of the operator or system matrix A.
• In many control and modeling contexts, p can be a scalar or a vector and may represent physical parameters, scheduling variables (in LPV systems), or other coefficients required to define the system at a given instance

t:

• t stands for time.
• This allows the operator or system matrix A to be explicitly dependent on time, accommodating time-varying systems or dynamics .

Example from Related Documentation

A typical function signature in such modeling frameworks is:

A(u, p, t)

where:

• u is the input or vector being multiplied,
• p is the parameters,
• t is the current time

This is consistent with common practices in control and simulation libraries, where functions for system dynamics or measurement are written as f(x, u, p, t), with p and t providing flexibility for parameterized and time-varying models

This are AI generated answers, someone with more practical experiance in this field could probably provide a better example.

To summarize: Yes, there is room for improvement for the documentation (pull requests welcome), but AI can already answer most of your questions. I use perplexity .

Thanks for sharing these interesting questions, and welcome in the community!

[WalterMadelim]

Thanks for your instructive answer.
But I guess using LinearAlgebra is also required? Otherwise the name I is undefined.
And one more issue is that

julia> LinearMap
ERROR: UndefVarError: `LinearMap` not defined in `Main`
Suggestion: check for spelling errors or missing imports.

julia> Linear # hit <tab>, I only see these names
LinearAlgebra         LinearAliasSpecifier
LinearIndices         LinearProblem
LinearSolve           LinearSolveAdjoint
LinearSolveFunction

[karei]

@WalterMadelim

But I guess using LinearAlgebra is also required?

Yes. using LinearAlgebra.

@WalterMadelim

And one more issue is that

] add LinearMaps to add this package. And then using LinearMaps.

[WalterMadelim]

Yes, there is a package.
But then I encouter

ERROR: All methods for the model function `f` had too few arguments. For example,
an ODEProblem `f` must define either `f(u,p,t)` or `f(du,u,p,t)`. This error
can be thrown if you define an ODE model for example as `f(u,t)`. The parameters
`p` are not optional in the definition of `f`! For more information on the required
number of arguments for the function you were defining, consult the documentation
for the `SciMLProblem` or `SciMLFunction` type that was being constructed.

For example, here is the no parameter Lorenz equation. The two valid versions
are out of place:

function lorenz(u,p,t)
  du1 = 10.0*(u[2]-u[1])
  du2 = u[1]*(28.0-u[3]) - u[2]
  du3 = u[1]*u[2] - 8/3*u[3]
  [du1,du2,du3]
 end
 u0 = [1.0;0.0;0.0]
 tspan = (0.0,100.0)
 prob = ODEProblem(lorenz,u0,tspan)

function lorenz!(du,u,p,t)
  du[1] = 10.0*(u[2]-u[1])
  du[2] = u[1]*(28.0-u[3]) - u[2]
  du[3] = u[1]*u[2] - 8/3*u[3]
 end
 u0 = [1.0;0.0;0.0]
 tspan = (0.0,100.0)
 prob = ODEProblem(lorenz!,u0,tspan)

Offending function: f
Methods:
# 1 method for callable object:
 [1] (L::LinearMap)(x::AbstractVector)
     @ K:\judepot1114\packages\LinearMaps\KV78T\src\LinearMaps.jl:141

Stacktrace:
 [1] isinplace(f::FunctionMap{…}, inplace_param_number::Int64, fname::String, iip_preferred::Bool; has_two_dispatches::Bool, isoptimization::Bool, outofplace_param_number::Int64)
   @ SciMLBase K:\judepot1114\packages\SciMLBase\cwnDi\src\utils.jl:280
 [2] isinplace (repeats 2 times)
   @ K:\judepot1114\packages\SciMLBase\cwnDi\src\utils.jl:246 [inlined]
 [3] LinearProblem(::FunctionMap{Float64, var"#3#4", Nothing, true}, ::Vector{Float64}; kwargs::@Kwargs{})
   @ SciMLBase K:\judepot1114\packages\SciMLBase\cwnDi\src\problems\linear_problems.jl:76
 [4] LinearProblem(::FunctionMap{Float64, var"#3#4", Nothing, true}, ::Vector{Float64})
   @ SciMLBase K:\judepot1114\packages\SciMLBase\cwnDi\src\problems\linear_problems.jl:70
 [5] top-level scope
   @ REPL[18]:1
Some type information was truncated. Use `show(err)` to see complete types.

[gdalle]

@ufechner7

This are AI generated answers, someone with more practical experiance in this field could probably provide a better example.

@ufechner7 did you at least check that the AI-generated answers could run and were correct? If not, what is the benefit of posting them here?

[karei]

I spent a couple of hours making this mindmap with the help of chatGPT, as I often struggle with linear system solving methods myself. I hope this mindmap will help you. Anyone is more than welcome to point out errors in the diagram or improve this diagram.

Download:
Solving linear system.drawio

[ufechner7]

Sorry, I was too much in the hurry. The ai renerated example doesn’t work. So some human expert is needed to provide a working example, and that is not me.

[stevengj]

@WalterMadelim

The A may be very large scale and sparse. I might want to use an iterative solution approach (e.g. CG, MINRES etc.), which is allegedly more appropriate than the direct (factorization) methods.

Always try a sparse-direct solver first. They have the advantage of being the most straightforward solution: they basically either work or run out of memory. You can try a couple of different options (e.g. the built in A \ b, or Pardiso.jl, or MUMPS.jl).

Is your matrix symmetric? Is it symmetric positive-definite (SPD)? If it is SPD, use sparse Cholesky.

Otherwise, you are stuck with iterative algorithms, which require a lot more experimentation to make them work well. If the matrix is SPD, you should use the conjugate-gradient algorithm (CG). Otherwise, I would start with GMRES(n) or BiCGStab(ℓ), and for both you’ll need to start by experimenting with the restarting period (n) or subspace size (ℓ), which trade off memory and computational cost for convergence rates. (If the matrix is symmetric but not PD, you could try e.g. MINRES or SYMMLQ too.)

If one of these iterative algorithms converges satisfactorily quickly, congratulations! Otherwise you’ll need to experiment with preconditioners, which are a problem-dependent art. Basically, a preconditioner can be thought of as a crude approximate inverse for A that can be applied quickly, and it depends heavily on where your matrix comes from. There are various generic choices you can try, e.g. Jacobi (the inverse of the diagonal or diagonal blocks of A), incomplete LU (sparse LU factors where you discard lots of elements), algebraic multigrid, and so on, but in many cases people use something more problem-specific. In some fields, a good preconditioner is known, while in others it is a research problem.

[ChrisRackauckas]

We’re continually improving the documentation. Here’s a pull request with your comments included:

Address comments of Discourse

Keep shooting and we’ll keep improving.

@WalterMadelim

And there seems to be multiple solver packages, I don’t know which one should I choose. e.g. I want to use MINRES, there is IterativeSolversJL_MINRES and KrylovJL_MINRES. Can someone tell me the status quo?

The solvers page goes into detail with recommended methods in different cases:

Linear System Solvers · LinearSolve.jl

@WalterMadelim

Yes, there is a package.
But then I encouter

LinearMaps is a wrong suggestion there. It should be SciMLOperators.jl. I’ll add this to that doc page, it should’ve had a separate section.

[WalterMadelim]

@stevengj

Is your matrix symmetric? Is it symmetric positive-definite (SPD)? If it is SPD, use sparse Cholesky.

But what if, say, I get a A, e.g.

A = some_procedure(some_arg); # being a Matrix{Float64}

How do the solver (e.g. the built-in \) know that A is symmetric? Does the solver has a procedure to detect?

Moreover, assume A is symmetric, e.g. assume I will call

A = LinearAlgebra.Symmetric(A)

At this time I don’t know whether A is positive definite (it may be, or may not be).
Then what will the solver react?

• Will it first start a procedure to determine whether A is positive definite? Is this procedure computationally intensive? What method will it use? Does it calculate its minimum eigenvalue and check its positiveness (we discussed this last time).
• Or will it attempt to start linear solving directly?

[stevengj]

@WalterMadelim

Does the solver has a procedure to detect?

No. Symmetric matrices and definite matrices don’t happen by accident. You should know your problem well enough to know whether you have these nice properties.

[ChrisRackauckas]

It’s up now: Getting Started with Solving Linear Systems in Julia · LinearSolve.jl let me know if this reads more clearly or if it needs further edits. I’m leaving the matrix-free part improvements to a further update since there’s a coming change to the SciMLOperators interface that I will want to get in before fully committing this doc to.

[danielwe]

@WalterMadelim

How do the solver (e.g. the built-in \) know that A is symmetric?

@stevengj

@WalterMadelim

Does the solver has a procedure to detect?

No.

This is incorrect. When A is a Matrix, the basic A \ b solver from LinearAlgebra.jl does check for symmetry/hermiticity, (unit) triangularity, and diagonality. This is used to select an appropriate matrix factorization. You can see how this is implemented in Julia 1.11 here: LinearAlgebra.jl/src/dense.jl at release-1.11 · JuliaLang/LinearAlgebra.jl · GitHub (It’s been refactored for Julia 1.12, but the logic is the same.)

However, it does not detect positive definiteness.

Symmetry/hermiticity (but not positive definiteness) is in fact detected for sparse matrices, see SparseArrays.jl/src/linalg.jl at release-1.11 · JuliaSparse/SparseArrays.jl · GitHub and SparseArrays.jl/src/sparsematrix.jl at release-1.11 · JuliaSparse/SparseArrays.jl · GitHub

The generic fallback for other AbstractMatrix subtypes does not check for symmetry/hermiticity, only triangularity and diagonality: LinearAlgebra.jl/src/generic.jl at release-1.11 · JuliaLang/LinearAlgebra.jl · GitHub

EDIT: fixed an incorrect link for sparse matrices

EDIT 2: struck out incorrect claims, see below

[WalterMadelim]

I left a comment on Github. The #600 Pull Request.

[liuyxpp]

I am looking forward to the last section, namely Using Matrix-Free Operators, which I believe will be extremely helpful for lots of people to get started with SciML.

[WalterMadelim]

Thanks for the links.

@danielwe

However, it does not detect positive definiteness.

I notice that the latest (dev) doc has been revised, compared with v1.11.5, about this point.

[danielwe]

EDIT: this is correct for factorization, but irrelevant for \, see below.

You’re right, I totally missed that branch when scanning the code, but up to and including 1.11 it does indeed check for positive definiteness: LinearAlgebra.jl/src/dense.jl at release-1.11 · JuliaLang/LinearAlgebra.jl · GitHub.

This behavior was inconsistently implemented (it would do the pos-def check when passed a regular Matrix that happened to be Hermitian, but not when passed a Hermitian{T,Matrix{T}} directly) and has been removed in 1.12: LinearAlgebra.jl/src/dense.jl at release-1.12 · JuliaLang/LinearAlgebra.jl · GitHub (PR here: Don't try `cholesky` for `Hermitian` in `factorize` by dkarrasch · Pull Request #54775 · JuliaLang/julia · GitHub).

[danielwe]

I need to issue a mea culpa here. In previous posts, I’ve been pointing to the structure identification in the factorization method, since I was convinced I had seen a specialized \ method for Matrix or at least Matrix{<:BlasFloat} that either went through factorization or mirrored its implementation. However, such a method does not exist. All \(::Matrix, ::Union{Matrix,Vector}) calls go straight to the generic method at LinearAlgebra.jl/src/generic.jl at release-1.11 · JuliaLang/LinearAlgebra.jl · GitHub. This method checks for triangular and diagonal structure, but nothing else—no symmetry/hermiticity or positive definiteness, regardless of Julia version; those checks are only performed when you call factorization directly.

For sparse matrices, however, symmetry/hermiticity (but not positive definiteness) is detected by \.