How to Solve Linear Equations in NumPy

Key Insights

• Use numpy.linalg.solve() as your default method for solving linear systems—it’s faster and more numerically stable than computing matrix inverses
• Never invert matrices just to solve equations; it’s computationally expensive and introduces unnecessary numerical errors
• For overdetermined systems (more equations than unknowns), use numpy.linalg.lstsq() to find the best approximate solution

Introduction

Linear equations form the backbone of scientific computing. Whether you’re analyzing electrical circuits, fitting curves to data, balancing chemical equations, or training machine learning models, you’ll eventually need to solve systems of linear equations. NumPy provides battle-tested tools that handle these computations efficiently, leveraging optimized LAPACK and BLAS libraries under the hood.

This article walks you through the practical methods for solving linear equations in NumPy, from the straightforward solve() function to handling edge cases with least squares. You’ll learn when to use each approach and why some common practices (like computing matrix inverses) should be avoided.

Understanding Linear Equations in Matrix Form

Any system of linear equations can be expressed in matrix form as Ax = b, where:

• A is the coefficient matrix
• x is the vector of unknowns
• b is the vector of constants

Consider this simple system:

2x + 3y = 8
4x - y = 2

In matrix form, this becomes:

import numpy as np

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

# Constants vector b
b = np.array([8, 2])

# We're solving for x where Ax = b

The coefficient matrix A contains the multipliers for each variable, arranged so that each row represents one equation and each column represents one variable. The vector b holds the right-hand side values.

This matrix representation isn’t just mathematical elegance—it’s how NumPy efficiently processes these problems using vectorized operations and optimized linear algebra routines.

Using numpy.linalg.solve() for Direct Solutions

The numpy.linalg.solve() function is your primary tool for solving linear systems. It uses LU decomposition internally to find the exact solution when one exists.

import numpy as np

# System of 3 equations with 3 unknowns:
# x + 2y + 3z = 14
# 2x + 5y + 3z = 18
# 3x + 6y + 9z = 42

A = np.array([
    [1, 2, 3],
    [2, 5, 3],
    [3, 6, 9]
], dtype=float)

b = np.array([14, 18, 42], dtype=float)

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

print(f"Solution: x={x[0]:.2f}, y={x[1]:.2f}, z={x[2]:.2f}")
# Output: Solution: x=1.00, y=2.00, z=3.00

# Verify the solution
print(f"Verification (should equal b): {A @ x}")
# Output: Verification (should equal b): [14. 18. 42.]

The function requires that matrix A be square and non-singular (invertible). If these conditions aren’t met, NumPy raises a LinAlgError.

Why use solve() instead of alternatives?

1. Numerical stability: It uses partial pivoting to minimize rounding errors
2. Performance: O(n³) complexity, optimized at the BLAS level
3. Memory efficiency: Doesn’t create intermediate matrices

Always verify your solution by computing A @ x and comparing it to b. Small floating-point differences are normal; use np.allclose() for robust comparisons.

Alternative Methods: Inverse and LU Decomposition

You’ll often see tutorials suggesting the inverse method: if Ax = b, then x = A⁻¹b. While mathematically correct, this approach has significant drawbacks.

import numpy as np

A = np.array([
    [1, 2, 3],
    [2, 5, 3],
    [3, 6, 9]
], dtype=float)

b = np.array([14, 18, 42], dtype=float)

# Method 1: Using solve() - RECOMMENDED
x_solve = np.linalg.solve(A, b)

# Method 2: Using matrix inverse - NOT RECOMMENDED
A_inv = np.linalg.inv(A)
x_inv = A_inv @ b

print(f"solve() result: {x_solve}")
print(f"inverse result: {x_inv}")

Both produce the same answer, but the inverse method:

• Requires twice the computation (computing inverse, then multiplying)
• Introduces additional numerical errors
• Uses more memory for storing the inverse matrix

When is the inverse useful? Only when you need to solve the same system with many different b vectors. Even then, consider LU decomposition instead.

For more control over the decomposition process, use SciPy’s LU factorization:

import numpy as np
from scipy import linalg

A = np.array([
    [1, 2, 3],
    [2, 5, 3],
    [3, 6, 9]
], dtype=float)

b = np.array([14, 18, 42], dtype=float)

# Compute LU decomposition once
lu, piv = linalg.lu_factor(A)

# Solve for multiple right-hand sides efficiently
b1 = np.array([14, 18, 42], dtype=float)
b2 = np.array([7, 9, 21], dtype=float)

x1 = linalg.lu_solve((lu, piv), b1)
x2 = linalg.lu_solve((lu, piv), b2)

print(f"Solution for b1: {x1}")
print(f"Solution for b2: {x2}")

