Moore-Penrose pseudoinverse

The formula

Every matrix A, square or not, singular or not, has exactly one matrix A⁺ satisfying the four Penrose conditions:

A·A⁺·A = A          A⁺·A·A⁺ = A⁺
(A·A⁺)ᵀ = A·A⁺      (A⁺·A)ᵀ = A⁺·A

Reading: A⁺ undoes A as far as A can be undone, wastes nothing doing it, and the two products are orthogonal projections — onto the column space and the row space respectively. Uniqueness is the point: there is one pseudoinverse, not a family. Written here for real matrices; for complex data — FRFs, phasors, S-parameters — replace every ᵀ with the conjugate transpose ᴴ, and everything below still holds ^[1].

The workhorse construction goes through the singular value decomposition. If A = U·Σ·Vᵀ, then:

A⁺ = V·Σ⁺·Uᵀ        where Σ⁺ flips each nonzero σᵢ to 1/σᵢ and leaves zeros alone

Reading: invert what's invertible, project away what isn't. This is how pinv actually computes it.

Two special cases cover most bench work. Tall matrix (more equations than unknowns), columns independent:

A⁺ = (Aᵀ·A)⁻¹·Aᵀ        x̂ = A⁺·b  minimizes  ‖A·x − b‖

Reading: the least-squares fit — the x̂ that makes the residual as small as the 2-norm allows.

Wide matrix (more unknowns than equations), rows independent:

A⁺ = Aᵀ·(A·Aᵀ)⁻¹        x̂ = A⁺·b  is the smallest x with A·x = b exactly

Reading: infinitely many solutions exist; this picks the one of minimum norm — no motion, current, or coefficient the data didn't ask for.

Where you meet it

Calibrating a multi-axis load cell on the test stand. A six-axis balance has a 6×6 sensitivity matrix, and you never characterize it with exactly six loadings — you apply thirty or forty known load cases and record all the channels. That's a tall system C·L = V, and the calibration matrix falls out of a pseudoinverse. Same story for accelerometer triads, magnetometer cal, and thermocouple polynomial fits: more measurements than coefficients, always.
Driving a robot arm through inverse kinematics. Resolved-rate control turns a commanded tool velocity into joint rates through q̇ = J⁺·ẋ. A 7-joint arm chasing a 6-DOF tool pose is underdetermined — the pseudoinverse picks the joint motion of minimum norm, which is why the elbow doesn't flail while the wrist does the work.
At the data review after a modal survey. Extracting mode shapes from a forest of FRF measurements, or fitting a state-space model to test data, is one overdetermined solve after another. When the reviewer asks "how did you fit this?", the honest answer usually contains a pseudoinverse.
Adjusting a survey or antenna-array network with a floating datum. Bjerhammar built the matrix version of this for geodesy: a network measured only by relative observations has a rank-deficient normal matrix, and the pseudoinverse supplies the free-network adjustment nobody's arbitrary datum choice can bias.

How it works

Picture the two failure modes of A·x = b. Overdetermined: b isn't in the column space of A, so no exact solution exists — A⁺·b gives the x whose image lands on the closest point that is, the orthogonal projection. Underdetermined: a whole affine family of solutions exists — A⁺·b gives the member closest to the origin. Both behaviors come from the same object, and Penrose proved in 1956 that A⁺·b is the minimum-norm least-squares solution in every case at once ^[7]. Small check you can do in your head: the single equation x + 2y = 5 has a line of solutions; A⁺·b = (1, 2), the point on that line nearest the origin. And fitting a line through (0, 1), (1, 1.9), (2, 3.1) gives intercept 0.95, slope 1.05 — identical to what any least-squares routine returns, because it is the same computation.

The SVD form tells you everything about when to trust it. Each singular value gets flipped to its reciprocal, so a singular value of 10⁻⁸ becomes a gain of 10⁸ on whatever noise lives in that direction. Three consequences, in the order they bite:

