A Visual Guide to Solvers

1. Direct solvers: matrix factorization

Triangular systems are easy to solve: find the single unknown in the first equation, then substitute it into the next. Direct methods work by factoring a general matrix into triangular pieces.

In the diagrams below, coloured squares show potentially nonzero entries and pale squares show structural zeros.

LU: record elimination steps

LU splits a row-permuted square matrix into lower and upper triangular factors.

Gaussian elimination removes entries below the diagonal. LU factorization records these steps as \(PA = LU\), where \(P\) reorders rows for numerical stability. Solving \(Ax = b\) reduces to two triangular solves: \(Ly = Pb\), then \(Ux = y\).

Use solve_lu() for general nonsingular square matrices. Full pivoting also permutes columns; the solver handles all permutations automatically.

QR: build perpendicular directions

QR splits a matrix into an orthogonal factor Q and an upper triangular factor R.

In \(A = QR\), the columns of \(Q\) are orthonormal, i.e., perpendicular unit vectors. Rotating into \(Q\)’s coordinates preserves lengths, leaving an upper-triangular system in \(R\). Pivoted variants also reorder columns.

Use solve_qr() for square systems or rectangular least squares. When there are more equations than unknowns, it finds the best fit by minimizing \(\Vert{}Ax - b\Vert{}_2^2\). Full column rank guarantees a unique solution.

Cholesky: use symmetry to halve the work

Cholesky splits a symmetric positive definite matrix into L and its transpose.

A symmetric positive definite (SPD) matrix satisfies \(A = A^\top\) and \(z^\top Az > 0\) for every nonzero \(z\). It factors as \(A = LL^\top\), where a single triangular factor serves as both halves—requiring roughly half the compute and memory of LU.

Use solve_chol(A, b, pivot = 0) for the standard \(LL^\top\) form. By default, solve_chol() computes \(LDL^\top\), which inserts a diagonal factor \(D\). While SPD input is the safe baseline for both, do not assume the default \(LDL^\top\) can solve arbitrary indefinite or singular symmetric systems.

2. Iterative solvers: improve a guess

Iterative methods start with an initial guess \(x_0\) and progressively update it. At each step, the residual \(r_k = b - Ax_k\) measures how far the current output is from \(b\).

For SPD matrices, solving \(Ax = b\) is equivalent to finding the bottom of a quadratic valley: \[ f(x) = \tfrac12x^\top Ax - b^\top x. \]

CG follows two steps from an initial guess to the centre of elliptical contours; residual norms decrease from 19.82 to approximately zero in this example.

The ellipses trace contours of equal height. Conjugate gradient (CG) steps along \(A\)-conjugate directions so it never undoes earlier progress, reaching this 2D minimum in two steps. In exact arithmetic, CG converges in at most \(n\) steps for an \(n\)-dimensional SPD system. In practice, conditioning and roundoff govern convergence, and the residual norm may not decrease at every step.

Despite its name, solve_cg() defaults to BiCGSTAB – specify type explicitly when you want CG or LSCG. Inputs are automatically converted from dense to sparse format. Because each step only requires matrix-vector products, iterative solvers avoid fill-in, where direct factors become much denser than \(A\).

A preconditioner reshapes the problem to speed up convergence—like rounding out a narrow, elongated valley. The default diagonal scaling (precond = 1) is a helpful baseline. Use tol and iter to control stopping criteria; always check convergence and the final residual.

3. Spectral methods: find key directions

Instead of solving for a single input \(x\), spectral methods ask how the machine \(A\) acts on space as a whole: which directions does it stretch the most, which does it squash near zero, and which pass through unchanged? Answering this reveals whether \(Ax = b\) is sensitive to noise, redundant, or reducible to a much smaller problem.

Arnoldi and Lanczos: project a large matrix onto a small one

Repeatedly multiplying a starting vector \(v\) by \(A\) explores a Krylov subspace: \[ \mathcal K_k(A,v) = \operatorname{span}\{v,Av,\ldots,A^{k-1}v\}. \]

arnoldi() builds an orthonormal basis \(Q\) for this subspace and projects \(A\) onto a smaller upper-Hessenberg matrix with zeros below the first subdiagonal.

Arnoldi builds a Krylov basis from a starting vector and projects a large sparse matrix onto a small upper-Hessenberg matrix.

