Linear Solvers

Solve Ax = b with dense and sparse solvers

Run AllReset

Solving linear systems Ax = b is at the heart of numerical computing — from fitting models to simulating physics. Equana offers two solver families, each suited to different problem sizes and matrix types:

Dense direct solvers — LU, QR, Cholesky, least-squares for small-to-moderate systems
Sparse solvers — iterative methods (GMRES, BiCGSTAB, PCG, LSQR) with convergence diagnostics, plus sparse direct factorizations for large sparse systems

This tutorial walks through both families, explains when to use which solver, and provides runnable examples. Run the cells in order — variables persist across the notebook.

For sparse matrix creation and manipulation, see the Sparse Matrices tutorial. For dense matrix basics, see the Working with Matrices tutorial.

Choosing the Right Solver

The right solver depends on your matrix properties and problem size:

Matrix Properties	Recommended Solver	Backend
Dense, small-to-moderate (n < 5000)	`A \ b` or `linsolve(A, b)`	NumWasm
Dense, overdetermined (m > n)	`lstsq(A, b)`	NumWasm
Dense, SPD, need factorization	`chol(A)` then back-substitution	NumWasm
Sparse, symmetric positive definite	`pcg(A, b)`	SciWasm
Sparse, general nonsymmetric	`gmres(A, b)` or `bicgstab(A, b)`	SciWasm
Sparse, rectangular / least-squares	`lsqr(A, b)`	SciWasm

Direct vs. iterative: Direct solvers (LU, QR, Cholesky) compute an exact answer in a fixed number of operations — O(n³) for dense matrices. Iterative solvers approximate the solution through successive refinement and are much faster for large sparse systems where most entries are zero.

Dense Direct Solvers (NumWasm)

Dense solvers compute the exact solution (up to floating-point precision) by factoring the matrix. They are best for systems up to a few thousand unknowns.

Backslash Operator and linsolve

The backslash operator A \ b, linsolve(A, b), and mldivide(A, b) are three equivalent ways to solve Ax = b directly:

Code [1]Run

# Create a 3x3 system
A = [2, 1, -1; -3, -1, 2; -2, 1, 2]
b = [8; -11; -3]

# Three equivalent ways to solve Ax = b
x1 = A \ b
x2 = linsolve(A, b)
x3 = mldivide(A, b)

# All give the same answer
x1

# Verify: A*x should equal b
norm(A * x1 - b)

LU Decomposition

LU decomposition factors A into a lower triangular matrix L, an upper triangular matrix U, and a permutation matrix P such that P·A = L·U. Once factored, solving for different right-hand sides is cheap (forward and back substitution):

Code [2]Run

# LU decomposition with partial pivoting
(L, U, P) = lu(A)

# Verify the factorization: P*A = L*U
norm(P * A - L * U)

QR Decomposition

