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

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.