The Moore-Penrose Pseudoinverse

Ordinary matrix inverses only exist for square, full-rank matrices, but real data almost never comes that neatly shaped. A regression with 500 daily returns and 5 factors has a $500 \times 5$ design matrix — not square, so it has no ordinary inverse at all. The Moore-Penrose pseudoinverse is the fix: a generalized inverse that exists for any matrix, square or not, full rank or not, and that produces the best possible approximate answer when an exact one doesn't exist.

An analogy: fitting a line through scattered points

You can't draw a single straight line that passes exactly through every point of a noisy scatter of stock returns against a factor — there are more data points (equations) than there is freedom in the line (two unknowns: slope and intercept), so no exact solution exists. But you can find the line that comes closest, in a precise least-squares sense, minimizing the total squared vertical distance from the points to the line. The pseudoinverse is the machinery that produces exactly this "closest possible" answer whenever a system is over- or under-determined, rather than giving up because no exact solution exists.

The mathematics, one symbol at a time

For a matrix $A$ that is tall (more rows than columns, an over-determined system) and full column rank, the pseudoinverse is

$$A^{+} = (A^T A)^{-1} A^T,$$

read as: transpose $A$, multiply it by $A$ to get a square matrix that (unlike $A$ itself) usually is invertible, invert that, then multiply by $A^T$ again. Applying $A^+$ to a target vector $b$ that has no exact solution to $Ax=b$ gives $x^* = A^+ b$, the vector $x^*$ that minimizes $\|Ax - b\|^2$ — the total squared error between what $Ax^*$ actually produces and the target $b$ you wanted. In words: when no $x$ can hit $b$ exactly, $A^+ b$ finds the $x$ that gets closest, in the sense of smallest total squared miss. For a square, invertible $A$, the pseudoinverse reduces exactly to the ordinary inverse $A^{-1}$ — it's a strict generalization, not a different concept.

Worked example 1: least-squares regression by hand

Three data points relate a factor $f$ to a return $r$: $(f,r) = (1, 2), (2, 3), (3, 5)$, and we fit $r = \beta f$ (no intercept, for simplicity). Here $A = \begin{pmatrix}1\\2\\3\end{pmatrix}$, $b=\begin{pmatrix}2\\3\\5\end{pmatrix}$. Then $A^TA = 1^2+2^2+3^2=14$, and $A^Tb = 1(2)+2(3)+3(5)=2+6+15=23$. So $\beta = A^+b = (A^TA)^{-1}A^Tb = 23/14 \approx 1.643$. This is exactly the ordinary least-squares slope — the pseudoinverse and OLS are the same computation. Checking fit: predicted values are $1.64, 3.29, 4.93$ against actual $2, 3, 5$ — close but not exact, exactly as expected since three points can't lie perfectly on one line through the origin.

Worked example 2: an under-determined system

Now the opposite case: fewer equations than unknowns, e.g. one equation $x_1 + x_2 = 4$ (a single constraint, two unknowns) — infinitely many solutions exist ($x_1=4,x_2=0$; $x_1=0,x_2=4$; and everywhere between). The pseudoinverse picks out the minimum-norm solution: the one closest to the origin among all valid solutions, which here is $x_1=x_2=2$ (by symmetry, splitting the constraint evenly minimizes $x_1^2+x_2^2$). This "smallest, most conservative solution among the valid ones" behavior is why the pseudoinverse is also used to regularize under-determined systems, like allocating weight across many correlated factors with no unique answer.

What this means in practice

Whenever numpy.linalg.lstsq or a regression library fits coefficients, it's computing (a numerically stable version of) the pseudoinverse under the hood — this is the actual mechanism behind ordinary least squares. It's essential when the design matrix isn't full rank (near-duplicate factors) or isn't square (more or fewer observations than parameters), both extremely common situations in factor models and signal construction.