14.3: Underdetermined Systems

Solving Underdetermined Systems of Linear Equations

An underdetermined system of linear equations is a system where there are more variables than equations. Such a system typically has infinitely many solutions, or no solutions at all (though the latter is less common when discussing "solutions"). When solutions exist, they can often be expressed in terms of free variables.

Consider a linear system in the form \(Ax=b\), where:

• A is an \(m \times n\) matrix, with \(m \lt n\) , fewer equations than variables.
• x is an \(n \times 1\) column vector of variables.
• b is an \(m \times1\) column vector of constants.

Since there are more variables than equations, the system has at least n−m free variables. This means we can express some variables in terms of others.

General Approach

1. Identify the rank of the matrix A: The rank of matrix A determines the number of linearly independent rows (or columns). For a system to have solutions, the rank of the augmented matrix \( \left [A \vert b \right ]  \) must be equal to the rank of A.
2. Find a particular solution: One way to find a particular solution is to use the pseudoinverse of A. The pseudoinverse, denoted \(A^{+}\), can be used to find the minimum-norm solution (the solution with the smallest Euclidean norm) when infinite solutions exist. This solution is given by \(x_p ​= A^{+} b\).
3. Find the null space of A: The null space of A consists of all vectors \(x_n\)​ such that \(A x_n​ = 0\). Any vector in the null space can be added to a particular solution to form another valid solution.
4. Construct the general solution: The general solution to an underdetermined system is the sum of a particular solution and any vector from the null space: \(x = x_p​+x_h​\), where \(x_h ​\in null(A)\).

In Python, libraries like the NumPy and SciPy provide powerful linear algebra tools to handle the pseudoinverse operation and for finding the null space.

NumPy also provides the numpy.linalg.lstsq function to compute a least-squares solution to a linear matrix equation. While primarily designed for overdetermined systems (where it finds the "best fit" solution), it can also be used for underdetermined systems. For underdetermined systems, lstsq returns the minimum-norm solution (similar to using the pseudoinverse).

When lstsq is used for an underdetermined system, it returns the solution that has the smallest Euclidean norm (i.e., the "shortest" vector \(x\) that satisfies the equations). This is often the particular solution \(x_p\)​.

Example \(\PageIndex{1}\)

Underdetermined System

Solution

Let's consider a simple underdetermined system:

\[\begin{align}
\begin{split}
x + 2 y - z & = 5 \\
3 x - y + 2 z & = 1
\end{split}
\end{align}\]

Here, we have 2 equations and 3 variables \( \left ( x, y, z \right ) \). The matrix A would be:

\begin{bmatrix}
1 & 2 & -1 \\
3 & -1 & 2
\end{bmatrix}

The vector b would be:

\begin{bmatrix}
5 \\
1
\end{bmatrix}

We can use NumPy's np.linalg.pinv to find the pseudoinverse and then a particular solution where scipy.linalg.null_space would provide the null space.

Example \(\PageIndex{2}\)

Underdetermined System

Python Example 

Compare the results from the pseudoinvers and the least-squares (lstsq)

# --- Example: Underdetermined System ---
import numpy as np
from scipy.linalg import null_space

# Define the coefficient matrix A and the constant vector b
A = np.array([[1, 2, -1],
              [3, -1, 2]])

b = np.array([5, 1])

print("Matrix A:\n", A)
print("\nVector b:\n", b)

# --- Step 1: Check for consistency (optional, but good practice) ---
# For an underdetermined system, if solutions exist, they are infinite.
# The rank of A must be equal to the rank of the augmented matrix [A|b]
# for a solution to exist.
rank_A = np.linalg.matrix_rank(A)
augmented_A = np.hstack((A, b.reshape(-1, 1))) # Create augmented matrix
rank_augmented_A = np.linalg.matrix_rank(augmented_A)

print(f"\nRank of A: {rank_A}")
print(f"Rank of augmented A: {rank_augmented_A}")

if rank_A != rank_augmented_A:
    print("\nNo solution exists for this system.")
