Numerical Solving for AI: Use Solve Without Explicit Inverses

An AI calibration step may need to recover two unknown parameters from two measured responses. The equations are small enough to write on a whiteboard, but production code still has to decide how to solve them and how to check the answer. A plausible-looking vector is not evidence that the data determined it reliably.

This lesson develops that decision with NumPy. You will solve a square linear system, compare numpy.linalg.solve with the explicit-inverse route, compute a residual, measure conditioning, and recognize the singular boundary where a unique parameter estimate does not exist.

Learning outcome: Use numpy.linalg.solve and explain why explicit inversion is usually avoided.

Lesson 22: Identity and inverse transformations showed what an inverse means mathematically. Lesson 24: Row-echelon form reveals structure showed how pivots and free variables expose what a system can determine. Here, we connect that exact algebra to finite-precision computation. The official Mathematics for Machine Learning companion places this work in Section 2.3 of the book PDF; the calibration scenario, numbers, code, visual, and assessments below are original.

A linear system is a parameter-recovery contract

Write the unknown parameter vector as

x = \begin{bmatrix}x_1\\x_2\end{bmatrix}\in\mathbb{R}^2.

The vector has two coordinates, so each row of a coefficient matrix $A$ must contain two coefficients. If two observations produce the right-hand side vector $b$ , a parameter-recovery system has the form

A x=b, \qquad A\in\mathbb{R}^{2\times 2}, \quad b\in\mathbb{R}^{2}, \quad x\in\mathbb{R}^{2}.

For a tiny calibration head, $x_1$ might be a scale applied to one signal and $x_2$ an offset-like correction applied to another. More generally, the same notation describes a small linear layer, a sensor correction, or a parameter recovery subproblem inside a larger AI pipeline. The rows of $A$ are the observed coefficient patterns; $b$ contains the responses those patterns should produce.

If $A$ is square and invertible, the algebra from the previous lesson gives

x=A^{-1}b.

That identity is useful for reasoning. It does not mean that a program should first construct $A^{-1}$ every time it needs $x$ .

Worked example: recover two calibration parameters

Consider the deliberately small system

A=\begin{bmatrix}3&1\\1&2\end{bmatrix}, \qquad b=\begin{bmatrix}8\\1\end{bmatrix}.

The two equations are

3x_1+x_2=8, \qquad x_1+2x_2=1.

Solve it by substitution before asking NumPy to do the same work. The first equation gives $x_2=8-3x_1$ . Substitute that expression into the second:

x_1+2(8-3x_1)=1 \Longrightarrow -5x_1=-15 \Longrightarrow x_1=3.

Then $x_2=8-3(3)=-1$ , so

x=\begin{bmatrix}3\\-1\end{bmatrix}.

Check both original equations, not only the rearranged one:

A x = \begin{bmatrix}3&1\\1&2\end{bmatrix} \begin{bmatrix}3\\-1\end{bmatrix} = \begin{bmatrix}8\\1\end{bmatrix}=b.

The same calculation in NumPy

import numpy as np

A = np.array([[3.0, 1.0],
              [1.0, 2.0]])
b = np.array([8.0, 1.0])

x_solve = np.linalg.solve(A, b)
x_inverse = np.linalg.inv(A) @ b

print(x_solve)
print(x_inverse)
print(np.linalg.norm(A @ x_solve - b))
print(np.linalg.norm(A @ x_inverse - b))

On the bundled NumPy 2.3.5 runtime, the verified output is:

[ 3. -1.]
[ 3. -1.]
2.220446049250313e-16
9.155133597044475e-16

The tiny nonzero values are floating-point residual norms. In exact arithmetic both residuals are zero, and both methods represent the same mathematical answer. The slight difference in the last digits is a reminder to compare floating-point results with a tolerance rather than with exact equality.

Why `solve` is usually the right intermediate

The direct mathematical expression $A^{-1}b$ hides two different computational tasks:

