Solving Linear Systems with NumPy

NumPy provides powerful tools for solving systems of linear equations, a foundational task in scientific computing, engineering, and machine learning. With functions like np.linalg.solve(), you can directly solve equations of the form Ax = b without manually computing matrix inverses.

Here's a quick overview of key concepts involved:

  • System of Equations: Represents multiple linear equations, which can be written in matrix form as Ax = b.
  • Matrix Form: A is the coefficient matrix, x is the unknown vector, and b is the output or result vector.
  • Efficient Solving: np.linalg.solve() is optimized and more stable than using the matrix inverse directly.

Understanding how to model and solve these systems using NumPy is essential for working with real-world data, fitting models, and implementing algorithms efficiently. Let’s dive in and see how it works in practice.

What is a System of Linear Equations?

A system of linear equations consists of two or more linear equations involving the same set of variables. These equations are solved simultaneously to find values that satisfy all of them at once.

For example, a simple system with two variables might look like:

2x + 3y = 8

3x + y = 5

Solving this system means finding values of x and y that make both equations true. Such systems are foundational in fields like data science, physics, economics, and engineering.

Representing the System in Matrix Form

A system of linear equations can be rewritten in a compact form using matrices. This is known as the matrix form of the system, expressed as Ax = b, where:

  • A is the coefficient matrix, containing the constants from the left-hand side of the equations.
  • x is the variable vector, representing the unknowns (e.g. x and y).
  • b is the output (result) vector, containing the constants from the right-hand side of the equations.

For the example:

2x + 3y = 8

3x + y = 5

This can be written in matrix form as Ax = b, where:

python
import numpy as np

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

# Output vector b
b = np.array([8, 5])

This form makes it easy to perform efficient computations using NumPy’s linear algebra functions, especially for larger systems with many variables.

Solving the System with NumPy

Once the system is written in matrix form Ax = b, you can solve for the unknown vector x using the np.linalg.solve() function in NumPy.

This method is preferred over computing the inverse manually because it’s faster and numerically more stable.

python
import numpy as np

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

# Output vector b
b = np.array([8, 5])

# Solve for x
x = np.linalg.solve(A, b)

print("Solution:", x)

How It Works:

  • np.linalg.solve(A, b) solves the matrix equation Ax = b directly for x using optimized numerical algorithms.
  • This avoids computing the inverse of A, which can be slower and less stable numerically.
  • It returns a solution vector that satisfies all the equations in the system.
  • If A is not square or has a zero determinant (i.e., it's singular), NumPy will raise a LinAlgError.

Output

Solution: [1. 2.]

This means the solution to the system is x = 1 and y = 2. NumPy’s solver works by applying efficient algorithms under the hood to find this solution directly.

💡 Tip: Only use np.linalg.solve() when A is a square matrix (same number of equations as unknowns) and has a non-zero determinant.

What If the Matrix Is Singular?

A matrix is singular if it is not invertible — typically because its rows or columns are linearly dependent. In such cases, the system of equations may have no solution or infinitely many solutions.

If you try to solve a system with a singular matrix using np.linalg.solve(), NumPy will raise a LinAlgError.

python
import numpy as np

# Singular matrix (second row is a multiple of the first)
A = np.array([[1, 2],
              [2, 4]])

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

try:
    x = np.linalg.solve(A, b)
    print("Solution:", x)
except np.linalg.LinAlgError as e:
    print("Error:", e)

How It Works:

  • NumPy detects that the matrix A is singular (its determinant is 0).
  • np.linalg.solve() cannot proceed and raises a LinAlgError.
  • This typically means the system has no unique solution — either inconsistent or dependent equations.
  • Always check for singularity if there's a possibility of linearly dependent rows or columns.

Output

Error: Singular matrix

💡 Tip: Use np.linalg.det(A) to check if the determinant is zero before attempting to solve.

Solving Overdetermined or Underdetermined Systems

Sometimes, systems have either more equations than unknowns (overdetermined) or more unknowns than equations (underdetermined). These systems may not have an exact solution.

NumPy’s np.linalg.lstsq() function computes the least-squares solution, which finds the best approximate solution that minimizes the sum of squared residuals.

python
import numpy as np

# Overdetermined system (3 equations, 2 unknowns)
A = np.array([[1, 1],
              [2, 1],
              [1, 3]])

b = np.array([2, 3, 4])

# Compute least-squares solution
x, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)

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

How It Works:

  • np.linalg.lstsq(A, b) finds the vector x that minimizes the squared differences between Ax and b.
  • This method works for both overdetermined and underdetermined systems.
  • residuals gives the sum of squared residuals (difference between actual and predicted).
  • rank shows the rank of matrix A, and s are singular values used internally.
  • Use this when no exact solution exists or system is not square.

Output

Least-squares solution: [0.71428571 1.14285714]
Residuals: [0.28571429]

💡 Tip: rcond=None uses NumPy’s default cutoff for small singular values; adjust if needed for numerical stability.

Frequently Asked Questions

How do I solve a system of linear equations using NumPy?

Use np.linalg.solve(A, b), where A is the coefficient matrix and b is the constants vector. This returns the solution x such that A @ x = b.

What happens if the matrix is singular or not square?

np.linalg.solve() requires a square and non-singular matrix. If the matrix is singular or rectangular, it raises a LinAlgError. You can use np.linalg.lstsq() to compute a least-squares solution in such cases.

Can I solve linear systems with multiple right-hand sides in NumPy?

Yes! If b is a 2D array with multiple columns, np.linalg.solve() solves for each right-hand side vector simultaneously, returning a solution matrix.

How do I interpret the output of np.linalg.lstsq()?

np.linalg.lstsq() returns a tuple where the first element is the solution vector minimizing the squared error, useful for overdetermined or singular systems.

Is there a way to check if a matrix is singular before solving?

You can compute the determinant with np.linalg.det(A). If it is zero or near zero, the matrix is singular or ill-conditioned, so np.linalg.solve() may fail.

What's Next?

Coming up next, we’ll explore the concept of the Matrix Rank in NumPy — an important property that tells us about the linear independence of rows or columns in a matrix. Understanding rank is essential for analyzing the solvability of linear systems and for gaining deeper insights into the structure of data in applied mathematics and machine learning.