14.5: A General Solution

Introduction to A General Solution

When dealing with systems of linear equations, we often encounter scenarios where the standard matrix inverse isn't applicable—either because the matrix is not square (meaning we have more equations than unknowns, or vice versa) or because it's singular. In such cases, finding a unique solution can be challenging, or even impossible in the traditional sense.

This document introduces the Moore-Penrose pseudo-inverse, a powerful generalization of the matrix inverse that allows us to find a "best-fit" or general solution to any system of linear equations, regardless of whether the matrix is square or invertible. It provides a robust method for tackling problems that are either overdetermined (more equations than unknowns) or underdetermined (fewer equations than unknowns), offering a unique solution that minimizes the error.

Let's explore how the Moore-Penrose pseudo-inverse provides a comprehensive approach to solving linear systems.

In NumPy, the Moore-Penrose pseudo-inverse of a matrix is computed using the numpy.linalg.pinv() function. This function provides a generalized inverse that exists for any matrix, including rectangular and singular matrices where a traditional inverse does not exist.

Key aspects of numpy.linalg.pinv():

• Generalized Inverse:
• Singular Value Decomposition (SVD):
• Handling Singular Values:
• Applications: The pseudo-inverse is widely used in various fields, including:
  • Least-squares regression: Finding the "best approximate solution" to over-determined or under-determined systems of linear equations.
  • Machine learning: In algorithms like linear regression or support vector machines.
  • Signal processing: For inverse problems and system identification.

Moore-Penrose Pseudo-Inverse

The Moore-Penrose pseudo-inverse (often denoted as \(A^{+}\)) is a generalization of the inverse matrix. It allows us to find a "best-fit" solution to a system of linear equations \(Ax = b\), even when the matrix A is not square or is singular (non-invertible).

For a system \(Ax = b\), the pseudo-inverse solution is given by \(x = A^{+}b\).

This solution has specific properties:

1. If a unique solution exists, the pseudo-inverse finds it.
2. If the system is overdetermined (more equations than unknowns, A is tall), the pseudo-inverse finds the least-squares solution that minimizes \(\parallel Ax - b \parallel^2\).
3. If the system is underdetermined (fewer equations than unknowns, A is wide), the pseudo-inverse finds the solution with the minimum Euclidean norm (\(\parallel x \parallel^2\)) among all possible solutions.

Examples:

Example \(\PageIndex{1}\)

Well-determined system with square and invertible A matrix

 In this case, the pseudo-inverse will behave like the regular inverse.

import numpy as np
print("--- Example 1: Well-determined System (Unique Solution) ---")
A1 = np.array([[2, 1],
               [1, 3]])
b1 = np.array([4, 7])

print("Matrix A1:\n", A1)
print("Vector b1:\n", b1)

# Calculate the pseudo-inverse of A1
A1_plus = np.linalg.pinv(A1)
print("Pseudo-inverse of A1:\n", A1_plus)

# Solve for x1
x1 = A1_plus @ b1
print("Solution x1 (A1+ @ b1):\n", x1)
print("Verification (A1 @ x1):\n", A1 @ x1)
print("-" * 50 + "\n")

--- Example 1: Well-determined System (Unique Solution) ---
Matrix A1:
 [[2 1]
 [1 3]]
Vector b1:
 [4 7]
Pseudo-inverse of A1:
 [[ 0.6 -0.2]
 [-0.2  0.4]]
Solution x1 (A1+ @ b1):
 [1. 2.]
Verification (A1 @ x1):
 [4. 7.]
--------------------------------------------------

Example \(\PageIndex{2}\)

Overdetermined System (more equations than unknowns)

 No exact solution exists, so the pseudo-inverse finds the least-squares solution.

import numpy as np
import matplotlib.pyplot as plt
print("--- Example 2: Overdetermined System (Least-Squares Solution) ---")
A2 = np.array([[1, 1],
               [2, 1],
               [3, 1]])
b2 = np.array([2, 3, 5]) # This system is inconsistent

print("Matrix A2:\n", A2)
print("Vector b2:\n", b2)

# Calculate the pseudo-inverse of A2
A2_plus = np.linalg.pinv(A2)
print("Pseudo-inverse of A2:\n", A2_plus)

