How to Find Moore Penrose Inverse

Question

I have a matrix: $$A= \begin{bmatrix} -1 & 0 & 1 & 2 \\ -1 & 1 & 0 & -1 \\ 0 & -1 & 1 & 3 \\ 0 & 1 & -1 & -3 \\ 1 & -1 & 0 & 1 \\ 1 & 0 & -1 & -2 \\ \end{bmatrix} $$ I want to find a Moore Penrose generalized inverse for it, but I only know how to do this using computer software. My attempt: I think the inverse is supposed to be of the form $$C'(CC')^{-1}(B'B)^{-1}B'$$ where $A$ has dimensions $m\times n$, $B$ has dimensions $m\times k$, $C$ has dimensions $k\times n$, and all three have rank $k$. I just don't understand how to actually find this inverse matrix.

Answer 1

After doing Singular-value decomposition of the matrix, it would be very easy to find pseudo inverse of the matrix.

Steps:

Step 1: Find the eigenvalues of $A^TA$ or $AA^T,$ whichever is of less size.

In this case, $A^TA$ would be of less size ($4\times 4$, rather than $6\times 6$).

Step 2: Extract the singular values and frame the $\Sigma$ matrix which is formed by framing a diagonal matrix with diagonal elements equal to the squareroot of each of these eigenvalues and extending it to the shape of $A$ by introducing additional zeroes.

The eigenvalues of $A^TA$ are $\sigma_1^2=34$, $\sigma_2^2=6$, $\sigma_3^2=0$, $\sigma_4^2=0$. Shape of $\Sigma$=Shape of A=$6\times4$.

$\therefore \Sigma=\left[\begin{matrix}\sigma_1 & 0 & 0 & 0\\0 & \sigma_2 & 0 & 0\\0 & 0 & \sigma_3 & 0\\0 & 0 & 0 & \sigma_4\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\end{matrix}\right]=\left[\begin{matrix}5.8311 & 0 & 0 & 0\\0 & 2.4495 & 0 & 0\\0 & 0 & 0.0000 & 0\\0 & 0 & 0 & 0.0000\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\end{matrix}\right]$.

Step 3: Frame V as the matrix with columns as linearly independent orthonormal vectors spacing the domain of $A^TA$.

Since the orthonormal eigenvectors span the domain of $A^TA$, let's find out the eigenvectors of $A^TA$.

The eigenvector corresponding to $\lambda_1=34$ is $v_1=\left[\begin{matrix}0.0648\\0.2593\\-0.3241\\-0.9075\end{matrix}\right]$.

The eigenvector corresponding to $\lambda_2=6$ is $v_2=\left[\begin{matrix}-0.8018\\0.5345\\0.2673\\0.0\end{matrix}\right]$.

The eigenvector corresponding to $\lambda_3=0$ is $v_3=\left[\begin{matrix}0.5773\\0.5773\\0.5773\\0.0\end{matrix}\right]$.

It can be noticed that an independent eigenvector for $\lambda_4=0$ does not exist. Therefore, the other independent vector that together helps in spanning the entire domain can be found by using Gram–Schmidt process. (It is basically taking a random vector and removing the components of other orthonormal vectors from it.)

Therefore, by taking a random vector $r_4=\left[\begin{matrix}1\\1\\1\\1\end{matrix}\right]$,

$v'_4=r_4-(r_4.v_1)v_1-(r_4.v_2)v_2-(r_4.v_3)v_3$ and $v_4=\frac{v'_4}{\left|v'_4\right|}$.

$\Longrightarrow v_4=\left[\begin{matrix}0.14\\0.5602\\-0.7001\\0.42\end{matrix}\right]$

$\therefore V=\left[\begin{matrix}v_1&v_2&v_3&v_4\end{matrix}\right]=\left[\begin{matrix}0.0648 & -0.8018 & 0.5774 & 0.14\\0.2593 & 0.5345 & 0.5774 & 0.5602\\-0.3241 & 0.2673 & 0.5774 & -0.7001\\-0.9075 & 0.0 & 0.0 & 0.42\end{matrix}\right]$.

Step 4: Frame U as the matrix with columns equal to the orthonormal eigenvectors of $AA^T$, and satisfy the equation A=$U\Sigma V^T$.

Let $U=\left[\begin{matrix}u_1&u_2&u_3&u_4&u_5&u_6\end{matrix}\right]$, where $u_1,u_2,u_3,u_4,u_5,u_6$ are the columns of $U$.

Now, $A=U\Sigma V^T \Longrightarrow AV=U\Sigma$

$\Longrightarrow \left[\begin{matrix}-2.2039 & 1.0691 & 0.0 & 0.0\\1.102 & 1.3363 & 0.0 & 0.0\\-3.3059 & -0.2672 & 0.0 & 0.0\\3.3059 & 0.2672 & 0.0 & 0.0\\-1.102 & -1.3363 & 0.0 & 0.0\\2.2039 & -1.0691 & 0.0 & 0.0\end{matrix}\right]=\left[\begin{matrix}\sigma_1 u_1& \sigma_2 u_2&0&0\end{matrix}\right]$.

$\Longrightarrow u_1=\frac{1}{\sqrt{34}}\left[\begin{matrix}-2.2039\\1.102\\-3.3059\\3.3059\\-1.102\\2.2039\end{matrix}\right]=\left[\begin{matrix}-0.378\\0.189\\-0.567\\0.567\\-0.189\\0.378\end{matrix}\right]$ and $u_2=\frac{1}{\sqrt{6}}\left[\begin{matrix}1.0691\\1.3363\\-0.2672\\0.2672\\-1.3363\\-1.0691\end{matrix}\right]=\left[\begin{matrix}0.4365\\0.5455\\-0.1091\\0.1091\\-0.5455\\-0.4365\end{matrix}\right]$.

To find $u_3,u_4,u_5$ and $u_6$, we again use Gram–Schmidt process.

By taking random vector as $rv=\left[\begin{matrix}1\\2\\3\\4\\5\\6\end{matrix}\right]$,

$u'_3=rv-(rv.u_1)u_1-(rv.u_2)u_2$ and $u_3=\frac{u'_3}{\left|u'_3\right|} \Longrightarrow u_3=\left[\begin{matrix}0.3884\\0.4272\\0.4272\\0.3884\\0.3884\\0.4272\end{matrix}\right]$.

Similarly by taking random vectors as $rv_4=\left[\begin{matrix}1\\0\\0\\0\\0\\0\end{matrix}\right]$,$rv_5=\left[\begin{matrix}1\\1\\0\\0\\0\\0\end{matrix}\right]$ and $rv_6=\left[\begin{matrix}1\\1\\1\\0\\0\\0\end{matrix}\right]$

$u'_4=rv_4-(rv_4.u_1)u_1-(rv_4.u_2)u_2-(rv_4.u_3)u_3$ and $u_4=\frac{u'_4}{\left|u'_4\right|} \Longrightarrow u_4=\left[\begin{matrix}0.7182\\-0.4631\\-0.4631\\0.0221\\0.0221\\0.2331\end{matrix}\right]$,

$u'_5=rv_5-(rv_5.u_1)u_1-(rv_5.u_2)u_2-(rv_5.u_3)u_3-(rv_5.u_4)u_4$ and $u_5=\frac{u'_5}{\left|u'_5\right|} \Longrightarrow u_5=\left[\begin{matrix}0.0\\0.5193\\-0.4435\\-0.6207\\0.342\\0.1774\end{matrix}\right]$

$u'_6=rv_6-(rv_6.u_1)u_1-(rv_6.u_2)u_2-(rv_6.u_3)u_3-(rv_6.u_4)u_4-(rv_6.u_5)u_5$ and $u_6=\frac{u'_6}{\left|u'_6\right|} \Longrightarrow u_6=\left[\begin{matrix}0.0\\0.0001\\0.2698\\-0.361\\-0.6309\\0.6316\end{matrix}\right]$

$\therefore U=\left[\begin{matrix}u_1&u_2&u_3&u_4&u_5&u_6\end{matrix}\right]$

$\Longrightarrow U=\left[\begin{matrix}-0.378 & 0.4365 & 0.3884 & 0.7182 & 0.0 & 0.0\\0.189 & 0.5455 & 0.4272 & -0.4631 & 0.5193 & 0.0001\\-0.567 & -0.1091 & 0.4272 & -0.4631 & -0.4435 & 0.2698\\0.567 & 0.1091 & 0.3884 & 0.0221 & -0.6207 & -0.361\\-0.189 & -0.5455 & 0.3884 & 0.0221 & 0.342 & -0.6309\\0.378 & -0.4365 & 0.4272 & 0.2331 & 0.1774 & 0.6316\end{matrix}\right]$

(In this case, alternative matrices for $U$ are possible. And any one amongst those can be considered to be $U$.)

Step 5: Frame the matrix $\Sigma^+$ as the matrix $\Sigma$ with its non-zero singular elements replaced by their corresponding reciprocals and then transposed.

$\therefore \Sigma^+=\left[\begin{matrix}\frac{1}{5.8311} & 0 & 0 & 0\\0 & \frac{1}{2.4495} & 0 & 0\\0 & 0 & 0.00 & 0\\0 & 0 & 0 & 0.00\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\end{matrix}\right]^T=\left[\begin{matrix}0.1715 & 0 & 0 & 0 & 0 & 0\\0 & 0.4082 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\end{matrix}\right]$

Step 6: Finally calculate the Moore-Penrose inverse by the relation $A^+=V \Sigma^+ U^T$.

$\Longrightarrow A^+=\left[\begin{matrix}-0.1471 & -0.1765 & 0.0294 & -0.0294 & 0.1765 & 0.1471\\0.0784 & 0.1274 & -0.049 & 0.049 & -0.1274 & -0.0784\\0.0686 & 0.049 & 0.0196 & -0.0196 & -0.049 & -0.0686\\0.0588 & -0.0294 & 0.0882 & -0.0882 & 0.0294 & -0.0588\end{matrix}\right]$

Answer 2

1. Compute SVD $A=U\Sigma V^T$
2. $A^+=V\Sigma^{+}U^{T}$ where $\Sigma^{+}$ can be computed by taking the inverse of the singular values, the rest of $\Sigma$ should be zero and then transpose.

Answer 3

Here is an alternative way to compute the Moore-Penrose inverse $A'$ of $A$ following Penrose's argument.

---

Let's first state Penrose's generalized inverse theorem for matrices with real entries in a way that is indicative of his proof:

Theorem: Let $A\in\mathcal{M}(m\times n;\mathbb{R})$. Then:

1. There is a unique solution $X\in\mathcal{M}(n\times m;\mathbb{R})$ to the following system of four matrix equations:

\begin{align*} AXA &= A\\ XAX &= X\\ (AX)^T &= AX\\ (XA)^T &= XA. \end{align*}

The unique solution $A'=X$ is by definition the Moore-Penrose inverse of $A$.

2. Either $A=0\in\mathcal{M}(m\times n;\mathbb{R})$, in which case $A'=0\in\mathcal{M}(n\times m;\mathbb{R})$, xor we have that for any $p\in\mathbb{Z}_{\geq0}$, for any $r\in\{1,2,...,p-1\}$ and for any $c_r, c_{r+1},...,c_p\in\mathbb{R}$, if $c_r\neq 0 $ and

$$c_r (A^TA)^r+c_{r+1} (A^TA)^{r+1} + \cdots + c_p (A^TA)^p = 0,$$

then

$$A'=-\dfrac{1}{c_r}(c_{r+1} I + c_{r+2} A^TA + \cdots + c_k (A^TA)^{k-r-1}) A^T.$$

---

Here are some comments regarding item 2 in the above theorem:

• That $0'=0$ follows from uniqueness.

• The fundamental reason behind the expression for $A'$ is that in his proof Penrose considers the iterates of $A^TA$; since the space of $n\times n$ matrices is finite dimensional eventually one obtains a linear relation between iterates.

• By Cayley-Hamilton it is guaranteed that there are such $p,r$ and $c_i$'s for $p\leq n$ (compared to the general bound $p\leq n^2+1$).

• Once there is a linear combination with coefficients $c_1,...,c_p$, the index $r$ is chosen to be the smallest index such that $c_r\neq 0 $.

• If $A\neq0$, $A^TA$ would be a symmetric matrix with at least one nonzero entry, hence it can not be nilpotent (see e.g. Is the zero matrix the only symmetric, nilpotent matrix with real values?); whence $r<p$ in this case.

---

Finally here is how one can use the above theorem to compute the Moore-Penrose inverse $A'$ for the given $A$. Put $S=A^TA$. A calculation gives

$$ S = \begin{pmatrix} 4 & -2 & -2 & -2 \\ -2 & 4 & -2 & -8 \\ -2 & -2 & 4 & 10 \\ -2 & -8 & 10 & 28 \\ \end{pmatrix}; $$

a further calculation gives that $S$ has rank $2$ (or note that the first two rows are linearly independent, the sum of the first three rows is zero and the sum of twice the first row thrice the second row and the fourth row is also zero). Consequently $S$ has nullity $2$ and thus $0$ is an eigenvalue of $S$ of algebraic multiplicity $2$.

Say the characteristic polynomial of $S$ is

$$\chi_S(\lambda) = \lambda^4 - c_3 \lambda^3 + c_2\lambda^2 - c_1 \lambda + c_0.$$

Here $c_3$ depends linearly on eigenvalues of $S$, $c_2$ depends bilinearly, and both $c_1$ and $c_0$ are linear combinations of products of three or four eigenvalues of $S$. Since $S$ has exactly two zero eigenvalues, $c_1=c_0=0$. Further it is known that

$$ c_3 = \operatorname{tr}(S),\quad\quad c_2 = \dfrac{\operatorname{tr}(S)^2-\operatorname{tr}(S^2)}{2}, $$

so that by a calculation we have

$$\chi_S(\lambda) = \lambda^4 - 40 \lambda^3 + 204\lambda^2.$$

Thus by Cayley-Hamilton we have a linear relation as needed in the second item in the above theorem:

$$204 S^2 -40 S^3 + S^4 = 0,$$

whence we have the following expression for the Moore-Penrose inverse:

$$A' = -\dfrac{1}{204}(-40 I + S) A^T. $$

A final calculation thus gives

$$ A' = -\dfrac{1}{204}\begin{pmatrix} 30 & 36 & -6 & 6 & -36 & -30 \\ -16 & -26 & 10 & -10 & 26 & 16 \\ -14 & -10 & -4 & 4 & 10 & 14 \\ -12 & 6 & -18 & 18 & -6 & 12 \\ \end{pmatrix}. $$

(This matches with the already made calculations; see e.g. WolframAlpha.)