Solving Ax = b
The central problem of linear algebra: finding the vector x that a matrix A maps to a given vector b.
The three outcomes
A linear system Ax = b has exactly one of three fates: a unique solution, no solution, or infinitely many. Which one occurs depends on the rank of A and whether b lies in the column space. A square full-rank matrix always gives a unique solution; otherwise the system may be inconsistent or underdetermined.
Existence and uniqueness
- a solution exists when b is in the column space of A, i.e. rank(A) = rank([A | b])
- the solution is unique when the columns of A are independent (full column rank)
- if a solution exists and the null space is nontrivial, there are infinitely many
- the general solution is any particular solution plus any null-space vector
How to solve in practice
Never invert the matrix. For a general square system, factor A = LU with partial pivoting and run two triangular solves. For symmetric positive-definite systems, use Cholesky. For overdetermined systems with no exact solution, use QR to find the least-squares answer. For very large sparse systems, use iterative methods. The choice is dictated by the size, structure, and conditioning of A.
Conditioning governs accuracy
Even a correctly implemented solver returns an inaccurate answer if the matrix is ill conditioned, because the problem itself amplifies data error. Always consider the condition number when interpreting a computed solution, and regularize when the matrix is nearly singular.
import numpy as np
A = np.array([[3.0, 2.0], [1.0, 2.0]])
b = np.array([7.0, 5.0])
x = np.linalg.solve(A, b) # uses LU with pivoting internally
print(x, np.allclose(A @ x, b))Every timestep of an implicit physics simulation reduces to solving a large linear system for the updated field, so the efficiency of the Ax = b solve sets the pace of the entire model run.