Pseudoinverse: choose the smallest least-squares solution

The Moore–Penrose pseudoinverse chooses the smallest input among all least-squares solutions. It works when a system has one exact solution, many exact solutions, or an unreachable target.

Its two criteria concern different spaces: the residual measures error in the output, while the solution norm measures input length. This lesson keeps both visible through small, synthetic examples.

Make two choices in order

Let A be a real m × n matrix and b an m-component target. Its pseudoinverse, written A⁺, has shape n × m. The vector x* = A⁺b follows this rule:

1. Minimize the least-squares residual norm ‖b − Ax‖₂.
2. Among all inputs that attain that minimum, choose the one with smallest ‖x‖₂.

This produces a unique x* for each b. If an exact solution exists, the minimum residual is zero. If A is square and invertible, A⁺ equals the ordinary inverse.

Here “smallest” means Euclidean norm in the chosen coordinates. It does not mean minimizing a sum of residual error and input size with an unspecified tradeoff. MIT's linear-algebra notes describe the minimum-norm least-squares choice.

Coordinate scales matter. A robot's joint-velocity vector may mix angular and linear rates, or joints with different limits. An unweighted Euclidean minimum does not automatically minimize energy, satisfy limits, or produce safe motion; those goals need their own physical model and constraints.

The Jacobian matrices lesson builds a concrete joint-velocity map for a two-link arm. Its changing rank shows why some tip velocities become unreachable at a particular configuration.

Choose one input from an exact solution family

Use the 2 × 3 matrix from the rank and null-space lesson:

Every exact solution has the form x = (2 − t, 1 − t, t). The parameter changes the input along null direction (−1, −1, 1), so Ax stays (2, 1).

Compare their squared lengths:

The squared term reaches its minimum at t = 1. Thus A⁺b = (1, 0, 1), with residual zero and norm √2 ≈ 1.4142. At t = 0, input (2, 1, 0) fits equally well but has norm √5 ≈ 2.2361.

For this matrix, the pseudoinverse is:

Multiplying it by (2, 1) gives ((4 − 1)/3, (−2 + 2)/3, (2 + 1)/3) = (1, 0, 1). This checks the matrix formula against the minimum found from the solution family.

Compare fit error with input length

Start at t = 1, then move Free parameter t to zero. The fitted output stays fixed, while the input bars and solution norm change. Choose minimum norm returns to A⁺b.

Interactive experiment

Compare equally good fits by input length

Every displayed input minimizes the least-squares residual for the selected example. Move along its null directions to see which input has minimum Euclidean norm.

Pseudoinverse exampleWide: exact solution familyTall: inconsistent targetRank one: dependent columnsZero map: every input ties on fit

A: 2 × 3

1	0	1
0	1	1

A⁺: 3 × 2

0.667	-0.333
-0.333	0.667
0.333	0.333

Displayed entries are rounded; calculations keep additional precision.

x = (2 − t, 1 − t, t)

Free parameter t

1.0

Input coordinates

Current input x

(1.000, 0.000, 1.000)

Pseudoinverse input A⁺b

(1.000, 0.000, 1.000)

Residual norm ‖b − Ax‖₂

0.0000

Solution norm ‖x‖₂

1.4142

Minimum solution norm

1.4142

Minimum norm choice

Yes

This is the minimum-norm least-squares input. Residual norm: 0.0000. Input norm: 1.4142.

Target and fitted output; residual is b − Ax

Output coordinate	b	Ax	b − Ax
1	2.000	2.000	0.000
2	1.000	1.000	0.000

Choose minimum normReset experiment

These are fixed analytic examples. The browser evaluates their solution families; it does not run a numerical SVD or choose a singular-value cutoff.

The bars show individual input coordinates, including all three coordinates of the wide example. They do not represent a two-dimensional projection of a three-dimensional vector. The table separately reports the target, fitted output, and signed residual b − Ax.

Every input displayed by the sliders is already a least-squares minimizer. The sliders explore null-space freedom within that set. For the tall example, the least-squares solution is unique, so there is no free slider.

The browser uses analytic pseudoinverses for four fixed matrices. Matrix entries displayed as decimals may be rounded. This experiment does not compute a general numerical SVD or select a singular-value threshold.

Fit an unreachable target

