SciPy - Solving Linear Systems

In linear algebra solving linear systems refers to finding the solution(s) to a system of linear equations. A linear system can be represented as follows −

A . x = b

Where −

• A is a matrix of coefficients).
• x is the vector of variables (unknowns)
• b is the vector of constants.

In SciPy solving linear systems is done using several methods depending on the type and properties of the matrix A. SciPy provides highly optimized functions to solve linear systems directly or via matrix decompositions such as LU, Cholesky, QR and others.

SciPy Function to Solve Linear Systems

The primary function in SciPy to solve a linear system is scipy.linalg.solve(). This function is used to compute the solution x of the system Ax=b where A is the matrix of coefficients and b is the vector of constants.

Syntax

Below is the syntax of the function scipy.linalg.solve() which is used to Solve Linear systems −

scipy.linalg.solve(
   a, 
   b, 
   sym_pos=False, 
   lower=False, 
   overwrite_a=False, 
   overwrite_b=False, 
   check_finite=True
)

• a(array_like): The coefficient matrix A which should be a square matrix i.e., number of rows = number of columns and in most cases non-singular i.e., invertible.
• b(array_like): The right-hand side vector or matrix i.e., b which must have the same number of rows as the matrix A.
• sym_pos(bool, optional): If True then it indicates that the matrix A is symmetric and positive-definite. In this case the function uses more efficient algorithms such as Cholesky decomposition.
• lower(bool, optional): If True then it indicates that the matrix A is lower triangular. This will reduce the computational cost.
• overwrite_a(bool, optional): If True then allows overwriting the input matrix A in memory to save space which perform faster computations but changes the original matrix.
• overwrite_b(bool, optional): If True then allows overwriting the input vector b in memory.
• check_finite(bool, optional): If True then the function checks whether A and b contain only finite numbers. This helps to avoid numerical issues but may incur some overhead.

This function returns the solution vector or matrix i.e., the values of the variables that satisfy the linear system Ax=b.

Example of Solving a Linear System

Following is the example of scipy.linalg.solve() function which is used to solve a system of linear equations with the shape 2 x 2 −

import numpy as np
from scipy.linalg import solve

# Define matrix A (2x2)
A = np.array([[3, 2],
              [1, 2]])

# Define vector b (2x1)
b = np.array([5, 5])

# Solve for x in Ax = b
x = solve(A, b)

print("Solution x:", x)

Following is the output of the above code −

Solution x: [0.  2.5]

Solving Linear Systems for Non-Square Matrices

In this example we have an over-determined system i.e., more equations than unknowns and we use the least-squares method to solve it −

import numpy as np
from scipy.linalg import solve

# Define the overdetermined matrix A (3x2 matrix)
A = np.array([[1, 2],
              [2, 3],
              [3, 4]])

# Define the right-hand side vector b
b = np.array([5, 6, 7])

# Solve using least-squares method
x = solve(A.T @ A, A.T @ b)

print("Least-squares solution x:", x)

Following is the output of the above code −

Least-squares solution x: [-3.  4.]

Solving a Singular System

Following is the example in which we have a singular matrix i.e., the system does not have a unique solution x + y = 1 and x + y = 2 −

import numpy as np
from scipy.linalg import solve

# Define a singular matrix A (rows are linearly dependent)
A = np.array([[1, 1],
              [1, 1]])

# Define the right-hand side vector b
b = np.array([1, 2])

# Attempt to solve the system A * x = b
try:
    x = solve(A, b)
    print("Solution x:", x)
except Exception as e:
    print(f"Error: {e}")

Following is the output of the above code −

Error: Matrix is singular.