Route	What it computes	Why that matters for one right-hand side
`np.linalg.solve(A, b)`	Factorizes $A$ and applies the factors to $b$	Computes only the vector requested.
`np.linalg.inv(A) @ b`	Builds every column of $A^{-1}$ , then multiplies by $b$	Computes much more information than this one solve needs.

For a single right-hand side, solve avoids materializing an entire inverse. For several right-hand sides stored as the columns of $B$ , use np.linalg.solve(A, B) so the same coefficient matrix can be factored once while all requested solutions are computed. An inverse is an object to keep only when the application genuinely needs the inverse itself; it is not a default substitute for a solve.

The current NumPy documentation for solve says that the coefficient matrix must be square and full rank, and documents that the routine uses a LAPACK _gesv implementation. The documentation for inv also warns that an ill-conditioned matrix can produce inaccurate results even when inversion does not raise an error. That is why the implementation choice and the diagnostics belong in the same mental model.

The runtime difference is hardware- and library-dependent, so it is not a portable benchmark. As an illustration, the following authoring-host run used a fixed random $64\times64$ matrix, five timing repeats, and twenty calls per repeat:

from timeit import repeat

rng = np.random.default_rng(7)
A_timing = rng.normal(size=(64, 64)) + 64 * np.eye(64)
b_timing = rng.normal(size=64)

solve_time = min(
    repeat(lambda: np.linalg.solve(A_timing, b_timing), repeat=5, number=20)
) / 20
inverse_time = min(
    repeat(lambda: np.linalg.inv(A_timing) @ b_timing, repeat=5, number=20)
) / 20

print(solve_time * 1e6, inverse_time * 1e6, inverse_time / solve_time)

The verified run printed approximately 15.2 µs, 52.2 µs, and 3.43× for the solve time, inverse-route time, and ratio. Treat those numbers as a reproducible shape of the comparison, not as a performance promise: the algorithmic reason to prefer solve is that it avoids computing an unused matrix, and a different BLAS library or matrix size can change the ratio.

A residual is a fit check, not a trust certificate

For the worked example, calculate the residual from the direct route:

r=b-A\hat{x} = \begin{bmatrix}8\\1\end{bmatrix} - \begin{bmatrix}8\\1\end{bmatrix} = \begin{bmatrix}0\\0\end{bmatrix}.

The floating-point program reports a norm near $10^{-15}$ instead of exactly zero. That is excellent evidence that the returned vector fits this $b$ . However, imagine that the coefficients or responses came from rounded sensor measurements. A solver can find a vector that fits those rounded values very closely while still moving substantially when the measurements are changed within their uncertainty.

Conditioning measures sensitivity

The condition number compares how much a matrix can amplify relative input changes. For an invertible matrix and the 2-norm,

\kappa_2(A)=\lVert A\rVert_2\,\lVert A^{-1}\rVert_2 =\frac{\sigma_{\max}(A)}{\sigma_{\min}(A)},

where $\sigma_{\max}$ and $\sigma_{\min}$ are the largest and smallest singular values. A condition number near $1$ indicates balanced directions. A large value means that some small change in the inputs can cause a much larger relative change in the solution. If the smallest singular value is zero, the matrix is singular and the condition number is infinite.

For a perturbed right-hand side, a useful bound has the form

\frac{\lVert\delta x\rVert}{\lVert x\rVert} \lesssim \kappa(A)\frac{\lVert\delta b\rVert}{\lVert b\rVert},

under the usual assumptions for a nonsingular system. The inequality is a warning about possible amplification, not a promise that every perturbation will reach the worst case.

Now use nearly dependent rows:

B=\begin{bmatrix}1&1\\1&1.0001\end{bmatrix}, \qquad x_{\text{true}}=\begin{bmatrix}2\\-1\end{bmatrix}, \qquad b=B x_{\text{true}}=\begin{bmatrix}1\\0.9999\end{bmatrix}.

The second row is almost a copy of the first. With the exact displayed $b$ , the solve recovers $[2,-1]^\mathsf{T}$ . Change only the second response by $10^{-6}$ , then solve again. NumPy reports