QR decomposes A into an orthogonal matrix Q (where Q'·Q = I) and an upper triangular matrix R. QR is numerically more stable than LU and is the basis for least-squares solvers:

Code [3]Run

# QR decomposition: A = Q*R
(Q, R) = qr(A)

# Q is orthogonal: Q'*Q should be the identity
Q' * Q

# Verify: A = Q*R
norm(A - Q * R)

Cholesky Factorization

For symmetric positive definite (SPD) matrices, Cholesky computes an upper triangular R such that R'·R = A. It requires half the work of LU and guarantees numerical stability:

Code [4]Run

# Create a symmetric positive definite matrix
B = [4, 2, 1; 2, 5, 3; 1, 3, 6]

# Cholesky: B = R'*R
R = chol(B)

# Verify the factorization
norm(R' * R - B)

Least-Squares Solutions

When a system is overdetermined (more equations than unknowns, m > n), no exact solution exists. lstsq(A, b) finds the x that minimizes ||Ax - b||₂:

Code [5]Run

# Overdetermined system: 5 equations, 3 unknowns
A_tall = [1, 0, 0; 0, 1, 0; 0, 0, 1; 1, 1, 0; 0, 1, 1]
b_tall = [1; 2; 3; 3.1; 4.9]

# Least-squares: minimize ||A*x - b||
x_ls = lstsq(A_tall, b_tall)
println(x_ls)

The residual, matrix rank, and singular values characterize how well-posed the fit is — rank and svd compute the latter two directly:

Code [6]Run

# Residual of the least-squares fit, plus rank and singular values of A
residual = norm(A_tall * x_ls - b_tall)
println(residual)     # ||A*x - b||

println(rank(A_tall)) # matrix rank

(U, s, V) = svd(A_tall)
println(s)            # singular values (descending)

Sparse Iterative Solvers (SciWasm)

For large sparse systems, direct solvers become impractical — a 10,000×10,000 dense LU costs ~667 billion operations. Iterative solvers approximate the solution through successive refinement, exploiting the sparsity to perform each iteration in O(nnz) time.

Available Solvers

Solver	Best For	Requires
`gmres(A, b)`	General nonsymmetric systems	Nothing special
`bicgstab(A, b)`	General nonsymmetric systems	Nothing special
`pcg(A, b)`	Symmetric positive definite (SPD)	A must be SPD
`lsqr(A, b)`	Rectangular / least-squares	Nothing special
`cgs(A, b)`	General nonsymmetric systems	Nothing special
`bicg(A, b)`	General nonsymmetric (uses A')	Nothing special
`minres(A, b)`	Symmetric indefinite	A must be symmetric
`tfqmr(A, b)`	General nonsymmetric (transpose-free)	Nothing special

All solvers support optional tolerance and iteration-limit parameters. With multiple outputs, they return convergence diagnostics:

(x, flag, relres, iter) = solver(A, b, ...)

Flag	Meaning
0	Converged to the requested tolerance
1	Did not converge within the maximum iterations

Setting Up a Test Problem

We build a classic tridiagonal Laplacian — symmetric positive definite and well-suited for benchmarking:

Code [7]Run

# Build a 100x100 SPD tridiagonal Laplacian: [-1, 2, -1]
# Assembled as (row, col, value) triplets — duplicates are summed
n = 100
rows = zeros(3 * n - 2)
cols = zeros(3 * n - 2)
vals = zeros(3 * n - 2)
idx = 1
for i in 1:n {
  rows[idx] = i
  cols[idx] = i
  vals[idx] = 2.0
  idx = idx + 1
  if i > 1 {
    rows[idx] = i
    cols[idx] = i - 1
    vals[idx] = -1.0
    idx = idx + 1
  }
  if i < n {
    rows[idx] = i
    cols[idx] = i + 1
    vals[idx] = -1.0
    idx = idx + 1
  }
}
T = sparse(rows, cols, vals, n, n)

# Random right-hand side
b_sp = rand(n)

# Matrix properties
println(issparse(T))     # true
println(nnz(T))          # 298 (3n - 2)

GMRES — General Systems

GMRES (Generalized Minimum Residual) works for any square nonsingular matrix. The optional restart parameter limits memory usage by restarting the Krylov subspace after restart iterations:

Code [8]Run

# Solve with GMRES (restart=30, tol=1e-10, maxit=200)
(x_gm, flag_gm, relres_gm, iter_gm) = gmres(T, b_sp, 30, 1e-10, 200)
flag_gm       # 0 = converged
relres_gm     # relative residual
iter_gm       # iterations used

BiCGSTAB — Nonsymmetric Systems

BiCGSTAB (Bi-Conjugate Gradient Stabilized) is another choice for general nonsymmetric systems. It uses less memory than GMRES (no restart needed) and often has smoother convergence:

Code [9]Run

# Solve with BiCGSTAB
(x_bi, flag_bi, relres_bi, iter_bi) = bicgstab(T, b_sp, 1e-10, 200)
flag_bi       # 0 = converged
relres_bi
iter_bi

PCG — Symmetric Positive Definite Systems

Preconditioned Conjugate Gradient is the most efficient solver when the matrix is symmetric and positive definite (SPD). It uses half the work per iteration of BiCGSTAB and is guaranteed to converge for SPD matrices. Many physical problems — diffusion, elasticity, Poisson — produce SPD systems:

Code [10]Run

# PCG — best choice when A is SPD
(x_cg, flag_cg, relres_cg, iter_cg) = pcg(T, b_sp, 1e-10, 200)
println(flag_cg)      # 0 = converged
println(relres_cg)
println(iter_cg)

# Verify the solution
res_cg = norm(reshape(full(spmm(T, x_cg)), n) - b_sp)
println("Residual: ${res_cg}")

LSQR — Least-Squares Problems

lsqr solves overdetermined or underdetermined sparse systems in the least-squares sense (minimizes ||Ax - b||₂). It works on rectangular sparse matrices where direct methods and square solvers cannot be applied:

Code [11]Run

# LSQR for a square system (also works on rectangular)
(x_lsqr, flag_lsqr, relres_lsqr, iter_lsqr) = lsqr(T, b_sp, 1e-10, 200)
println(flag_lsqr)     # 0 = converged

CGS — Conjugate Gradient Squared

CGS is another Krylov method for general nonsymmetric systems. It converges faster than BiCGSTAB in some cases but can have irregular convergence behavior:

Code [12]Run

# Solve with CGS
(x_cgs, flag_cgs, relres_cgs, iter_cgs) = cgs(T, b_sp, 1e-10, 200)
flag_cgs
iter_cgs

Sparse Direct Solvers

For moderate-sized sparse systems where iterative convergence is problematic, sparse direct solvers factorize the matrix. On desktop, Equana uses SuiteSparse (UMFPACK/CHOLMOD) for native-speed factorization:

Code [13]Run

# Sparse Cholesky factorization (A must be SPD)
L = spchol(T)

# Sparse LU factorization with pivoting
(L_lu, U_lu, p_lu) = splu(T)

# Direct solve: automatically selects LU or Cholesky
x_direct = spsolve(T, b_sp)
res_direct = norm(reshape(full(spmm(T, x_direct)), n) - b_sp)
println("Direct-solve residual: ${res_direct}")

Direct or Iterative? A Head-to-Head

For the same 100×100 Laplacian, compare the direct factorization against PCG:

Code [14]Run

# Direct: factorize + solve
tic()
x_d = spsolve(T, b_sp)
t_direct = toc()

# Iterative: PCG (SPD system)
tic()
(x_i, flag_i, relres_i, iter_i) = pcg(T, b_sp, 1e-10, 500)
t_iter = toc()

diff = norm(x_d - x_i)
println("Direct:    ${t_direct}s")
println("Iterative: ${t_iter}s   (${iter_i} iterations)")
println("|x_direct - x_pcg| = ${diff}")

Rules of thumb:

Solve once, moderate size — direct (spsolve): no tuning, exact up to rounding.
Many right-hand sides — factorize once (spchol/splu), reuse the factors.
Very large systems — iterative (pcg for SPD, gmres otherwise): memory stays O(nnz), and you control the accuracy/time trade-off via tol and maxit.
Ill-conditioned systems — watch relres and flag; if convergence stalls, a direct solve or reordering (amd, symrcm) may help.

Quick Reference

Dense Solvers (NumWasm)

Function	Description
`A \ b` / `linsolve(A, b)` / `mldivide(A, b)`	Solve Ax = b directly
`lstsq(A, b)`	Least-squares solution (overdetermined systems)
`(L, U, P) = lu(A)`	LU decomposition with partial pivoting
`(Q, R) = qr(A)`	QR decomposition
`R = chol(A)`	Cholesky factorization (SPD matrices)

Sparse Iterative Solvers (SciWasm)

Function	System Type	Key Parameter
`gmres(A, b, restart, tol, maxit)`	General	`restart` controls memory
`bicgstab(A, b, tol, maxit)`	General	Less memory than GMRES
`pcg(A, b, tol, maxit)`	SPD only	Fastest for SPD
`lsqr(A, b, tol, maxit)`	Rectangular	Least-squares
`cgs(A, b, tol, maxit)`	General	Fast, irregular convergence
`bicg(A, b, tol, maxit)`	General	Requires A'
`minres(A, b, tol, maxit)`	Symmetric	Works for indefinite
`tfqmr(A, b, tol, maxit)`	General	Transpose-free

Sparse Direct Solvers

Function	Description
`spsolve(A, b)`	Auto-select LU or Cholesky and solve
`spchol(A)`	Sparse Cholesky factorization (SPD)
`(L, U, p) = splu(A)`	Sparse LU with pivoting
`sprank(A)`	Structural rank

Summary

This tutorial covered Equana's solver families:

Dense direct solvers (NumWasm) — linsolve, mldivide, the backslash operator A \ b, plus LU, QR, and Cholesky factorizations, and lstsq for overdetermined systems
Sparse iterative solvers (SciWasm) — gmres, bicgstab, pcg, lsqr, cgs, bicg, minres, and tfqmr with convergence diagnostics (flag, relres, iter)
Sparse direct solvers — spsolve, spchol, splu for moderate-sized sparse systems

Next Steps

Explore the Sparse Matrices tutorial for creating and manipulating sparse matrices
Try the Finite Element Method tutorial for a complete FEM workflow with solvers
See the Linear Elasticity tutorial for a structural mechanics solve

Workbench

Clear

No variables in workbench

Next Steps

Sparse Matrices

Create, inspect, and manipulate sparse matrices

Finite Element Method

See solvers in a complete FEM workflow