LU decomposition splits A into lower and upper triangular matrices. Once computed, solving for any b vector becomes a fast back-substitution operation.

Handling Special Cases

Real-world problems don’t always give you nice square systems with unique solutions. Here’s how to handle the edge cases.

Overdetermined Systems (More Equations Than Unknowns)

When you have more equations than variables, there’s typically no exact solution. Use least squares to find the best approximation that minimizes the sum of squared residuals.

import numpy as np

# Fitting a line y = mx + c to 5 data points
# This gives us 5 equations but only 2 unknowns (m and c)

x_data = np.array([1, 2, 3, 4, 5])
y_data = np.array([2.1, 3.9, 6.2, 7.8, 10.1])

# Build the coefficient matrix for y = mx + c
# Each row: [x_i, 1] corresponding to m*x_i + c*1 = y_i
A = np.column_stack([x_data, np.ones(len(x_data))])

# Solve using least squares
solution, residuals, rank, singular_values = np.linalg.lstsq(A, y_data, rcond=None)

m, c = solution
print(f"Best fit line: y = {m:.3f}x + {c:.3f}")
# Output: Best fit line: y = 2.000x + 0.060

# The residuals tell us how well the fit matches
print(f"Sum of squared residuals: {residuals[0]:.4f}")

The lstsq() function returns four values: the solution, residuals, matrix rank, and singular values. The rcond=None parameter uses machine precision for determining rank.

Singular and Near-Singular Matrices

When A is singular (determinant is zero), no unique solution exists. NumPy will raise a LinAlgError. Check the condition number to detect near-singular matrices:

import numpy as np

A = np.array([
    [1, 2],
    [2, 4]  # This row is just 2x the first row - singular!
], dtype=float)

# Check condition number before solving
cond = np.linalg.cond(A)
print(f"Condition number: {cond}")  # Will be very large (inf for singular)

# For near-singular matrices, use pseudoinverse
A_pinv = np.linalg.pinv(A)

A condition number above 10¹⁰ suggests numerical instability. Consider reformulating your problem or using higher precision arithmetic.

Practical Application: Electrical Circuit Analysis

Let’s apply these techniques to a real problem: analyzing a resistor network using Kirchhoff’s laws.

Consider a circuit with three loops and unknown currents I₁, I₂, and I₃. Applying Kirchhoff’s voltage law to each loop gives us:

import numpy as np

# Circuit parameters
R1, R2, R3, R4, R5 = 10, 20, 30, 15, 25  # Resistances in ohms
V1, V2 = 12, 9  # Voltage sources in volts

# Kirchhoff's voltage law equations:
# Loop 1: R1*I1 + R2*(I1-I2) = V1
# Loop 2: R2*(I2-I1) + R3*I2 + R4*(I2-I3) = 0
# Loop 3: R4*(I3-I2) + R5*I3 = -V2

# Rearranged into matrix form:
# (R1+R2)*I1 - R2*I2 + 0*I3 = V1
# -R2*I1 + (R2+R3+R4)*I2 - R4*I3 = 0
# 0*I1 - R4*I2 + (R4+R5)*I3 = -V2

A = np.array([
    [R1 + R2, -R2, 0],
    [-R2, R2 + R3 + R4, -R4],
    [0, -R4, R4 + R5]
], dtype=float)

b = np.array([V1, 0, -V2], dtype=float)

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

print("Circuit Analysis Results")
print("-" * 30)
print(f"I1 = {currents[0]*1000:.2f} mA")
print(f"I2 = {currents[1]*1000:.2f} mA")
print(f"I3 = {currents[2]*1000:.2f} mA")

# Verify power balance
power_supplied = V1 * currents[0] - V2 * currents[2]
power_dissipated = (R1 * currents[0]**2 + 
                   R2 * (currents[0] - currents[1])**2 +
                   R3 * currents[1]**2 +
                   R4 * (currents[1] - currents[2])**2 +
                   R5 * currents[2]**2)

print(f"\nPower supplied: {power_supplied*1000:.2f} mW")
print(f"Power dissipated: {power_dissipated*1000:.2f} mW")
print(f"Energy conserved: {np.isclose(power_supplied, power_dissipated)}")

This example demonstrates the complete workflow: translating a physical problem into matrix form, solving it, and verifying the results make physical sense.

Summary and Best Practices

Here’s your quick reference for choosing the right method:

Performance tips for large systems:

1. Use dtype=float64 explicitly—mixed types trigger conversions
2. Ensure arrays are contiguous in memory with np.ascontiguousarray()
3. For sparse matrices, switch to scipy.sparse.linalg solvers
4. Profile before optimizing; NumPy’s defaults are already highly optimized

Avoid computing matrix inverses for solving equations. Use solve() or decomposition methods instead. Your code will be faster, more stable, and more maintainable.