Moore-Penrose pseudoinverse solves the least squares problem (SVD framework)

Question

I am a computer science researcher who has to learn some numerical linear algebra for my work. I have been struggling with the SVD and Moore-Penrose pseudoinverse as of late. I am trying to solve some problems to get more comfortable with what should probably be routine manipulations.

First of all, I have gone through similar questions on Stack Exchange but I believe they were more general and are not equivalent. I am working in the framework where $A^{\dagger} = V\Lambda^{\dagger}U^T$. So, basically, I am using SVD's. The matrix $A$ of course is identified with $U\Lambda V^T$

---

Problem

Consider the matrix equation $Ax=y$, where $A\in R^{m\times n}$. The corresponding least squares problem is to find a least squares solution $x_{\text{LS}}$ that minimizes the Euclidean norm of the residual, i.e.,

$$\|Ax_{\text{LS}}-y\| = \min_{x \in \Bbb R^n} \|Ax-y\| = \min_{z \in \mbox{Ran}(A)}\|z-y\|$$

a) Show that $A^{\dagger}y$ is a least-squares solution and satisfies the normal equation $A^TAx=A^Ty$. Why is this solution special?

b) Show that $\ker(A^TA) = \ker(A)$.

c) Use the above results to deduce that $x \in \Bbb R^n$ is a least-squares solution if and only if it satisfies the normal equation.

---

Help on any or all of these parts is appreciated. I'd also appreciate links to relevant posts. Like I said, I've read similar questions but did not understand them as they were in a more general framework.

Edit: I have solved b). It didn't depend on a) as I had initially thought and it is pretty straightforward to solve, see eg. here: Prove that for a real matrix $A$, $\ker(A) = \ker(A^TA)$

Edit: I realize that part a) might be more involved than I had expected... Assuming parts a) and b), can someone help me with part c)?

Answer 1

Answer to a) and b) follows from Why does $\operatorname{null}(A) = \operatorname{null}(A^TA)$, intuitively? and Why does SVD provide the least squares and least norm solution to $Ax=b$? as poited out in comments.

Let $\tilde{f}$ being a fixed minimizer of $$Q(f)=\|Af-g\|_K^2,$$ in which $$ \qquad A\in \mathbb{R}^{m\times n},\quad f\in \mathbb{R}^{n}=H,\quad g\in \mathbb{R}^{m}=K.$$

Let $h\in H$, and note that $$Q(f+h)=Q(f)+\langle Af-g,Ah\rangle_K+\langle Ah,Af-g\rangle_K+\|Ah\|_K^2\,\quad \tag{1},$$ in which $\langle\cdot,\cdot\rangle$ denotes the inner product. In particular $$Q(f+h)=Q(f),\qquad f\in H, \quad h\in \ker(A).$$

This and $(1)$ means that $$Q(\tilde{f})\leq Q(\tilde{f}+h)=Q(\tilde{f})+2\langle A^T(A\tilde{f}-g),h\rangle_H+\|Ah\|_K^2,\qquad h\in H.$$

Since $\tilde{f}$ is a critical point of $Q$, it follows that $$\nabla Q(\tilde{f})=2A^T(A\tilde{f}-g)=0,$$ and $(A\tilde{f}-g)\in \ker(A^*)=range(A)^\perp$.

You can find more details, in a more general case, in Theorem 1.1 of The mathematics of computerized tomography. You can also find related results searching for "\(f=A^+g\)" on SearchOnMath.

Answer 2

Here we show that $x=A^+b$, where $A^+$ is the Moore-Penrose pseudo inverse $A^+$ of a matrix $A$ solves the ordinary least square problem (OLS): $$x=\operatorname{arg}.\min\|Ax-b\|_2$$ where $A\in M(\mathbb{C},m, n)$ and $b\in\mathbb{C}^n$, and that if $x'$ is another solution, then $\|A^+b\|_2\leq\|x'\|_2$ (This is one reason that makes $A^+$ special).

Notice that if $b\in \operatorname{range}(A)$, then $x_*=A^+b$ is a solution to $Ax=b$, and has the smallest $\|\;\|_2$-norm.

---

First some linear algebra facts:

Using simple linear algebra, we now that if $R$ is a linear subspace of $\mathbb{C}^n$, then