lanczos() exploits symmetry to produce a tridiagonal projection using a shorter recurrence. By default, sanic uses reorthogonalisation to keep the basis accurate.

Lanczos projects a symmetric sparse matrix onto a small tridiagonal matrix using a Krylov basis.

The eigenvalues of these projections approximate the most prominent parts of \(A\)’s spectrum rather than every eigenvalue. Set iter below nrow(A) (the default) to build a smaller projection. Both functions convert dense inputs to sparse format and return the projection, basis Q, and optional eigenvalues—not the eigenvectors of \(A\). For symmetric inputs, arnoldi() automatically switches to Lanczos.

4. Which should I try?

Treat these as starting points—size, fill-in, rank, and conditioning all matter. For a quick default, solve2() automatically picks a direct solver: QR for rectangular matrices, Cholesky for symmetric matrices (falling back to LU if needed), and LU otherwise.

5. One system, several routes

Consider the SPD system

\[ A = \begin{bmatrix} 4 & 1 \\ 1 & 3 \end{bmatrix}, \qquad b = \begin{bmatrix} 1 \\ 2 \end{bmatrix}, \qquad x = \begin{bmatrix} 1/11 \\ 7/11 \end{bmatrix}. \]

LU, QR, Cholesky, and CG all recover \(x\) via different routes:

A <- matrix(c(4, 1, 1, 3), nrow = 2)
b <- c(1, 2)
fits <- cbind(
  LU = solve_lu(A, b),
  QR = solve_qr(A, b),
  Cholesky = solve_chol(A, b),
  CG = solve_cg(sparsify(A), b, type = "CG", tol = 1e-10)
)
fits
#>            [,1]       [,2]       [,3]       [,4]
#> [1,] 0.09090909 0.09090909 0.09090909 0.09090909
#> [2,] 0.63636364 0.63636364 0.63636364 0.63636364

# Check how closely each result reproduces b.
r <- A %*% fits - matrix(b, nrow = 2, ncol = ncol(fits))
sqrt(colSums(r^2)) / sqrt(sum(b^2))
#> [1] 0.000000e+00 2.220446e-16 0.000000e+00 0.000000e+00

A small relative residual is reassuring, though for ill-conditioned matrices it does not guarantee a small error in \(x\). Here, all four methods match the exact answer, c(1, 7) / 11.

For a rectangular example, we fit an intercept and slope to noisy observations:

X <- cbind(intercept = 1, time = 0:4)
y <- c(1, 2, 2, 4, 4)
cbind(
  QR = solve_qr(X, y),
  LSCG = solve_cg(sparsify(X), y, type = "LSCG", tol = 1e-10)
)
#>      [,1] [,2]
#> [1,]  1.0  1.0
#> [2,]  0.8  0.8

A nonzero residual is expected here – we want the best fit, not an exact solve. On a tiny problem like this, direct factorization is effortless; on a sparse matrix with millions of rows, an iterative solver shines.

In short: direct methods factor the matrix, iterative methods refine a guess, and spectral methods reveal the geometry. See ?solve2, ?solve_cg, and ?arnoldi for full parameter details.

Further reading

• Trefethen, L. N., & Bau, D. (1997). Numerical Linear Algebra. SIAM. — A concise, geometric introduction to QR, SVD, conditioning, and Krylov subspaces. https://doi.org/10.1137/1.9780898719574.
• Shewchuk, J. R. (1994). An Introduction to the Conjugate Gradient Method Without the Agonizing Pain. — A classic visual guide to quadratic valleys, contours, and conjugate directions.
• Golub, G. H., & Van Loan, C. F. (2013). Matrix Computations (4th ed.). Johns Hopkins University Press. — The standard reference for matrix factorizations and spectral algorithms.
• Halko, N., Martinsson, P.-G., & Tropp, J. A. (2011). Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions. SIAM Review, 53(2), 217–288. — A foundational introduction to randomized low-rank approximation and probabilistic algorithms. https://doi.org/10.1137/090771806.
• Guennebaud, G., Jacob, B., et al. Eigen C++ Library Documentation. — Implementation details for the C++ solvers underlying sanic. https://libeigen.gitlab.io/.
• Qiu, Y., & Mei, J. (2026). RSpectra: Solvers for Large-Scale Eigenvalue and SVD Problems. R package. — Efficient partial eigenvalue and singular value decompositions using Krylov methods. https://CRAN.R-project.org/package=RSpectra.