Normal Equations vs Pseudoinverse

In [3]:
import numpy as np
import numpy.linalg as la

Here's a simple overdetermined linear system, which we'll solve using both the normal equations and the pseudoinverse:

In [4]:
A = np.random.randn(5, 3)
b = np.random.randn(5)

Normal Equations

Solve $Ax\cong b$ using the normal equations:

In [5]:
x1 = la.solve(A.T@A, A.T@b)
x1
Out[5]:
array([ 0.15751782,  0.37883556,  1.02377709])

Pseudoinverse

Solve $Ax\cong b$ using the pseudoinverse:

In [6]:
U, sigma, VT = la.svd(A)
print(U)
print(sigma)
print(VT)
[[ 0.40189265  0.43776205 -0.32197721 -0.68693466 -0.26701711]
 [-0.60148796 -0.48057566 -0.36776893 -0.32823601 -0.40529791]
 [-0.50534567  0.25317047  0.34332919 -0.50930065  0.55069808]
 [-0.46931436  0.7148832  -0.11893232  0.37004742 -0.34293739]
 [-0.03262382  0.04751987 -0.79313105  0.15511173  0.5861408 ]]
[ 2.48597339  2.0073724   1.09629545]
[[-0.16519593 -0.85534195  0.49101981]
 [ 0.90601342 -0.3283188  -0.26710758]
 [ 0.38967935  0.40074545  0.82918821]]
In [7]:
Sigma_inv = np.zeros_like(A.T)
Sigma_inv[:3,:3] = np.diag(1/sigma)
Sigma_inv
Out[7]:
array([[ 0.40225692,  0.        ,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.49816367,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.91216287,  0.        ,  0.        ]])
In [10]:
pinv = VT.T @ Sigma_inv @ U.T
x2 = pinv @ b
x2
Out[10]:
array([ 0.15751782,  0.37883556,  1.02377709])
In [9]:
la.norm(x1-x2)
Out[9]:
4.7428748402675471e-16
In [ ]: