numpy.linalg.solve: Solve Ax = b in Python

numpy.linalg.solve(a, b) solves the linear system ax = b for x when a is a square, full-rank matrix. It is the function to use instead of computing an inverse and multiplying. This page is for researchers and analysts who need to solve a system of equations in Python, read the output shapes correctly (the rules changed in NumPy 2.0) and diagnose the errors it raises.

What numpy.linalg.solve is used for in research

Any time a model reduces to a set of simultaneous linear equations, solve returns the exact solution of the well-determined case. Typical research uses are the normal equations of a regression, steady-state or balance equations, finite-difference systems and calibration equations. According to the NumPy reference, the solution is computed with the LAPACK routine _gesv, and a must be square and of full rank, meaning all rows (equivalently, all columns) are linearly independent.

If the system is over- or under-determined or rank-deficient, the NumPy documentation points to numpy.linalg.lstsq for the least-squares best solution instead. For the matrix-product side of the same library see matrix multiplication in NumPy, and for vector and matrix magnitudes see numpy.linalg.norm.

Install or access

numpy.linalg.solve ships with NumPy; there is nothing separate to install. Install NumPy from the official project (numpy.org) with pip install numpy or your conda distribution, then import it:

import numpy as np
print(np.__version__)

Check the version, because the handling of b changed in NumPy 2.0 (see the shape table below).

Step-by-step: solve a system of equations

Signature: numpy.linalg.solve(a, b). The parameter a is the coefficient matrix with shape (..., M, M); b is the right-hand side with shape (M,) or (..., M, K).

1. Write the system as a matrix equation. For 2x + y – z = 8, -3x – y + 2z = -11 and -2x + y + 2z = -3, the coefficient matrix is the 3 x 3 array of the coefficients and b is the vector of right-hand sides.
2. Build the arrays (floats avoid integer surprises):

   import numpy as np
   a = np.array([[2, 1, -1],
                 [-3, -1, 2],
                 [-2, 1, 2]], dtype=float)
   b = np.array([8, -11, -3], dtype=float)

3. Call solve:

   x = np.linalg.solve(a, b)
   print(x)   # [ 2.  3. -1.]

4. Check the residual. The NumPy reference example verifies a solution with np.allclose:

   np.allclose(a @ x, b)   # True

The two-equation example in the NumPy reference (x0 + 2*x1 = 1 and 3*x0 + 5*x1 = 2, with a = [[1, 2], [3, 5]] and b = [1, 2]) returns array([-1., 1.]). The 3 x 3 system above was run on NumPy 2.5.1 for this page and returned 2, 3 and -1.

Several right-hand sides at once

Pass b as an (M, K) matrix to solve K systems that share one a. Here np.linalg.solve(a, np.column_stack([b, 2*b])) returns shape (3, 2), one solution per column.

Output and how to read it

The returned x has shape (..., M) if b has shape (M,), and (..., M, K) if b has shape (..., M, K), where the leading dimensions are broadcast between a and b. The NumPy reference states that since version 2.0, b is treated as a single column vector only if it is exactly 1-dimensional; in every other case it is treated as a stack of (M, K) matrices. Before 2.0 it was treated as a stack of (M,) vectors when b.ndim == a.ndim - 1.

Shape of a	Shape of b	Shape of x	Reading
(M, M)	(M,)	(M,)	One system, 1-D solution
(M, M)	(M, K)	(M, K)	K systems, one solution per column
(N, M, M)	(N, M, K)	(N, M, K)	A stack of N matrices, each solved against its own (M, K) block

A practical check: a (M, 1) column returns (M, 1), not (M,) (confirmed by running np.linalg.solve(a, b.reshape(3, 1)).shape on NumPy 2.5.1, which gives (3, 1)). If downstream code expects a flat vector, use x.ravel() or pass a 1-D b. When porting older code that passed a stacked b of shape (N, M) against a of shape (N, M, M), re-check shapes under NumPy 2.0 or later, because that case is no longer read as a stack of vectors.

solve does not report how close to singular the matrix is. A solution can still be inaccurate for an ill-conditioned a; compare with np.linalg.cond(a) (the example above has a condition number of about 45) and always run the residual check.

Troubleshooting

Error	Cause	Fix
LinAlgError: Singular matrix	The NumPy reference says LinAlgError is raised if a is singular. Example: [[1, 2], [2, 4]], whose second row is twice the first.	Use np.linalg.lstsq(a, b, rcond=None) for a least-squares answer, or remove the redundant equation.
LinAlgError: Last 2 dimensions of the array must be square	a is not square (the documented “not square” case); reproduced with a 2 x 3 array.	Over- or under-determined system: use lstsq.
ValueError about mismatched core dimension	The length of b does not match M (reproduced with a length-2 b against a 3 x 3 a).	Make b have M rows.

For the singular example above, lstsq returned [0.2, 0.4], the minimum-norm least-squares solution. If you are fitting a model with nonlinear terms rather than solving a linear system, see SciPy least_squares. The Python-side routines call into optimised linear-algebra libraries; background is in BLAS.

Free alternatives and related functions

Function	Use when	Source note
numpy.linalg.solve	Square, full-rank a; exact solution	NumPy reference
numpy.linalg.lstsq	Non-square, or rank-deficient a; least-squares solution	Named in the NumPy solve notes
scipy.linalg.solve	You want structure hints or extra options: signature solve(a, b, lower=False, overwrite_a=False, overwrite_b=False, check_finite=True, assume_a=None, transposed=False)	SciPy reference; NumPy lists it under See also

The SciPy documentation describes assume_a as a hint about the matrix structure (diagonal, tridiagonal, banded, upper or lower triangular, symmetric, hermitian, positive definite or general); when it is None the structure is detected automatically. SciPy raises LinAlgError for a singular matrix, ValueError for a size mismatch or non-square input, and LinAlgWarning for an ill-conditioned one. Setting check_finite=False skips the check for NaN and infinity and can crash or return meaningless output on such input, so leave it on unless the data are already validated. If you set assume_a explicitly, make sure it is true of your matrix, as the function relies on it. For GPU arrays with a NumPy-like interface see CuPy.

For MATLAB users, the equivalent concerns (backslash versus an explicit inverse) are covered in MATLAB inverse matrix.

FAQ

What does numpy.linalg.solve do?

It computes the exact solution x of ax = b for a square, full-rank coefficient matrix a, using LAPACK _gesv.

What if my matrix is not square?

solve raises LinAlgError. Use numpy.linalg.lstsq for the least-squares solution.

Did anything change in NumPy 2.0?

Yes. b is a column vector only if exactly 1-D; otherwise it is a stack of (M, K) matrices.

Is there a SciPy version?

Yes, scipy.linalg.solve, with assume_a, check_finite and overwrite options.

Sources

• NumPy documentation, numpy.linalg.solve, read 2026-11-08.
• SciPy documentation, scipy.linalg.solve, read 2026-11-08.
• Outputs labelled as run were produced by the page author with NumPy 2.5.1.