Near-singular is worse than singular. An exactly-zero singular value gets projected away cleanly; a tiny one gets amplified enormously. This is why library routines take a cutoff — NumPy's pinv(A, rcond=...) treats singular values below a threshold as zero ^[10]. Where you set that threshold is an engineering judgment about your noise floor, not a math question.
The pseudoinverse is discontinuous in rank. pinv(diag(1, ε)) has an entry 1/ε — 100 at ε = 10⁻², 10⁸ at ε = 10⁻⁸ — then snaps to 0 at ε = 0 exactly ^[1]. A robot arm approaching a kinematic singularity rides this cliff: joint-rate commands blow up just before the Jacobian loses rank. The standard fix is damped least squares, Jᵀ·(J·Jᵀ + λ²·I)⁻¹, which trades a little tracking error for bounded joint rates.
Don't form Aᵀ·A when the problem is ill-conditioned. The normal-equations shortcut squares the condition number — a matrix with κ = 10⁵ becomes a solve at κ = 10¹⁰, and in double precision you've spent half your digits before starting. QR or SVD works on A directly and doesn't pay that tax ^[1].

The mistake people make is treating A⁺ as if it were A⁻¹ with better manners. It isn't. A⁺·A ≠ I unless the columns are independent, and A·A⁺ ≠ I unless the rows are; (A·B)⁺ ≠ B⁺·A⁺ in general ^[9], which trips up everyone who chains transformations symbolically. And "minimum norm" is a mathematical preference, not a physical one — the smallest-norm joint motion or current distribution is only the right answer if all your unknowns share units and comparable scales. Mixed units, mixed channel noise? Scale or weight the system first, or the pseudoinverse quietly optimizes something you never meant to ask for.

If A happens to be square and nonsingular, A⁺ = A⁻¹ exactly. The pseudoinverse never contradicts the inverse; it extends it to the 99 percent of test matrices that were never square to begin with.

History

The idea is exactly as old as it sounds and got ignored for exactly the reason you'd guess. E. H. Moore — head of mathematics at the University of Chicago from the school's opening — worked out what he called the "general reciprocal" of an arbitrary matrix during the 1910s and announced it in a two-page abstract in the Bulletin of the American Mathematical Society in 1920 ^[2]^[3]. The full treatment waited until General Analysis, published in 1935, three years after his death ^[2]. Almost nobody read it: Moore wrote in a private symbolic notation so idiosyncratic that a centennial review politely called it "illegible for all but very dedicated readers" ^[6]. The result sat unused for three decades.

It was reinvented twice. Arne Bjerhammar, a Swedish geodesist at the Royal Institute of Technology in Stockholm, built reciprocals of rectangular matrices in 1951 to run least-squares adjustments of survey networks — an engineer deriving the tool because the job demanded it ^[1]^[8]. Then in 1955 Roger Penrose, at the time a Cambridge graduate student in his early twenties, published "A generalized inverse for matrices," defining the inverse by the four equations above and proving there is exactly one solution for any matrix ^[4]^[5]. The clean axiomatic framing is why the subject finally caught fire. A year later Penrose delivered the payoff engineers actually use — the proof that this inverse produces best approximate solutions of linear systems ^[7] — and, in the same 1956 volume of the same journal, Richard Rado published a note connecting Penrose's inverse back to Moore's forgotten general reciprocal ^[6]^[11]. The hyphenated name stuck. Moore had been dead twenty-three years when it did.

Related tools

/tools/strain-gauge-bridge — one channel of the multi-channel calibration whose full sensitivity matrix a pseudoinverse extracts
/tools/loop-4-20ma — the two-point exactly-determined scaling; add more calibration points than coefficients and you're in pseudoinverse territory
/tools/snr-enob — the acquisition noise floor that decides where to set the singular-value cutoff
/tools/wheatstone-bridge — the front end producing the voltages your overdetermined calibration solve consumes

Sources

Written by HE in our own words from the cited sources — engineering judgment included, your stamp still required. All entries →