Numerical stability of a matrix inversion problem

121 Views Asked by At

Suppose, we have a matrix X with dimension m and n. Now lets say, we are interested to find the inverse of the covariance matrix: $X^T X$. The straightforward way to do it is to compute the matrix $X^T X$ and perform the inversion operation $(X^T X)^{-1}$. But we also know from SVD decomposition

$$X= USV^T \\ X^TX = VS^TUU^TSV^T=VS^TSV^T=VS^2V^T \\ \text{since } (X^T X)^{-1}(X^T X) = I \\ \text{and }(V(S^2)^{-1}V^T)(VS^2V^T )= V(S^2)^{-1}S^2V^T = VV^T = I$$

we can alternatively calculate $V(S^2)^{-1}V^T$ to find the inverse of $X^T X$. But I have found experimentally that for m>>n, both yield the same answer but for m << n (e.g., m=100, n=1000), this is not the case. I have attached a python code snippet for this also.

Can someone explain why and if there is any remedy for it? Thanks in advance!

import numpy as np


def inverse_delta(x):
  u, s, vh = np.linalg.svd(x,  full_matrices=False)
  v = vh.T
  inverse = np.linalg.inv(x.T@x)
  inverse_tilde = ([email protected](1/(s**2)))@v.T
  print(np.linalg.norm(inverse-inverse_tilde))


inverse_delta(x = np.random.randn(1000, 100))
3.448652345497296e-17
inverse_delta(x = np.random.randn(100, 1000))
8797061188866344.0