• (A) The orthogonal projection $P_R:\mathbb{C}^n\rightarrow R$ defined for each $\mathbf{y}\in \mathbb{C}^n$ as a vector $P_r\mathbf{y}\in R$ such that $$\|P_R\mathbf{y}-\mathbf{y}\|_2\leq \|\mathbf{z}-\mathbf{y}\|_2,\qquad \mathbf{z}\in R$$ exists and is unique.
• (B) Furthermore, for each $\mathbf{y}\in\mathbb{C}^n$, $P_R\mathbf{y}$ is the unique vector that satisfies $$\mathbf{y}- P_R\mathbf{y}\perp R$$ that is, $\langle \mathbf{y}- P_R\mathbf{y},\mathbf{z}\rangle=0$ for all $\mathbf{z}\in R$ where $\langle\cdot,\cdot\rangle$ is the standard inner product in $\mathbb{C}^n$ (hence the name orthogonal projection).

---

Existence of solutions to the OLS

Consider the linear space $R:=\{A\mathbf{w}:\mathbf{w}\in\mathbb{C}^m\}$, that is, the range of $A$. The geometrical digression made above shows that there is exist a unique vector $\hat{\mathbf{y}}\in R$ such that $$\widehat{\mathbf{y}}=\operatorname{arg.min}_{\mathbf{z}\in R}\|\mathbf{y}-\mathbf{z}\|_2$$ namely the orthogonal projection of $\mathbf{y}$ onto $R$, $\widehat{\mathbf{y}}=P_R\mathbf{y}$. In particular, this shows that there is at least one $\boldsymbol{\omega}\in\mathbb{C}^m$ such that $\widehat{\mathbf{y}}=A\boldsymbol{\omega}$. In fact, for any solution $\mathbf{w}_H$ to the homogeneous equation $A\mathbf{w}=\mathbf{0}$ we have that $\widehat{\mathbf{y}}=A(\boldsymbol{\omega}+\mathbf{w}_H)$.

The problem is now reduced to finding an $\boldsymbol{\omega}\in \mathbb{C}^m$ for which $$\begin{align} \widehat{\mathbf{y}}=A\boldsymbol{\omega}\tag{1}\label{one} \end{align}$$