condition number: 40002.000074915224
original solution: [ 2. -1.]
perturbed solution: [ 1.99 -0.99]
relative input change: 7.071421374297804e-07
relative solution change: 0.006324555320519327
perturbed residual norm: 0.0

The input changed by about $7.1\times10^{-7}$ relative, while the parameter vector changed by about $6.3\times10^{-3}$ relative. That is roughly a four-order-of-magnitude condition number expressing itself as sensitivity. The perturbed solution still has zero residual for the perturbed data: it solved the new equations correctly, but the data did not identify a stable parameter vector.

Static stability map

A small residual can hide a sensitive solution

The residual checks the equation that was solved. The condition number warns how much the solution can move when the input data move.

The same solve pipeline can produce a tiny residual in all three cases. Conditioning determines whether that fit is also a trustworthy parameter estimate.

Exact quantities represented in the stability map
Case	Condition number	Relative input → solution change	Residual and interpretation
Well-conditioned	2.62	1.2 × 10⁻⁷ → 2.0 × 10⁻⁷	Residuals around 10⁻¹⁵; both routes agree for this well-conditioned example.
Ill-conditioned	4.0002 × 10⁴	7.1 × 10⁻⁷ → 6.3 × 10⁻³	Residual 0 for the perturbed data, but a 7.1 × 10⁻⁷ input change moves x by 6.3 × 10⁻³.
Singular	∞	Not applicable	κ₂(A) is infinite; solve and inverse raise a singular-matrix error instead of returning a unique answer.

Read the two diagnostics together: residual asks “does this x fit b?”; conditioning asks “how much would x change if b or A were slightly different?”

The visual makes the two questions separate. The residual asks whether the candidate fits the equations currently in memory. The condition number asks whether a small change in the equations could move that candidate. An AI calibration system needs both answers.

Boundary case: singular means “no unique solve”

Consider

C=\begin{bmatrix}1&2\\2&4\end{bmatrix}.

The second row is twice the first, so $C$ has rank $1$ , not $2$ . With $b=[3,6]^\mathsf{T}$ , the equations are consistent but have infinitely many solutions: every $(x_1,x_2)$ satisfying $x_1+2x_2=3$ works. With $b=[3,7]^\mathsf{T}$ , they are inconsistent because the second equation asks for twice the first left-hand side but gives a different right-hand side.

Neither situation has a unique vector for np.linalg.solve to return. In the bundled runtime, both np.linalg.solve(C, [3, 6]) and np.linalg.inv(C) raise numpy.linalg.LinAlgError: Singular matrix. The NumPy solve documentation specifies the full-rank requirement and points to numpy.linalg.lstsq when a least-squares solution is appropriate for a non-square or rank-deficient problem. That is a different question from solving an invertible square system, so do not silently replace one method with the other.

A disciplined NumPy workflow for AI systems

For a small parameter-recovery subproblem, a defensible workflow is:

Check the contract. Confirm that the dimensions of $A$ , $b$ , and the intended parameter vector agree.
Solve directly. Use np.linalg.solve(A, b) for a square, full-rank system; use a method designed for least squares when the problem is not well-determined.
Measure the residual. Compute np.linalg.norm(A @ x - b) and compare it with the scale of b, not with an arbitrary absolute threshold alone.
Inspect conditioning. Compute np.linalg.cond(A) and compare the possible sensitivity with measurement precision and the number of digits you plan to use downstream.
Treat failure as information. A singular error, an enormous condition number, or a residual that is large relative to b should trigger a data, feature, scaling, or model-design investigation.

In a calibration or parameter-recovery task, this workflow turns a vector of weights into an auditable result. A tiny residual can support the claim that the implementation solved the supplied equations. A moderate condition number and a perturbation check provide evidence that those equations constrain the parameters strongly enough for the intended use. Neither diagnostic proves that the broader AI system is accurate; they answer narrower implementation questions that can be tested.

The next lesson steps away from numerical algorithms and asks which algebraic operations are valid over the number system used by common AI tensors. For now, keep the practical rule: solve the system you have, check the residual you got, and ask conditioning whether a nearby system would tell a different parameter story.