I was sitting in my living room chair at 3:30 AM, staring at the ceiling, mulling the meaning of life. I didn’t reach any conclusions so I figured I’d put together a Python language demo of computing the Moore-Penrose pseudo-inverse using three different techniques. If A is the source matrix, then

1. pinv = inv(At * A) * At
2. pinv = V * Sinv * Ut (from SVD)
3. pinv = Rinv * Qt (from QR)

In the first technique, At is the transpose of A, and inv() is any standard matrix inverse function (there are many, such as LUP inverse or Cholesky inverse). This technique is not a true Moore-Penrose pseudo-inverse, which has many rigorous math conditions. This technique is sometimes called the pseudo-inverse using normal equations.

Because At * A is symmetric, positive definite (if you add a small conditioning constant to the diagonal), it’s inverse can be computed using the relatively simple Cholesky decomposition. NumPy does not have a dedicated Cholesky inverse function, so for my demo I used the np.linalg.inv() function, which I believe uses SVD — but I’m not entirely. sure because the documentation is a bit vague.

In the second technique, U, s, Vh is the result of SVD decomposition, and V is the transpose of Vh, Sinv is vector s converted to a diagonal matrix with elements 1/s, and Ut (or Uh) is the transpose of U.

In the third technique, Q, R is the result of QR decomposition, and Rinv is the standard matrix inverse of R (which has a nice shortcut technique), and Qt is the transpose of Q (which is the inverse of Q).

To validate the three techniques, I started by computing the pseudo-inverse using the built-in np.linalg.pinv() function.

The output of my demo is:

A =
[[1.0000 4.0000 2.0000]
 [6.0000 0.0000 3.0000]
 [7.0000 2.0000 1.0000]
 [5.0000 9.0000 8.0000]]

np.linalg pinv(A) =
[[ 0.0047  0.0370  0.1331 -0.0317]
 [ 0.1306 -0.1946  0.1310  0.0239]
 [-0.1158  0.2113 -0.2393  0.1046]]

scratch inv(At * A) * At =
[[ 0.0047  0.0370  0.1331 -0.0317]
 [ 0.1306 -0.1946  0.1310  0.0239]
 [-0.1158  0.2113 -0.2393  0.1046]]

scratch V * Sinv * Ut =
[[ 0.0047  0.0370  0.1331 -0.0317]
 [ 0.1306 -0.1946  0.1310  0.0239]
 [-0.1158  0.2113 -0.2393  0.1046]]

scratch Rinv * Qt =
[[ 0.0047  0.0370  0.1331 -0.0317]
 [ 0.1306 -0.1946  0.1310  0.0239]
 [-0.1158  0.2113 -0.2393  0.1046]]

End pseudo-inverse demo

Each of the three from-scratch techniques has pros and cons — it would take far too long to explain them. Briefly, the standard inverse is good for small matrices. For large matrices the SVD technique is most stable in theory (but I’m skeptical about that claim). For medium size matrices the QR technique is good.

Note that for the first inv(At * A) * At technique, you should add a small constant to the diagonal elements of At * A so that its inverse will always exist.

Good fun at 3:30 AM in my living room (except it’s 4:30 AM now).

To supplement this blog post, I googled about for inverse color photography and found some attractive images. Left and Center by Julia Kuzmenko. Center by Geoffrey Jones.

I could never work in the Arts. In machine learning, success and failure are usually objective and well defined. But in the Arts, everything subjective.

Demo program.

# pinv_explore.py
# three techniques to compute Moore-Penrose pinv

import numpy as np

np.set_printoptions(precision=4, suppress=True,
  floatmode='fixed')

A = np.array([
 [1.0, 4.0, 2.0],
 [6.0, 0.0, 3.0],
 [7.0, 2.0, 1.0],
 [5.0, 9.0, 8.0]])

print("\nA = "); print(A)

def my_pinv1(A):
  # inv(At * A) * At
  At = np.transpose(A)
  tmp1 = np.matmul(At, A)
  tmp2 = np.linalg.inv(tmp1)
  tmp3 = np.matmul(tmp2, At)
  return tmp3

def my_pinv2(A):
  # V * Sinv * Ut
  U, s, Vh = np.linalg.svd(A,
    full_matrices=False)
  Sinv = np.diag(s)
  for i in range(len(s)):
    Sinv[i,i] = 1.0 / Sinv[i,i]
  Ut = np.transpose(U)
  V = np.transpose(Vh)
  tmp = np.matmul(V, Sinv)
  return np.matmul(tmp, Ut)

def my_pinv3(A):
  # Rinv * Qt
  Q, R = np.linalg.qr(A)
  Qt = np.transpose(Q)
  Rinv = np.linalg.inv(R)
  return np.matmul(Rinv, Qt)

Apinv = np.linalg.pinv(A)
print("\nnp.linalg pinv(A) = "); print(Apinv)

B = my_pinv1(A)
print("\nscratch inv(At * A) * At = "); print(B)

C = my_pinv2(A)
print("\nscratch V * Sinv * Ut = "); print(C)

D = my_pinv3(A)
print("\nscratch Rinv * Qt = "); print(D)

print("\nEnd pseudo-inverse demo ")