Low-Rank Approximation (I): Pseudo-Inverse
Perhaps many readers, like myself, have a feeling toward low-rank matrix approximation that is both familiar and unfamiliar. The familiarity comes from the fact that the concept and significance of low-rank approximation are not hard to grasp, and with fine-tuning techniques such as LoRA — which are built on low-rank approximation — now flourishing everywhere, the idea has already become deeply ingrained in our minds through sheer exposure. However, the material covered under "low-rank approximation" is extremely broad, and in papers on the topic one often runs into unfamiliar yet impressive tricks, which leaves us with a sense of half-understanding, half-unfamiliarity.
For this reason, in this series I will try to systematically organize the theoretical content related to low-rank matrix approximation, so as to fill in the gaps in our understanding. In this first post, we will introduce a relatively simple concept within the low-rank approximation family — the pseudo-inverse.
An Optimization Perspective
The pseudo-inverse, also called the "generalized inverse," is, as the name suggests, a "generalized inverse matrix." It is in fact a generalization of the notion of an "inverse matrix" to matrices that are not invertible. more
We know that for a matrix equation $\boldsymbol{A}\boldsymbol{B}=\boldsymbol{M}$, if $\boldsymbol{A}$ is square and invertible, then we can directly obtain $\boldsymbol{B}=\boldsymbol{A}^{-1}\boldsymbol{M}$. But what if $\boldsymbol{A}$ is not invertible, or is not even square? In that case we may well be unable to find any $\boldsymbol{B}$ satisfying $\boldsymbol{A}\boldsymbol{B}=\boldsymbol{M}$, and if we still want to make progress, the usual move is to convert this into an optimization problem:
\begin{equation}\mathop{\text{argmin}}_\boldsymbol{B} \Vert \boldsymbol{A}\boldsymbol{B} - \boldsymbol{M}\Vert_F^2\label{eq:loss-ab-m}\end{equation}
where $\boldsymbol{A}\in\mathbb{R}^{n\times r}, \boldsymbol{B}\in\mathbb{R}^{r\times m}, \boldsymbol{M}\in\mathbb{R}^{n\times m}$ — note that the field here is $\mathbb{R}$, indicating that this series will focus on real matrices — and $\Vert\cdot\Vert_F$ is the $F$ norm of the matrix (Frobenius norm), used to measure the distance between the matrix $\boldsymbol{A}\boldsymbol{B} - \boldsymbol{M}$ and the all-zero matrix, defined as
\begin{equation}\Vert \boldsymbol{M}\Vert_F = \sqrt{\sum_{i=1}^n\sum_{j=1}^m M_{i,j}^2}\end{equation}
In plain terms, this changes the task from finding an exact inverse matrix to minimizing the squared error between $\boldsymbol{A}\boldsymbol{B}$ and $\boldsymbol{M}$. Since the theme of this series is low-rank approximation, in what follows we assume $r \ll \min(n,m)$; the machine learning interpretation is that we are reconstructing the full matrix $\boldsymbol{M}$ from a low-dimensional, lossy input matrix $\boldsymbol{A}$ together with a linear transformation $\boldsymbol{B}$.
When $m=n$ and we take $\boldsymbol{M}$ to be the identity matrix $\boldsymbol{I}_n$, we obtain a result that depends only on $\boldsymbol{A}$, which we denote as
\begin{equation}\boldsymbol{A}^{\dagger} = \mathop{\text{argmin}}_\boldsymbol{B} \Vert \boldsymbol{A}\boldsymbol{B} - \boldsymbol{I}_n\Vert_F^2\label{eq:loss-ab-m-b}\end{equation}
Its role is analogous to the inverse matrix of $\boldsymbol{A}$, so it is called the "(right) pseudo-inverse" of $\boldsymbol{A}$. Similarly, if what is given is a $\boldsymbol{B}$ matrix, we can instead treat $\boldsymbol{A}$ as the optimization variable, obtaining the "(left) pseudo-inverse" of $\boldsymbol{B}$:
\begin{equation}\boldsymbol{B}^{\dagger} = \mathop{\text{argmin}}_\boldsymbol{A} \Vert \boldsymbol{A}\boldsymbol{B} - \boldsymbol{I}_n\Vert_F^2\end{equation}
Norms
Before proceeding further, let us fill in some background on the $F$ norm. The norm of a vector is presumably already familiar to many readers, the classic example being the "$p$-norm": for $\boldsymbol{x}=(x_1,x_2,\cdots,x_m)$, its $p$-norm is defined as
\begin{equation}\Vert \boldsymbol{x}\Vert_p = \sqrt[\uproot{10}p]{\sum_{i=1}^m |x_i|^p}\end{equation}
Among the $p$-norms, the most common case is $p=2$, which is simply what we usually call the vector's magnitude, also known as the "Euclidean norm." If we omit the subscript and just write $\Vert \boldsymbol{x}\Vert$, it basically defaults to $p=2$.
The norm of a matrix is a bit more complex — there are at least two different but commonly used norms. One of them is the $F$ norm already mentioned in the previous section, which is computed by simply flattening the matrix into a vector:
\begin{equation}\Vert \boldsymbol{M}\Vert_F = \Vert \text{vec}(\boldsymbol{M})\Vert_2 = \sqrt{\sum_{i=1}^n\sum_{j=1}^m M_{i,j}^2}\end{equation}
We'll introduce other matrix norms as we encounter them. Because there is more than one kind of matrix norm, the subscript ${}_F$ in $\Vert \boldsymbol{M}\Vert_F$ usually cannot be omitted, to avoid confusion. The $F$ norm comes from treating the matrix as a vector and simply borrowing the vector-norm definition, which inspires us to try carrying more vector operations over to matrices — for instance, the inner product:
\begin{equation}\langle \boldsymbol{P}, \boldsymbol{Q} \rangle_F = \langle \text{vec}(\boldsymbol{P}), \text{vec}(\boldsymbol{Q}) \rangle = \sum_{i=1}^n\sum_{j=1}^m P_{i,j} Q_{i,j}\end{equation}
This is called the $F$ inner product (Frobenius inner product) of matrices $\boldsymbol{P},\boldsymbol{Q}$, where $\boldsymbol{P},\boldsymbol{Q}\in\mathbb{R}^{n\times m}$, and it can be expressed using the trace operation on matrices:
\begin{equation}\langle \boldsymbol{P}, \boldsymbol{Q} \rangle_F = \text{Tr}(\boldsymbol{P}^{\top} \boldsymbol{Q})\end{equation}
This can be proved directly from the definitions of matrix multiplication and the trace (readers are encouraged to try it themselves). When $\boldsymbol{P},\boldsymbol{Q}$ is a product of several matrices multiplied together, converting it into the equivalent trace form often helps simplify things. For example, we can use this to prove that orthogonal transformations preserve the $F$ norm: suppose $\boldsymbol{U}$ is an orthogonal matrix; using $\Vert \boldsymbol{M}\Vert_F^2 = \langle \boldsymbol{M}, \boldsymbol{M} \rangle_F$ together with the relation between the $F$ inner product and the trace, we obtain
\begin{equation}\Vert \boldsymbol{U}\boldsymbol{M}\Vert_F^2 = \text{Tr}((\boldsymbol{U}\boldsymbol{M})^{\top} \boldsymbol{U}\boldsymbol{M})= \text{Tr}(\boldsymbol{M}^{\top} \boldsymbol{U}^{\top}\boldsymbol{U}\boldsymbol{M})=\text{Tr}(\boldsymbol{M}^{\top} \boldsymbol{M}) = \Vert \boldsymbol{M}\Vert_F^2\end{equation}
Matrix Derivatives
Back to the main thread: for an optimization objective, the ideal outcome is naturally to obtain an analytical solution via differentiation, and $\eqref{eq:loss-ab-m}$ happens to allow exactly that! This conclusion can be "eyeballed" directly: $\boldsymbol{A}\boldsymbol{B}-\boldsymbol{M}$ is a linear function of $\boldsymbol{B}$, so $\Vert \boldsymbol{A}\boldsymbol{B}-\boldsymbol{M}\Vert_F^2$ is a quadratic function of $\boldsymbol{A}$ or $\boldsymbol{B}$, and a quadratic function's minimum has an analytical solution.
To find the derivative of $\mathcal{L}=\Vert \boldsymbol{A}\boldsymbol{B}-\boldsymbol{M}\Vert_F^2$ with respect to $\boldsymbol{B}$, we first need the derivative of $\mathcal{L}$ with respect to $\boldsymbol{E}=\boldsymbol{A}\boldsymbol{B}-\boldsymbol{M}$, then the derivative of $\boldsymbol{E}$ with respect to $\boldsymbol{B}$, and finally combine them via the chain rule, i.e.
\begin{equation}\frac{\partial \mathcal{L}}{\partial B_{i,j}} = \sum_{k,l}\frac{\partial \mathcal{L}}{\partial E_{k,l}}\frac{\partial E_{k,l}}{\partial B_{i,j}} \end{equation}
By definition $\mathcal{L}=\Vert \boldsymbol{E}\Vert_F^2 = \sum_{i,j} E_{i,j}^2$, it's clear that among all the squared terms in the sum, only when $(i,j)=(k,l)$ does the derivative with respect to $E_{k,l}$ fail to vanish, so the derivative of $\mathcal{L}$ with respect to $E_{k,l}$ equals the derivative of $E_{k,l}^2$ with respect to $E_{k,l}$, namely $\frac{\partial \mathcal{L}}{\partial E_{k,l}}=2E_{k,l}$. Next, from the definition of matrix multiplication we have
\begin{equation}E_{k,l} = \left(\sum_{\alpha} A_{k,\alpha}B_{\alpha,l}\right) - M_{k,l}\end{equation}
Similarly, only when $(\alpha,l)=(i,j)$ does the derivative of the above expression with respect to $B_{i,j}$ produce a nonzero result $A_{k,i}$, so we can write $\frac{\partial E_{k,l}}{\partial B_{i,j}}=A_{k,i}\delta_{l,j}$, where $\delta_{l,j}$ is the Kronecker delta, used to encode the condition $l=j$. Putting the pieces together, we get
\begin{equation}\frac{\partial \mathcal{L}}{\partial B_{i,j}} = 2\sum_{k,l}E_{k,l}A_{k,i}\delta_{l,j} = 2\sum_k E_{k,j}A_{k,i}\end{equation}
If we adopt the convention that the gradient of a scalar with respect to a matrix has the same shape as the matrix itself, then we can write
\begin{equation}\frac{\partial \mathcal{L}}{\partial \boldsymbol{B}} = 2\boldsymbol{A}^{\top}(\boldsymbol{A}\boldsymbol{B}-\boldsymbol{M})\end{equation}
Although the derivation took some effort, the result is fortunately quite intuitive: intuitively, $\frac{\partial \mathcal{L}}{\partial \boldsymbol{B}}$ should be the product of $2(\boldsymbol{A}\boldsymbol{B}-\boldsymbol{M})$ and $\boldsymbol{A}$ (by analogy with scalar differentiation), and since we've agreed that the shape of $\frac{\partial \mathcal{L}}{\partial \boldsymbol{B}}$ matches that of $\boldsymbol{B}$ (i.e., $r\times m$), we need to find a way to multiply $2(\boldsymbol{A}\boldsymbol{B}-\boldsymbol{M})\in\mathbb{R}^{n\times m}$ and $\boldsymbol{A}\in\mathbb{R}^{n\times r}$ together to produce a result of shape $r\times m$ — and the only way to do so is the expression on the right-hand side above. By the same reasoning, we can quickly write down
\begin{equation}\frac{\partial \mathcal{L}}{\partial \boldsymbol{A}} = 2(\boldsymbol{A}\boldsymbol{B}-\boldsymbol{M})\boldsymbol{B}^{\top}\end{equation}
Basic Result
Now that we've derived $\frac{\partial \mathcal{L}}{\partial \boldsymbol{B}}$ and $\frac{\partial \mathcal{L}}{\partial \boldsymbol{A}}$ separately, setting them to zero lets us solve for the corresponding optimal solution:
\begin{align} 2\boldsymbol{A}^{\top}(\boldsymbol{A}\boldsymbol{B}-\boldsymbol{M})=0\quad\Rightarrow\quad (\boldsymbol{A}^{\top} \boldsymbol{A})^{-1}\boldsymbol{A}^{\top}\boldsymbol{M} = \mathop{\text{argmin}}_\boldsymbol{B} \Vert \boldsymbol{A}\boldsymbol{B} - \boldsymbol{M}\Vert_F^2 \\ 2(\boldsymbol{A}\boldsymbol{B}-\boldsymbol{M})\boldsymbol{B}^{\top}=0\quad\Rightarrow\quad \boldsymbol{M}\boldsymbol{B}^{\top}(\boldsymbol{B} \boldsymbol{B}^{\top})^{-1} = \mathop{\text{argmin}}_\boldsymbol{A} \Vert \boldsymbol{A}\boldsymbol{B} - \boldsymbol{M}\Vert_F^2 \end{align}
Substituting into $\boldsymbol{M}=\boldsymbol{I}_n$, we get
\begin{align}\boldsymbol{A}^{\dagger} = (\boldsymbol{A}^{\top} \boldsymbol{A})^{-1}\boldsymbol{A}^{\top} \label{eq:p-inv-a}\\ \boldsymbol{B}^{\dagger} = \boldsymbol{B}^{\top}(\boldsymbol{B} \boldsymbol{B}^{\top})^{-1}\label{eq:p-inv-b}\end{align}
If $\boldsymbol{A}$ or $\boldsymbol{B}$ is an invertible square matrix, it is easy to show that the pseudo-inverse equals the ordinary inverse, i.e., $\boldsymbol{A}^{\dagger}=\boldsymbol{A}^{-1},\boldsymbol{B}^{\dagger}=\boldsymbol{B}^{-1}$. Furthermore, from the above we can also verify:
1. $(\boldsymbol{A}^{\dagger})^{\dagger}=\boldsymbol{A},(\boldsymbol{B}^{\dagger})^{\dagger}=\boldsymbol{B}$, i.e., the pseudo-inverse of the pseudo-inverse equals itself. This means that while serving as an approximate inverse, the pseudo-inverse also fully preserves its own information.
2. $\boldsymbol{A}\boldsymbol{A}^{\dagger}\boldsymbol{A}=\boldsymbol{A},\boldsymbol{B}^{\dagger}\boldsymbol{B}\boldsymbol{B}^{\dagger}=\boldsymbol{B}^{\dagger}$, i.e., although $\boldsymbol{A}\boldsymbol{A}^{\dagger},\boldsymbol{B}^{\dagger}\boldsymbol{B}$ cannot become the identity matrix $\boldsymbol{I}$, it plays the role of an identity matrix as far as $\boldsymbol{A},\boldsymbol{B}^{\dagger}$ is concerned.
Incidentally, the pseudo-inverse of a matrix is actually a very broad concept with many different forms. What we've introduced here is in fact the most common one, the "Moore–Penrose inverse"; besides this there is also the "Drazin inverse," the "Bott–Duffin inverse," and others, but I am not familiar with these, so I won't go into them — interested readers can consult the Wikipedia entry on "generalized inverse."
The General Case
But we're not done yet. The derivation of $\eqref{eq:p-inv-a},\eqref{eq:p-inv-b}$ hinges on a crucial assumption: that the corresponding $\boldsymbol{A}^{\top} \boldsymbol{A}$ or $\boldsymbol{B} \boldsymbol{B}^{\top}$ is invertible. What if it isn't?
Let's take $\boldsymbol{A}^{\dagger}$ as an example, and suppose $\boldsymbol{A}^{\top} \boldsymbol{A}$ is not invertible. This means the rank of $\boldsymbol{A}$ is less than $r$, and the best we can do is find a maximal linearly independent set of $s < r$ column vectors forming a matrix $\boldsymbol{A}_s\in\mathbb{R}^{n\times s}$, after which $\boldsymbol{A}$ can be expressed as $\boldsymbol{A}_s \boldsymbol{P}$, where $\boldsymbol{P}\in\mathbb{R}^{s\times r}$ maps from $\boldsymbol{A}_s$ to $\boldsymbol{A}$. In this case,
\begin{equation}\mathop{\text{argmin}}_\boldsymbol{B} \Vert \boldsymbol{A}\boldsymbol{B} - \boldsymbol{I}_n\Vert_F^2 = \mathop{\text{argmin}}_\boldsymbol{B} \Vert \boldsymbol{A}_s \boldsymbol{P}\boldsymbol{B} - \boldsymbol{I}_n\Vert_F^2\end{equation}
If we still denote the optimal solution of $\boldsymbol{B}$ as $\boldsymbol{A}^{\dagger}$, then all we can determine is
\begin{equation}\boldsymbol{P}\boldsymbol{A}^{\dagger} = \boldsymbol{A}_s^{\dagger} = (\boldsymbol{A}_s^{\top} \boldsymbol{A}_s)^{-1}\boldsymbol{A}_s^{\top}\end{equation}
Since we've already assumed that $\boldsymbol{A}_s$ is a maximal linearly independent set, $\boldsymbol{A}_s^{\top} \boldsymbol{A}_s$ must be invertible, so the above expression is well-defined. However, going from $\boldsymbol{A}^{\dagger}$ to $\boldsymbol{P}\boldsymbol{A}^{\dagger}$ is a dimension-reducing process, which means there exist multiple $\boldsymbol{A}^{\dagger}$ such that $\boldsymbol{P}\boldsymbol{A}^{\dagger} = \boldsymbol{A}_s^{\dagger}$ — that is, at this point the optimal solution to the objective $\eqref{eq:loss-ab-m-b}$ is not unique. In other words, when $\boldsymbol{A}^{\top} \boldsymbol{A}$ is not invertible, we cannot pin down a unique pseudo-inverse $\boldsymbol{A}^{\dagger}$ from the objective $\eqref{eq:loss-ab-m-b}$ alone.
One possible approach is to add extra constraints on $(\boldsymbol{A}^{\dagger})^{\dagger}=\boldsymbol{A}$ or $\boldsymbol{A}^{\dagger}\boldsymbol{A}\boldsymbol{A}^{\dagger}=\boldsymbol{A}^{\dagger}$, so that combined with $\boldsymbol{P}\boldsymbol{A}^{\dagger} = \boldsymbol{A}_s^{\dagger}$ this would uniquely determine $\boldsymbol{A}^{\dagger}$. However, this patchwork approach feels a bit too ad hoc; in fact, we can handle this problem in a more elegant and unified way using a neat trick. The issue arises because when $\boldsymbol{A}^{\top} \boldsymbol{A}$ is not invertible, the optimal solution of the objective function $\eqref{eq:loss-ab-m-b}$ is not unique. We can add a regularization term to make it unique, solve for the result, and then let the weight of the regularization term tend to zero:
\begin{equation}\boldsymbol{A}^{\dagger} = \lim_{\epsilon\to 0}\,\mathop{\text{argmin}}_\boldsymbol{B} \Vert \boldsymbol{A}\boldsymbol{B} - \boldsymbol{I}_n\Vert_F^2 + \epsilon\Vert \boldsymbol{B}\Vert_F^2\end{equation}
Here $\epsilon > 0$, and $\epsilon\to 0$ denotes tending to zero from the positive side. From the above we can solve for
\begin{equation}\boldsymbol{A}^{\dagger} = \lim_{\epsilon\to 0}\,(\boldsymbol{A}^{\top} \boldsymbol{A} + \epsilon \boldsymbol{I}_r)^{-1}\boldsymbol{A}^{\top}\label{eq:a-pinv-lim}\end{equation}
When $\epsilon > 0$, $\boldsymbol{A}^{\top} \boldsymbol{A} + \epsilon \boldsymbol{I}_r$ must be invertible (readers are invited to prove this), so the above expression is well-defined. Since when $\epsilon\to 0$ the regularization term becomes negligible, the limit above must exist. Note that we're talking about the existence of the overall limit here: when $\boldsymbol{A}^{\top} \boldsymbol{A}$ is not invertible, the limit $\lim\limits_{\epsilon\to 0}\,(\boldsymbol{A}^{\top} \boldsymbol{A} + \epsilon \boldsymbol{I}_r)^{-1}$ itself does not exist (the result would blow up to infinity) — it's only after multiplying by $\boldsymbol{A}^{\top}$ and then taking the overall limit that we get a well-behaved result.
What advantages does formula $\eqref{eq:a-pinv-lim}$, as a general extension of the pseudo-inverse, offer? First, we already have the expression $\eqref{eq:p-inv-a}$ for $\boldsymbol{A}^{\dagger}$ when $\boldsymbol{A}^{\top} \boldsymbol{A}$ is invertible, and formula $\eqref{eq:a-pinv-lim}$, as its generalization, has the theoretical elegance of being both intuitive and formally consistent with it. Second, this formal consistency also means that the properties of $\boldsymbol{A}^{\dagger}$ — such as $(\boldsymbol{A}^{\dagger})^{\dagger}$ — that hold when $\boldsymbol{A}^{\top} \boldsymbol{A}$ is invertible continue to be preserved, so that when discussing $\boldsymbol{A}^{\dagger}$ we can almost entirely set aside the question of whether $\boldsymbol{A}^{\top} \boldsymbol{A}$ is invertible.
Numerical Computation
Of course, formula $\eqref{eq:a-pinv-lim}$ as it stands is only a formal definition; if we were to use it directly for numerical computation, we would need to take a sufficiently small $\epsilon$ and compute $(\boldsymbol{A}^{\top} \boldsymbol{A} + \epsilon \boldsymbol{I}_r)^{-1}$ accordingly, which would inevitably run into severe numerical instability. To arrive at a stable way of computing this, we make use of the fact that real symmetric matrices can always be orthogonally diagonalized (the spectral theorem), and decompose $\boldsymbol{A}^{\top} \boldsymbol{A}$ as follows:
\begin{equation}\boldsymbol{A}^{\top} \boldsymbol{A} = \boldsymbol{U}\boldsymbol{\Lambda} \boldsymbol{U}^{\top}\end{equation}
where $\boldsymbol{U}$ is an orthogonal matrix and $\boldsymbol{\Lambda}=\text{diag}(\lambda_1,\lambda_2,\cdots,\lambda_r)$ is a diagonal matrix of eigenvalues; since $\boldsymbol{A}^{\top} \boldsymbol{A}$ is positive semi-definite, its eigenvalues are always non-negative. Using this decomposition, we have
\begin{equation}\begin{aligned} (\boldsymbol{A}^{\top} \boldsymbol{A} + \epsilon \boldsymbol{I}_r)^{-1}\boldsymbol{A}^{\top} =&\, (\boldsymbol{U}\boldsymbol{\Lambda} \boldsymbol{U}^{\top} + \epsilon \boldsymbol{I}_r)^{-1} \boldsymbol{A}^{\top} \\ =&\, [\boldsymbol{U}(\boldsymbol{\Lambda} + \epsilon \boldsymbol{I}_r) \boldsymbol{U}^{\top}]^{-1}\boldsymbol{A}^{\top} \\ =&\, \boldsymbol{U}(\boldsymbol{\Lambda} + \epsilon \boldsymbol{I}_r)^{-1} \boldsymbol{U}^{\top}\boldsymbol{A}^{\top} \end{aligned}\end{equation}
For $(\boldsymbol{\Lambda} + \epsilon \boldsymbol{I}_r)^{-1}$ we have
\begin{equation}(\boldsymbol{\Lambda} + \epsilon \boldsymbol{I}_r)^{-1} = \text{diag}\Big((\lambda_1 + \epsilon)^{-1},(\lambda_2 + \epsilon)^{-1},\cdots,(\lambda_r + \epsilon)^{-1}\Big)\end{equation}
If $\lambda_i > 0$, then $\lim\limits_{\epsilon\to 0}\,(\lambda_i + \epsilon)^{-1}=\lambda_i^{-1}$ is a finite result, which poses no obstacle to computation; the problem occurs when $\lambda_i = 0$, in which case $\lim\limits_{\epsilon\to 0}\,(\lambda_i + \epsilon)^{-1}=\lim\limits_{\epsilon\to 0}\,\epsilon^{-1}\to\infty$. However, we know that as $\epsilon\to 0$ the influence of the regularization term must vanish, so we can conclude that the limit $\eqref{eq:a-pinv-lim}$ must not blow up to infinity. Therefore, if $\lambda_i=0$ exists, the factor $\boldsymbol{U}^{\top}\boldsymbol{A}^{\top}$ multiplying it on the right-hand side must somehow cancel out the infinity coming from $\lim\limits_{\epsilon\to 0}\, \epsilon^{-1}$. And the only thing that can cancel such an infinity is "multiplying by 0," i.e., $\lim\limits_{\epsilon\to 0}\, \epsilon^{-1}\times 0 = 0$.
In other words, if $\lambda_i=0$, then the factor that $\boldsymbol{U}^{\top}\boldsymbol{A}^{\top}$ contributes to multiplying $(\lambda_i+\epsilon)^{-1}$ must be $0$. Given this, since "0 times anything is 0," it turns out that the actual value taken by $(\lambda_i+\epsilon)^{-1}$ when $\lambda_i=0$ doesn't actually matter — we can simply set it to 0. This gives us a general and simple method for computing $\boldsymbol{A}^{\dagger}$:
\begin{equation}\boldsymbol{A}^{\dagger} = \boldsymbol{U}\boldsymbol{\Lambda}^{\dagger}\boldsymbol{U}^{\top}\boldsymbol{A}^{\top}, \quad \boldsymbol{A}^{\top} \boldsymbol{A} = \boldsymbol{U}\boldsymbol{\Lambda} \boldsymbol{U}^{\top}\end{equation}
where $\boldsymbol{\Lambda}^{\dagger}$ means that diagonal entries equal to zero are left unchanged, while nonzero entries are replaced by their reciprocal.
Some readers might wonder: since "0 times anything is 0," why should the entries of $\lambda_i$ that equal zero be left unchanged? Couldn't we just pick some other arbitrary value? Actually, picking any other value here wouldn't change the final result either, but since we've adopted the notation $\boldsymbol{\Lambda}^{\dagger}$, we should keep it consistent with formula $\eqref{eq:a-pinv-lim}$, i.e., it should agree with the result of directly substituting the diagonal matrix $\boldsymbol{\Lambda}$ into formula $\eqref{eq:a-pinv-lim}$:
\begin{equation}\boldsymbol{\Lambda}^{\dagger} = \lim_{\epsilon\to 0}\,(\boldsymbol{\Lambda}^{\top} \boldsymbol{\Lambda} + \epsilon \boldsymbol{I}_r)^{-1}\boldsymbol{\Lambda}^{\top} = \text{diag}(\kappa_1,\kappa_2,\cdots,\kappa_n),\quad \kappa_i = \left\{\begin{aligned}\lambda_i^{-1}, \,\,\lambda_i\neq 0 \\ 0, \,\,\lambda_i=0 \end{aligned} \right. \end{equation}
Summary
In this post, we introduced the pseudo-inverse from the perspective of low-rank approximation — an extension of the concept of the inverse matrix to non-square or non-invertible square matrices, allowing us to more effectively analyze and solve general matrix equations.
Translated automatically with claude-sonnet-5; all equations are reproduced verbatim from the source. Copyright remains with the original author.