14.4: Overdetermined Systems

Solving Overdetermined Systems of Linear Equations

An overdetermined system of linear equations is a system where there are more equations than variables. This means \(m \gt n\) in the equation \(Ax=b\), where A is an \(m \times n\) matrix.

Characteristics of Overdetermined Systems

• More Equations than Variables: The defining characteristic is that you have more constraints (equations) than unknowns (variables).
• No Exact Solution (Generally): In most real-world scenarios, an overdetermined system will not have an exact solution that satisfies all equations simultaneously. This is because the equations often represent measurements or observations that contain some level of error or inconsistency.
• Inconsistency: The lines (or planes, hyperplanes) represented by the equations typically do not intersect at a single common point.

The Least Squares Solution

Since an exact solution is usually not possible, the goal for overdetermined systems is to find the "best approximate" solution. This is where the least squares method comes in. The least squares solution is the vector x that minimizes the sum of the squares of the residuals. A residual is the difference between the observed value (from b) and the value predicted by the model (\(Ax\)).

Mathematically, minimize the Euclidean norm of the residual vector: \(\min_x​ \parallel Ax − b \parallel^2\)

This minimization problem has a unique solution (if A has full column rank), given by the normal equations: \(A^TAx = A^Tb\)

The solution x can then be found by: \(x = (A^T A)^{−1} A^T b\)

The term \((A^T A)^{−1} A^T\) is also known as the pseudoinverse of A (specifically, the left pseudoinverse, \(A^{+}\)), which is particularly useful when A has full column rank. So, \(x = A^{+} b\).

Why Least Squares?

The least squares approach is widely used because:

• It provides a single, unique solution even when an exact solution doesn't exist.
• It is computationally efficient.
• It has strong statistical foundations (e.g., it corresponds to the maximum likelihood estimator under assumptions of normally distributed errors).

Solving in Python with NumPy

NumPy's numpy.linalg.lstsq function is the standard and most robust way to solve overdetermined systems in Python. It uses a numerically stable method (Singular Value Decomposition - SVD) to compute the least squares solution.

The lstsq function returns:

1. \(x\): The least-squares solution.
2. residuals: The sum of squared residuals, or an empty array if the rank of A is less than the number of variables or if b is one dimensional.
3. rank: The effective rank of matrix A.
4. s: The singular values of A.

The \(x\) returned by lstsq is the vector that minimizes \(\parallel Ax − b\parallel ^2\). While \(Ax\) will not be exactly equal to b, it will be the closest possible approximation in a least-squares sense.

Example \(\PageIndex{1}\)

Overdetermined System 

import numpy as np

# Define an example overdetermined system:
# x + y = 2
# 2x - y = 1
# x + 2y = 4
A_over = np.array([[1, 1],
                   [2, -1],
                   [1, 2]])

b_over = np.array([2, 1, 4])

print("\nMatrix A (Overdetermined):\n", A_over)
print("\nVector b (Overdetermined):\n", b_over)

# --- Solving Overdetermined System using np.linalg.lstsq ---
# For overdetermined systems, lstsq finds the least-squares solution.
print("\n--- Solving Overdetermined System using np.linalg.lstsq ---")
x_lstsq_over, residuals_over, rank_over, s_over = np.linalg.lstsq(A_over, b_over, rcond=None)

print("\nLeast-squares solution from np.linalg.lstsq (Overdetermined):\n", x_lstsq_over)
print("Residuals (sum of squared residuals):", residuals_over)
print("Rank of A_over from lstsq:\n", rank_over)
print("Singular values of A_over:\n", s_over)

# Verify the solution (A @ x should be close to b, but not necessarily equal)
print("\nVerification (A_over @ x_lstsq_over):\n", A_over @ x_lstsq_over)
print("\nOriginal b_over:\n", b_over)
print("\nDifference (A_over @ x_lstsq_over - b_over):\n", A_over @ x_lstsq_over - b_over)

# The residuals array returned by lstsq is the sum of squared residuals,
# or an empty array if rank(A) < n or if b is 1-D.
# To get individual residuals:
individual_residuals_over = A_over @ x_lstsq_over - b_over
print("\nIndividual residuals (A_over @ x - b_over):\n", individual_residuals_over)
print("Sum of squared individual residuals:", np.sum(individual_residuals_over**2))

Matrix A (Overdetermined):
 [[ 1  1]
 [ 2 -1]
 [ 1  2]]

Vector b (Overdetermined):
 [2 1 4]

--- Solving Overdetermined System using np.linalg.lstsq ---

Least-squares solution from np.linalg.lstsq (Overdetermined):
 [1.11428571 1.31428571]
Residuals (sum of squared residuals): [0.25714286]
Rank of A_over from lstsq:
 2
Singular values of A_over:
 [2.64575131 2.23606798]

Verification (A_over @ x_lstsq_over):
 [2.42857143 0.91428571 3.74285714]

Original b_over:
 [2 1 4]

Difference (A_over @ x_lstsq_over - b_over):
 [ 0.42857143 -0.08571429 -0.25714286]

Individual residuals (A_over @ x - b_over):
 [ 0.42857143 -0.08571429 -0.25714286]
Sum of squared individual residuals: 0.25714285714285734