Now, as the range of the linear transformation $\mathbf{w}\mapsto A\mathbf{w}$ is generated by the columns of $A$, statement (B) shows that $\widehat{\mathbf{y}}=A\boldsymbol{\omega}$ is the unique vector in $\mathbb{C}^n$ satisfying $$\begin{align} A^*(\mathbf{y}-A\boldsymbol{\omega})=\mathbf{0}\tag{2}\label{two} \end{align}$$ where $A^*$ is the conjugate transpose of $A$ (when $A\in M(\mathbb{R},m.n)$ then $A^*=A^\intercal$). Equivalently, any $\boldsymbol{\omega}\in\mathbb{C}^m$ satisfying \eqref{one} (which we have already established exists) satisfies the equation $$\begin{align} A^* A\boldsymbol{\omega}=A^*\mathbf{y}\tag{2'}\label{twop} \end{align}$$

• If the rank of $A$ is $m$ (i.e. $n\geq m$ and the columns of $A$ are linearly independent), then $A^*A$ is a $m\times m$ invertible matrix and so, there is a unique solution $$\boldsymbol{\omega}=(A^* A)^{-1}A^*\mathbf{y}$$ Conversely, if $A^* A$ is invertible, then the rank of $A$ is $m$.

• In any other case, there will be infinitely many solutions to \eqref{two} (still, $\widehat{\mathbf{y}}=P_R\mathbf{y}$ is unique) since then $\operatorname{rank}(A^* A)<m$ and so, the homogeneous equation $A^* A\mathbf{w}=\mathbf{0}$ and thus, the homogeneous equation $A\mathbf{w}=\mathbf{0}$ will have infinitely many solutions.

Solutions to \eqref{twop} with different properties (minimal norm for example)can be obtained by appealing to generalized inverse matrices.

---

Minimal solution to the OLS problem and the Moore-Penrose pseudoinverse:

Recall that for any $A\in M(\mathbb{C},m, n)$, matrix $X\in M(\mathbb{C}, n, m)$ is a generalized inverse matrix of $A$ if $$AXA=X$$ Generalized inverse matrices exists and, when $A$ is invertible, there is only one such a matrix, namely $X=A^{-1}$.

Amongst all generalized inverse matrices os a matrix $A$, there is very special one, denoted by $A^+$, which satisfies $$\begin{align} AXA&=A\\ XAX&=X\\ (AX)^*&=AX\\ (XA)^*&=XA\end{align}$$

In terms of the singular value decomposition (SVD), if $A= VDU$, where $V$, $U$ are othogonal matrices and $D$ has zeroes entries except on the main diagonal (which contains the nonnegative square roots of the matrix $A^*A$), then is is easy to see that $U^*D^+V^*$ is the Moore-Penrose pseudo inverse $A^+$ of $A$. It can also be shown that $$A^+=(A^* A)^+A^*$$

There are different methods to show that $A^+$ exists and is unique! This is the Moore-Penrose inverse of $A$. The most efficient method to estimate $A^+$ is through eigenvalue methods, for example the singular value decomposition. Many computer packages (R, Octave, etc) estimate $A^+$ this way (usually coded as pinv(A))

The most important property of $A^+$ is that following:

Theorem: the minimal norm solution to the problem $$\begin{align} \operatorname{min}_{\mathbf{x}\in\mathbb{R}^n}\|A\mathbf{x}-\mathbf{b}\|_2\tag{3}\label{three} \end{align} $$ is $A^+\mathbf{b}$. That is,
$$A^+\mathbf{b}=\operatorname{arg.min}\{\|x\|_2: x=\operatorname{arg.min}\|Ax-b\|_2\} $$

Here is a short proof: Suppose $r=\operatorname{rank}(A)$, and let $A=VDU$ be a singular decomposition of $A$, so that $D\in M(\mathbb{C},m, n)$ has diagonal $\sigma_\geq\sigma_2\geq\ldots\geq\sigma_r>0$. $V$ and $U$ are orthogonal matrices that preserve the distances in $\mathbb{C}^n$ and $\mathbb{C}^m$ respectively (or $\mathbb{R}^n$ and $\mathbb{R}^m$ if you are dealing with real spaces only), that is $$\| Vy\|_2=\|y\|_2,\qquad \|Ux\|_2=\|x\|_2$$ for all $y\in\mathbb{C}^n$ and $x\in\mathbb{C}^m$. Thus $$\begin{align} \rho=\inf_x\|Ax-b\|_2&=\inf_x\|VDU-b\|_2=\inf_x\|V^*(VDUx-b)\|_2\\ &=\inf_x\|DUx-V^*b\|_2=\inf_{y}\|Dy-c\|_2 \end{align} $$ Notice that $$\|Dy-c\|^2_2=\sum^r_{k=1}|\sigma_ky_k-c_k|^2+\sum^m_{k=r+1}|c_k|^2$$ This is minimize when $y_k=c_k/\sigma_k$ for $1\leq k\leq r$ and by letting the remaining $y_{r+1},\ldots,y_n$ to be arbitrary. Therefore $$\rho^2=\sum^m_{k=r+1}|c_k|^2$$ Amongs all vectors $y$ that minimize \eqref{three}, the one with $y_{r+1}=\ldots y_m=0$ has minimum Euclidean norm. Clearly, this vector is given by $$ y=D^+c$$ where $D^+\in M(\mathbb{R},n, m)$ is the Moore-Penrose pseudo inverse of $D$, and is given by the matrix that has the $n$-vector $[\sigma^{-1}_1,\ldots,\sigma^{-1}_r,0,\ldots,0]$ ins the main diagonal and zeroes everywhere else.

Putting this together, the minimal norm solution to the minimum square problem is $$x=U^*y=U^*D^+c=U^*D^+V^*b=A^+b$$

---

Summary:

• As we can see, if $(A^* A)^g$ is a generalized inverse of $A^*A$, then $$\boldsymbol{\omega}=(A^* A)^gA^* \mathbf{y}$$ is a solution to \eqref{one}. Unless $A^* A$ is invertible, there are infinitely many generalized inverse of $(A^* A)$ and thus, many solution to \eqref{one}. (The kernel $A$ provides the infinitely many solutions)

• The particular choice $(A^* A)^+$ however provides the minimal norm solution $\boldsymbol{\omega}$ to \eqref{one}, namely $$\widehat{\boldsymbol{\omega}}=(A^* A)^+A^*\mathbf{y}=A^+\mathbf{y}$$ This is typically the solution that is implemented in OLS algorithms.

• When $A$ is of full rank, $(A^* A)^+=(A^*A)^{-1}$ and in this case $\hat{\boldsymbol{\omega}}=A^+\mathbf{y}=(A^* A)^{-1}A^*\mathbf{y}$ is the unique solution to \eqref{one}

---

Refereces:

1. Albert, A. Regression and the Moore-Penrose pseudo inverse, AP New York, 1972.
2. Rao, C. R. and Mitra, S. K. Generalized Inverse Matrices and its Applications, John Wiley and Sons, New York, 1971
3. Cheney, W. and Kincaid, D. Numerical Analysis, Brooks/Cole Publishing, Pacific Groove California, 1991.