Choose the tall 3 × 2 matrix T = [1, 0; 0, 1; 1, 1] and target b = (2, 1, 4). Its output always has the form (x₁, x₂, x₁ + x₂). The target's third coordinate differs from 2 + 1, so no exact solution exists.

The pseudoinverse gives x* = (7/3, 4/3). Its fitted output and residual are:

The residual norm is 1/√3 ≈ 0.5774. The input norm is √65/3 ≈ 2.6874. These numbers answer different questions and need not match.

Check r against T's two columns. Both dot products are zero: −1/3 + 1/3 = 0. The fitted output is the orthogonal projection of b onto the column space, as explained in projections and least squares.

T has full column rank, so its least-squares input is unique. The minimum-norm rule needs no further choice in this case.

Handle dependent columns and zero maps

Choose R = [1, 2; 2, 4], whose columns are dependent. For b = (2, 4), all exact solutions are x = (2 − 2t, t). Their squared norms satisfy:

The minimum occurs at t = 0.8, giving x* = (0.4, 0.8). Its dot product with null direction (−2, 1) is −0.8 + 0.8 = 0.

This illustrates the general geometry. The minimum-norm least-squares input lies in the row space of A, which is orthogonal to Null(A). Every other least-squares input is x* + z for a null vector z, so the Pythagorean relation gives ‖x* + z‖₂² = ‖x*‖₂² + ‖z‖₂².

The zero map supplies a useful boundary case. With A = 0 and b = (2, 1), every input produces zero output and residual (2, 1). All inputs tie on fit, but A⁺ = 0 selects x* = 0, the unique input with zero norm.

Read the two projector products

The products involving A and A⁺ have specific geometric meanings:

• AA⁺, an m × m matrix, orthogonally projects onto the column space of A.
• A⁺A, an n × n matrix, orthogonally projects onto the row space of A.

For the wide example, AA⁺ = I₂ because its column space is all of ℝ². But A⁺A is a 3 × 3 projector that sends null direction (−1, −1, 1) to zero. Applying A and then A⁺ cannot recover an input component that A discarded.

These products satisfy P² = P and Pᵀ = P. Applying an orthogonal projection twice has the same effect as applying it once. MIT's pseudoinverse notes connect both products to the fundamental subspaces.

For real matrices, four conditions uniquely identify the Moore–Penrose pseudoinverse:

1. AA⁺A = A.
2. A⁺AA⁺ = A⁺.
3. (AA⁺)ᵀ = AA⁺.
4. (A⁺A)ᵀ = A⁺A.

SciPy documents these defining conditions. They remain valid for rectangular and rank-deficient matrices; a universal two-sided identity is not required.

Reverse the nonzero singular-value directions

The singular value decomposition writes A = UΣVᵀ. The orthogonal matrices describe input and output directions. The nonnegative singular values in Σ describe the gains along paired directions.

The pseudoinverse reverses each direction with a nonzero gain:

Σ⁺ has the transposed shape of Σ. Replace each nonzero singular value σ with 1/σ, while each exact zero keeps a zero entry. This sends output components outside the column space to zero and chooses no null-space input component. NumPy describes this construction.

For full column rank, one algebraic formula is A⁺ = (AᵀA)⁻¹Aᵀ. For full row rank, it is A⁺ = Aᵀ(AAᵀ)⁻¹. Each formula requires its stated rank condition; neither provides a general recipe for a rank-deficient matrix.

For numerical work, use a suitable least-squares or pseudoinverse routine that applies a matrix factorization. Forming these inverses by hand is useful for the small exact examples here. MIT's inverse lecture develops the full-rank formulas.

Treat a cutoff as a model choice

A small nonzero singular value has a large reciprocal. Suppose D = diag(1, 0.001) and b = (1, 1). Its exact pseudoinverse gives x* = (1, 1000) and zero residual.

If a cutoff discards 0.001, the effective matrix becomes D_cut = diag(1, 0). Its pseudoinverse gives x_cut = (1, 0). Against the original D, that input has residual (0, 1) and residual norm 1.

This truncated result solves the minimum-norm least-squares problem for the modified matrix. It is not the true pseudoinverse solution of the original full-rank D. It trades away a weak direction that would require a large input.

NumPy specifies its cutoff relative to the largest singular value, while SciPy documents absolute and relative cutoff terms. Choose thresholds with coordinate scaling, numerical precision, and measurement uncertainty in mind.

