Principle:Haifengl Smile Linear System Solving

Overview

Linear System Solving is the process of finding a vector $x$ that satisfies the equation $Ax=b$, where $A$ is a known matrix and $b$ is a known vector. This is one of the most fundamental problems in scientific computing, appearing in regression, optimization, physics simulation, signal processing, and machine learning.

Smile provides two classes of solvers:

• Direct solvers that first decompose $A$ (LU, QR, Cholesky, SVD) and then solve via substitution. These are exact (up to floating-point precision) and have predictable cost.
• Iterative solvers (Biconjugate Gradient) that approximate the solution through successive refinement. These are efficient for large sparse systems where direct decomposition is too expensive.

Theoretical Basis

Direct Solving via Decomposition

Direct methods transform the original system into a sequence of easily solvable triangular systems.

LU Solve

Given $PA=LU$, solving $Ax=b$ proceeds in three steps:

1. Permute: $Pb$ (reorder $b$ according to the pivot vector)
2. Forward substitution: Solve $Ly=Pb$ where $L$ is lower triangular
3. Back substitution: Solve $Ux=y$ where $U$ is upper triangular

$\begin{array}{rl}{y}_i& =(Pb{)}_i-\sum _{j=1}^{i-1}{l}_{ij}{y}_j,i=1,2,\dots ,n\\ x_i& =\frac{y_i-\sum _{j=i+1}^n{u}_{ij}{x}_j}{u_{ii}},i=n,n-1,\dots ,1\end{array}$

Complexity: $O(n^2)$ for the solve phase (after $O(n^3)$ factorization).

LAPACK routine: dgetrs / sgetrs

Cholesky Solve

Given $A=L{L}^T$ (SPD matrix), solving $Ax=b$:

1. Forward substitution: Solve $Ly=b$
2. Back substitution: Solve $L^{T}x=y$

Complexity: $O(n^2)$ for the solve, $O(\frac{1}{3}{n}^3)$ for factorization.

Advantage over LU: Approximately 2x faster because Cholesky factors are half the size and no pivoting is required.

LAPACK routine: dpotrs / spotrs

QR Solve (Least Squares)

For overdetermined systems ($m>n$), QR solves the least squares problem:

$\underset{x}{\min }\Vert Ax-b{\Vert }_2$

Given $A=QR$:

1. Apply $Q^T$: Compute $Q^{T}b$ using dormqr/sormqr
2. Triangular solve: Solve $Rx=(Q^{T}b{)}_{1:n}$ using dtrtrs/strtrs

The solution minimizes the residual $\Vert Ax-b{\Vert }_2$ and is numerically more stable than the normal equations approach $A^{T}Ax=A^{T}b$.

LAPACK routines: dormqr + dtrtrs / sormqr + strtrs

SVD Solve (Pseudoinverse)

For rank-deficient or ill-conditioned systems, SVD provides the most robust solution:

$x=V{\mathrm{\Sigma }}^{-1}{U}^{T}b$

where ${\mathrm{\Sigma }}^{-1}$ is formed by inverting only the nonzero (above threshold) singular values. Small singular values (${\sigma }_i\le \text{rcond}$) are treated as zero, providing a regularized solution.

1. Compute $U^{T}b$ (project $b$ onto left singular vectors)
2. Divide by singular values: $({\mathrm{\Sigma }}^{-1}{U}^{T}b{)}_i=(U^{T}b{)}_i/{\sigma }_i$
3. Compute $V({\mathrm{\Sigma }}^{-1}{U}^{T}b)$ (expand in right singular vector basis)

Advantages: Works for any matrix (including singular and rank-deficient); provides the minimum-norm solution when the system is underdetermined.

Iterative Solving

Biconjugate Gradient Method (BiCG)

For large sparse systems where direct decomposition is impractical, the Biconjugate Gradient method iteratively refines an approximate solution. Given $Ax=b$:

1. Start with initial guess $x_0$
2. Compute residual $r_0=b-A{x}_0$
3. Iterate: update $x_k$ using search directions $p_k$ derived from the residual

Key operations per iteration:

• One matrix-vector product $Ap$
• One transpose matrix-vector product $A^{T}pp$
• One preconditioner solve

Convergence criteria (controlled by itol parameter):

itol	Criterion	Description
1	$\Vert Ax-b\Vert /\Vert b\Vert <\text{tol}$	Relative residual norm
2	$\Vert A^{-1}(Ax-b)\Vert /\Vert A^{-1}b\Vert <\text{tol}$	Preconditioned residual
3	$\Vert x_{k+1}-x_k{\Vert }_2<\text{tol}$	Solution change (L2)
4	$\Vert x_{k+1}-x_k{\Vert }_{\mathrm{\infty }}<\text{tol}$	Solution change (Linf)

Preconditioning: The default preconditioner is Jacobi (diagonal scaling), which divides each equation by its diagonal element. Better preconditioners (ILU, SSOR) can be supplied via the Preconditioner interface.

Choosing the Right Solver

Solver	Matrix Type	System Type	When to Use
LU	Square, non-singular	$Ax=b$ (exact)	General-purpose; multiple right-hand sides
Cholesky	SPD	$Ax=b$ (exact)	Covariance systems, optimization; 2x faster than LU
QR	Overdetermined ($m>n$)	Least squares	Regression, fitting; more stable than normal equations
SVD	Any (including rank-deficient)	Least squares / min-norm	Ill-conditioned systems; need rank information
BiCG	Large, sparse	$Ax=b$ (iterative)	When decomposition is too expensive; only needs $Av$ and $A^{T}v$

Multiple Right-Hand Sides

All direct solvers in Smile support solving $AX=B$ where $B$ is a matrix (multiple right-hand sides simultaneously). The matrix $B$ is overwritten in-place with the solution $X$.

This is more efficient than solving individual systems because:

• The decomposition is computed once
• LAPACK can optimize memory access patterns for the block solve

Numerical Stability Considerations

• LU with pivoting: Partial pivoting ensures stability for most practical matrices, but can amplify errors for pathological cases.
• QR: Always backward stable; preferred when numerical accuracy is paramount.
• Cholesky: The most stable option for SPD matrices because it exploits positive definiteness.
• SVD: The most robust option overall; can handle singular and near-singular matrices gracefully through thresholding.
• BiCG: Convergence depends on the condition number $\kappa (A)$ and the quality of the preconditioner. May stagnate for highly ill-conditioned systems.

Relationship to the Pipeline

Construction --> Arithmetic --> Decomposition --> Solving --> Result Extraction
                                                    ^
                                                    |
                                             (this principle)

Linear system solving is the fourth stage of the pipeline. It consumes decomposition results (LU, QR, Cholesky, SVD records) produced by the decomposition stage and produces solution vectors or matrices.

Related

• Implementation:Haifengl_Smile_Decomposition_Solvers
• Principle:Haifengl_Smile_Matrix_Decomposition

Knowledge Sources

• Smile
• LAPACK Users' Guide
• Templates for Linear Systems

Domains

Linear_Algebra, Numerical_Computing, Optimization