else:
    print("\nSolutions exist for this system (potentially infinite).")

    # --- Solution Method 1: Using np.linalg.pinv (Pseudoinverse) ---
    # The pseudoinverse (Moore-Penrose inverse) gives the minimum-norm solution.
    print("\n--- Solving using np.linalg.pinv (Pseudoinverse) ---")
    A_pinv = np.linalg.pinv(A)
    x_particular_pinv = A_pinv @ b

    print("\nPseudoinverse of A:\n", A_pinv)
    print("\nParticular solution (minimum-norm) from pinv:\n", x_particular_pinv)

    # Verify the particular solution
    print("\nVerification (A @ x_particular_pinv):\n", A @ x_particular_pinv)

    # --- Finding the null space of A (applicable to both methods for general solution) ---
    # The null space contains all vectors x_n such that A @ x_n = 0.
    # Any linear combination of these basis vectors can be added to the particular solution.
    # scipy.linalg.null_space returns an orthonormal basis for the null space.
    null_space_basis = null_space(A)

    print("\nBasis for the null space of A:\n", null_space_basis)

    # --- Construct the general solution ---
    # The general solution is x = x_particular + c1*v1 + c2*v2 + ...
    # where v1, v2, ... are the basis vectors of the null space, and c1, c2, ... are arbitrary scalars.

    print("\nGeneral Solution Form (using pinv's particular solution):")
    print("x = x_particular_pinv + c * null_space_basis_vector")
    # Assuming one basis vector for simplicity, adjust if null_space_basis has more columns
    if null_space_basis.shape[1] > 0:
        print(f"x = {x_particular_pinv} + c * {null_space_basis[:, 0]}")
    else:
        print("Null space is trivial (only the zero vector).")


    # Example: Let's pick an arbitrary scalar 'c' (e.g., c = 10)
    if null_space_basis.shape[1] > 0:
        c = 10
        another_solution = x_particular_pinv + c * null_space_basis[:, 0]
        print(f"\nExample solution with c = {c}:\n", another_solution)

        # Verify this new solution
        print("\nVerification (A @ another_solution):\n", A @ another_solution)

        # Example: Let's pick an arbitrary scalar 'c' (e.g., c = -5)
        c = -5
        yet_another_solution = x_particular_pinv + c * null_space_basis[:, 0]
        print(f"\nExample solution with c = {c}:\n", yet_another_solution)

        # Verify this new solution
        print("\nVerification (A @ yet_another_solution):\n", A @ yet_another_solution)


    # --- Solution Method 2: Using np.linalg.lstsq ---
    # For underdetermined systems, lstsq returns the minimum-norm solution.
    print("\n--- Solving using np.linalg.lstsq ---")
    x_lstsq, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)

    print("\nSolution from np.linalg.lstsq (minimum-norm):\n", x_lstsq)
    print("Residuals (should be close to zero for consistent systems):\n", residuals)
    print("Rank of A from lstsq:\n", rank)
    print("Singular values of A:\n", s)

    # Verify the lstsq solution
    print("\nVerification (A @ x_lstsq):\n", A @ x_lstsq)

    # Note: x_particular_pinv and x_lstsq should be identical for consistent underdetermined systems.
    print("\nComparison of solutions:")
    print(f"Solution from pinv: {x_particular_pinv}")
    print(f"Solution from lstsq: {x_lstsq}")
    print(f"Are they close? {np.allclose(x_particular_pinv, x_lstsq)}")

Matrix A:
 [[ 1  2 -1]
 [ 3 -1  2]]

Vector b:
 [5 1]

Rank of A: 2
Rank of augmented A: 2

Solutions exist for this system (potentially infinite).

--- Solving using np.linalg.pinv (Pseudoinverse) ---

Pseudoinverse of A:
 [[ 0.20481928  0.22891566]
 [ 0.3253012  -0.04819277]
 [-0.14457831  0.13253012]]

Particular solution (minimum-norm) from pinv:
 [ 1.25301205  1.57831325 -0.59036145]

Verification (A @ x_particular_pinv):
 [5. 1.]

Basis for the null space of A:
 [[-0.32929278]
 [ 0.5488213 ]
 [ 0.76834982]]

General Solution Form (using pinv's particular solution):
x = x_particular_pinv + c * null_space_basis_vector
x = [ 1.25301205  1.57831325 -0.59036145] + c * [-0.32929278  0.5488213   0.76834982]

Example solution with c = 10:
 [-2.03991575  7.06652625  7.09313675]

Verification (A @ another_solution):
 [5. 1.]

Example solution with c = -5:
 [ 2.89947595 -1.16579325 -4.43211055]

Verification (A @ yet_another_solution):
 [5. 1.]

--- Solving using np.linalg.lstsq ---

Solution from np.linalg.lstsq (minimum-norm):
 [ 1.25301205  1.57831325 -0.59036145]
Residuals (should be close to zero for consistent systems):
 []
Rank of A from lstsq:
 2
Singular values of A:
 [3.75807206 2.42423068]

Verification (A @ x_lstsq):
 [5. 1.]

Comparison of solutions:
Solution from pinv: [ 1.25301205  1.57831325 -0.59036145]
Solution from lstsq: [ 1.25301205  1.57831325 -0.59036145]
Are they close? True