# Solve for x2 (least-squares solution)
x2 = A2_plus @ b2
print("Solution x2 (A2+ @ b2 - Least-Squares):\n", x2)
print("Verification (A2 @ x2) - Note the error:\n", A2 @ x2)
print("Residual (||A2 @ x2 - b2||):\n", np.linalg.norm(A2 @ x2 - b2))
print("-" * 50 + "\n")
x = np.linspace(1,2,200)
y1 = -A2[0,0]*x + b2[0]
y2 = -A2[1,0]*x + b2[1]
y3 = -A2[2,0]*x + b2[2]
plt.figure(figsize=(8, 6))
plt.plot(x, y1, label='Line 1: $x_1 + x_2 = 2$')
plt.plot(x, y2, label='Line 2: $2 x_1 + x_2 =  3$')
plt.plot(x, y3, label='Line 3: $3 x_1 + x_2 =  5$')
plt.scatter(x2[0,],x2[1,], color='red', label='Least Squares Solution')
plt.legend()
plt.title('Overdetermined System Visualization')
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.grid()
plt.show()

--- Example 2: Overdetermined System (Least-Squares Solution) ---
Matrix A2:
 [[1 1]
 [2 1]
 [3 1]]
Vector b2:
 [2 3 5]
Pseudo-inverse of A2:
 [[-5.00000000e-01 -2.78472776e-16  5.00000000e-01]
 [ 1.33333333e+00  3.33333333e-01 -6.66666667e-01]]
Solution x2 (A2+ @ b2 - Least-Squares):
 [1.5        0.33333333]
Verification (A2 @ x2) - Note the error:
 [1.83333333 3.33333333 4.83333333]
Residual (||A2 @ x2 - b2||):
 0.408248290463863
--------------------------------------------------

Overdetermined system with three equations and two unknowns shown with least squares solution.

Example \(\PageIndex{3}\)

Underdetermined System (fewer equations than unknowns)

Many solutions exist, the pseudo-inverse finds the one with minimum norm.

import numpy as np
print("--- Example 3: Underdetermined System (Minimum-Norm Solution) ---")
A3 = np.array([[1, 2, 3],
               [4, 5, 6]])
b3 = np.array([7, 8])

print("Matrix A3:\n", A3)
print("Vector b3:\n", b3)

# Calculate the pseudo-inverse of A3
A3_plus = np.linalg.pinv(A3)
print("Pseudo-inverse of A3:\n", A3_plus)

# Solve for x3 (minimum-norm solution)
x3 = A3_plus @ b3
print("Solution x3 (A3+ @ b3 - Minimum-Norm):\n", x3)
print("Verification (A3 @ x3):\n", A3 @ x3)
print("Norm of x3 (||x3||):\n", np.linalg.norm(x3))
print("-" * 50 + "\n")

--- Example 3: Underdetermined System (Minimum-Norm Solution) ---
Matrix A3:
 [[1 2 3]
 [4 5 6]]
Vector b3:
 [7 8]
Pseudo-inverse of A3:
 [[-0.94444444  0.44444444]
 [-0.11111111  0.11111111]
 [ 0.72222222 -0.22222222]]
Solution x3 (A3+ @ b3 - Minimum-Norm):
 [-3.05555556  0.11111111  3.27777778]
Verification (A3 @ x3):
 [7. 8.]
Norm of x3 (||x3||):
 4.482476167543176
--------------------------------------------------

Example \(\PageIndex{4}\)

Singular (non-invertible square) Matrix

Even though a square and singular matrix system, does not have a regular inverse, the pseudo-inverse can still find a solution (least-squares if inconsistent, minimum-norm if consistent with multiple solutions).

import numpy as np
print("--- Example 4: Singular Square Matrix ---")
A4 = np.array([[1, 2],
               [2, 4]]) # This matrix is singular (rows are linearly dependent)
b4 = np.array([3, 6]) # This system is consistent (multiple solutions)

print("Matrix A4:\n", A4)
print("Vector b4:\n", b4)

# Calculate the pseudo-inverse of A4
A4_plus = np.linalg.pinv(A4)
print("Pseudo-inverse of A4:\n", A4_plus)

# Solve for x4
x4 = A4_plus @ b4
print("Solution x4 (A4+ @ b4 - Minimum-Norm):\n", x4)
print("Verification (A4 @ x4):\n", A4 @ x4)
print("Norm of x4 (||x4||):\n", np.linalg.norm(x4))
print("-" * 50 + "\n")

--- Example 4: Singular Square Matrix ---
Matrix A4:
 [[1 2]
 [2 4]]
Vector b4:
 [3 6]
Pseudo-inverse of A4:
 [[0.04 0.08]
 [0.08 0.16]]
Solution x4 (A4+ @ b4 - Minimum-Norm):
 [0.6 1.2]
Verification (A4 @ x4):
 [3. 6.]
Norm of x4 (||x4||):
 1.3416407864998734
--------------------------------------------------

Interpretation of Results:

1. For well-determined systems, x will be the exact unique solution.
2. For overdetermined systems, x will be the least-squares approximation. The residual (\(\parallel Ax - b \parallel\)) will indicate how well the solution fits.
3. For underdetermined systems, x will be one of infinitely many solutions, specifically the one with the smallest Euclidean norm.