# 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 [ ]: