Solving Linear Systems and the Matrix Inverse

The General Linear System

A system of linear equations with $m$ equations and $n$ unknowns can be written compactly as:

Here $\mathbf{A}$ is an $m \times n$ matrix of coefficients, $\mathbf{x}$ is the $n$-dimensional unknown vector, and $\mathbf{b}$ is the $m$-dimensional right-hand side vector.

The shape of $\mathbf{A}$ determines what kind of solution, if any, exists:

• Square ($m = n$, invertible): a unique solution exists, $\mathbf{x} = \mathbf{A}^{-1}\mathbf{b}$.
• Overdetermined ($m > n$): more equations than unknowns, typically no exact solution. Use least squares.
• Underdetermined ($m < n$): fewer equations than unknowns, infinitely many solutions. Choose the one with minimum norm.

When the System Is Square and Invertible

If $\mathbf{A}$ is square and all its columns are linearly independent (full rank), then $\mathbf{A}^{-1}$ exists and the solution is unique.

The classic Gaussian elimination algorithm solves such systems by performing row operations that transform $\mathbf{A}$ into an upper triangular matrix, then performing back-substitution. This is what numpy.linalg.solve does under the hood, and it is the method you should use in practice:

import numpy as np


#  3x +  y +  z = 12
#   x + 2y + 3z = 14
#   x +  y -  z =  2

A = np.array([[3, 1, 1],
              [1, 2, 3],
              [1, 1, -1]])

b = np.array([12, 14, 2])

x = np.linalg.solve(A, b)
print(f"Solution: x={x[0]:.2f}, y={x[1]:.2f}, z={x[2]:.2f}")

# Verify
print("Verification A @ x:", A @ x)
print("Matches b:", np.allclose(A @ x, b))

PROMPT

Use numpy to solve a 3x3 linear system with np.linalg.solve and verify the solution with np.allclose.

LU Decomposition: Breaking A into Pieces

For repeated solves with the same $\mathbf{A}$ but different $\mathbf{b}$ values, it is efficient to decompose $\mathbf{A}$ once:

where $\mathbf{L}$ is lower triangular and $\mathbf{U}$ is upper triangular. Once you have $\mathbf{L}$ and $\mathbf{U}$, solving $\mathbf{A}\mathbf{x} = \mathbf{b}$ requires two forward/back substitutions instead of full Gaussian elimination. This matters in portfolio optimization, where you solve the same system structure repeatedly with different target returns.

SciPy provides this:

import scipy.linalg as la

A = np.array([[3, 1, 1],
              [1, 2, 3],
              [1, 1, -1]])

P, L, U = la.lu(A)
print("L:\n", L)
print("U:\n", U)
print("P @ L @ U:\n", P @ L @ U)

PROMPT

Use scipy.linalg to do LU decomposition of a 3x3 matrix and verify that P @ L @ U reconstructs the original.

Overdetermined Systems: Least Squares

When you have more equations than unknowns, as in regression with many observations, an exact solution generally does not exist. The goal shifts to finding the $\mathbf{x}$ that makes $\mathbf{A}\mathbf{x}$ as close as possible to $\mathbf{b}$, in the least-squares sense:

The solution to this problem is given by the normal equations:

If $\mathbf{A}^\top\mathbf{A}$ is invertible, the solution is:

However, computing $(\mathbf{A}^\top\mathbf{A})^{-1}$ directly is numerically unstable, rounding errors compound. The QR decomposition is the preferred approach:

where $\mathbf{Q}$ is orthogonal ($\mathbf{Q}^\top\mathbf{Q} = \mathbf{I}$) and $\mathbf{R}$ is upper triangular. The solution is then $\mathbf{x} = \mathbf{R}^{-1}\mathbf{Q}^\top\mathbf{b}$, computed via back-substitution, which is far more numerically stable.

import numpy as np

# Overdetermined: 10 observations, 2 unknowns
np.random.seed(42)
x_obs = np.linspace(0, 5, 10)
y_obs = 2.5 * x_obs + 1.5 + np.random.randn(10) * 0.5

# Build A: [x, 1] for each observation
A = np.column_stack([x_obs, np.ones(10)])
b = y_obs

# Solve least squares
x_ls, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)
print(f"Least squares solution: slope={x_ls[0]:.3f}, intercept={x_ls[1]:.3f}")

# Compare with np.polyfit
slope, intercept = np.polyfit(x_obs, y_obs, 1)
print(f"numpy.polyfit: slope={slope:.3f}, intercept={intercept:.3f}")

PROMPT

Use numpy to fit a noisy linear relationship with np.linalg.lstsq and compare the result to np.polyfit.

Underdetermined Systems: Minimum Norm

When there are fewer equations than unknowns, infinitely many solutions exist. The minimum norm solution selects the solution with the smallest Euclidean norm $\|\mathbf{x}\|$. For the underdetermined case with full row rank, this is:

In practice, numpy.linalg.lstsq handles all three cases, square, overdetermined, and underdetermined, returning the minimum norm solution when applicable.

Financial Applications

Linear regression is a least-squares problem. In the regression article, you fit $Y = \beta_0 + \beta_1 X + \varepsilon$ by solving the normal equations. The $x$-values are the observations of the independent variable, and the $y$-values are the observations of the dependent variable. The solution gives you the regression coefficients.

Portfolio optimization constraints often produce linear systems. When you set up the problem of finding weights that achieve a target return with minimum variance, the first-order conditions produce a linear system in the weights and the Lagrange multiplier for the return constraint.

Option pricing in the binomial model requires solving a linear system that propagates prices backward through the tree.

import numpy as np

# Simulate a linear regression problem explicitly
np.random.seed(0)
n = 100
X = np.random.randn(n, 3)
y = X @ [1.5, -0.8, 2.1] + np.random.randn(n) * 0.5

# Build A = X, b = y, solve via normal equations
A = X
b = y

# Using lstsq (stable)
beta, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)

# Compare to the known true coefficients [1.5, -0.8, 2.1]
print(f"Estimated: {beta}")
print(f"True:      [1.5, -0.8, 2.1]")

PROMPT

Use numpy to fit a 3-feature linear regression with np.linalg.lstsq on simulated data and print the estimated coefficients.

References

Strang, Introduction to Linear Algebra, 5th edition, Wellesley-Cambridge Press, Chapters 2–3. ISBN: 978-0980232776.

Golub & Van Loan, Matrix Computations, 4th edition, Johns Hopkins University Press, Chapter 3. ISBN: 978-1421407944.