Solving Linear Systems

A system of linear equations can be written in matrix form as Ax = b, where A is the coefficient matrix, x is the vector of unknowns, and b is the right-hand side vector. For example, two equations with two unknowns can be expressed as a 2×2 matrix equation. Solving the system means finding x such that A multiplied by x equals b. NumPy's np.linalg.solve(A, b) handles this efficiently and accurately.

Mathematically, the solution to Ax = b is x = A⁻¹b. But computing np.linalg.inv(A) @ b is slower and less accurate than calling np.linalg.solve(A, b) directly. solve() uses LU decomposition internally, which requires fewer operations and accumulates less floating-point error. The rule of thumb: if you are solving for x in Ax=b, always use solve, never inv.

np.linalg.solve(A, b) requires that A is square and non-singular. If A is singular (determinant is zero), the function raises np.linalg.LinAlgError: Singular matrix. A is singular when rows are linearly dependent — for example, one equation is just a multiple of another. Always check that your system is well-posed before calling solve. If A is not square, use np.linalg.lstsq() instead.

np.linalg.solve internally applies LU decomposition: it factors A into a lower-triangular matrix L and an upper-triangular matrix U such that A = LU. Solving then becomes two simpler triangular solves (forward and back substitution), which are far more efficient than Gaussian elimination repeated from scratch. This is the same algorithm used by MATLAB, Fortran LAPACK, and virtually every scientific computing library.

np.linalg.solve(A, B) can handle a matrix B with multiple columns, solving for each column simultaneously. If B has shape (n, k), the result X has shape (n, k), where each column of X is the solution for the corresponding column of B. This is useful in physics simulations and machine learning where you need to solve the same coefficient matrix with many different right-hand sides.

When the system is over-determined (more equations than unknowns), no exact solution exists. np.linalg.lstsq(A, b) finds the x that minimises the sum of squared residuals ||Ax - b||² — this is exactly what linear regression computes. It returns the solution x, residuals, rank, and singular values. The rcond=None argument suppresses a deprecation warning and sets the cutoff for treating singular values as zero.

The normal equations of ordinary least squares are A.T @ A @ x = A.T @ b, which is a square system solvable with np.linalg.solve. This gives the same result as lstsq when A.T @ A is well-conditioned. In practice, lstsq is preferred because it uses SVD and is stable even when columns are nearly collinear, but the normal equations help you understand the math behind regression.

Before solving, it is good practice to check the condition number with np.linalg.cond(A). A condition number close to 1 means the system is well-conditioned and the solution is reliable. A very large condition number (e.g., 1e12) means the matrix is nearly singular: small changes in b produce huge changes in x, and the solution may be numerically unreliable. In that case, consider regularisation (ridge regression) or using a more robust solver.

For very large systems where A has mostly zero entries (a sparse matrix), storing the full matrix wastes memory and using dense np.linalg.solve is slow. SciPy provides scipy.sparse and scipy.sparse.linalg.spsolve that exploit sparsity for huge systems like finite element models or network problems. Understanding when your matrix is sparse is an important performance consideration in scientific computing.

Once you have the solution vector x from np.linalg.solve(A, b), always verify it by computing A @ x and checking it against b using np.allclose(A @ x, b). In floating-point arithmetic, you rarely get exact equality — allclose checks within a small tolerance. The residual ||Ax - b|| should be near machine epsilon times the scale of the problem for a well-conditioned system.

Linear solvers appear throughout machine learning: linear regression solves the normal equations, Gaussian processes solve kernel matrix systems, and Kalman filters solve state-update equations. By understanding np.linalg.solve and lstsq, you can implement these algorithms from scratch and debug them when results are unexpected. NumPy's linalg module gives you the same numerical tools used in production scientific libraries.

Test your understanding of Data Analysis concepts from this lesson.

In this lesson you learned: np.linalg.solve(A, b) efficiently solves square linear systems using LU decomposition, np.linalg.lstsq() finds the least-squares solution for over-determined systems, and the condition number measures numerical stability — large values warn of near-singular matrices. Next up we explore eigenvalues and the SVD, which underpin PCA and many ML algorithms.