Minimum norm itself does not guarantee a modest input: (1, 1000) is already the smallest least-squares input for the unmodified D. Limiting amplification requires an explicit change to the problem, such as truncation, damping, or constraints.

The numerical inverse kinematics lesson applies damping to a robot arm’s local position problem. Take one update, measure the remaining error, and repeat with the Jacobian at the new pose.

Verify the examples with exact fractions

This standard-library example uses analytic pseudoinverses and Fraction arithmetic. It verifies all four Moore–Penrose identities for the wide matrix, then calculates the wide and tall solutions. It does not implement a numerical SVD.

from fractions import Fraction as F
from math import sqrt

def transpose(matrix):
    return tuple(zip(*matrix))

def matmul(left, right):
    if len(left[0]) != len(right):
        raise ValueError("Inner dimensions must match")
    return tuple(tuple(sum(a * b for a, b in zip(row, column))
                       for column in transpose(right)) for row in left)

def matvec(matrix, vector):
    return tuple(sum(a * b for a, b in zip(row, vector)) for row in matrix)

def norm_squared(vector):
    return sum(value * value for value in vector)

def display(vector):
    return "(" + ", ".join(str(value) for value in vector) + ")"

A = ((1, 0, 1), (0, 1, 1))
P = ((F(2, 3), F(-1, 3)),
     (F(-1, 3), F(2, 3)),
     (F(1, 3), F(1, 3)))
AP, PA = matmul(A, P), matmul(P, A)
assert matmul(AP, A) == A
assert matmul(PA, P) == P
assert transpose(AP) == AP
assert transpose(PA) == PA
print("Moore-Penrose checks: passed")

wide_x = matvec(P, (2, 1))
print("Wide x:", display(wide_x))
for t in [0, 1, 2]:
    x = (2 - t, 1 - t, t)
    assert matvec(A, x) == (2, 1)
    print(f"t={t}: squared norm={norm_squared(x)}")

T = transpose(A)
T_plus = transpose(P)
b = (2, 1, 4)
tall_x = matvec(T_plus, b)
fitted = matvec(T, tall_x)
residual = tuple(target - value for target, value in zip(b, fitted))
print("Tall x:", display(tall_x))
print("Tall residual b - T x:", display(residual))
print(f"Tall residual norm: {sqrt(float(norm_squared(residual))):.4f}")

Expected output:

Moore-Penrose checks: passed
Wide x: (1, 0, 1)
t=0: squared norm=5
t=1: squared norm=2
t=2: squared norm=5
Tall x: (7/3, 4/3)
Tall residual b - T x: (-1/3, -1/3, 1/3)
Tall residual norm: 0.5774

The identities and residual components use exact rational arithmetic. Only the final square root converts to a rounded decimal.

Try it yourself

Exercise 1. For the wide example, compare t = 0 and t = 2. Compute both inputs, fitted outputs, and squared input norms. Does their equal residual identify either input as the pseudoinverse choice?

Show solution: fit and input length are separate

The inputs are (2, 1, 0) and (0, −1, 2). Both produce (2, 1), so both have residual zero. Both have squared norm 5.

The pseudoinverse chooses t = 1, with input (1, 0, 1) and squared norm 2. Matching the best residual alone leaves the null-space freedom unresolved.

Exercise 2. Let D = diag(1, 0.01) and b = (2, 1). Find the true pseudoinverse solution. Then discard the second singular value and find the truncated result and its residual against the original D.

Show solution: identify which matrix is being inverted

The true pseudoinverse is diag(1, 100), giving x* = (2, 100) with residual zero. Discarding 0.01 gives D_cut⁺ = diag(1, 0) and x_cut = (2, 0).

For the original D, Dx_cut = (2, 0), so b − Dx_cut = (0, 1). The truncated input has a smaller norm, but it solves the modified problem and no longer minimizes the original residual.

Sources and further study

• MIT: Key Ideas in Linear Algebra, page 15, for minimum-norm least squares and the two projectors.
• MIT: Left and Right Inverses; Pseudoinverse, for full-rank formulas and the fundamental subspaces.
• SciPy: pinv, for the Moore–Penrose conditions and numerical cutoff controls.
• NumPy: pinv, for SVD construction and singular-value thresholds.