import numpy as np
import pandas as pd
import math
import scipy as sp
import matplotlib.pyplot as plt
import seaborn; seaborn.set()
from sklearn.linear_model import LinearRegression

## random number generator
from numpy.random import default_rng
rng = np.random.default_rng(seed=42)

## X is n x p, iid N(0,1)
n = 20; p=3
X = np.random.normal(size=(n,p))

## get svd and check basic dimensions
tsvd = np.linalg.svd(X)
print("the length of tsvd is:",len(tsvd))
print("dim of first: ",tsvd[0].shape)
print("dim of second: ",tsvd[1].shape)
print("dim of third: ",tsvd[2].shape)

the length of tsvd is: 3
dim of first:  (20, 20)
dim of second:  (3,)
dim of third:  (3, 3)

## get U and check that U'U = I
U = tsvd[0]
U.T @ U
print("U'U, should be identity:\n",(U.T @ U)[:3,:3])
print("diagonals: ",np.diag(U.T @ U))

U'U, should be identity:
 [[ 1.00000000e+00 -3.51581677e-16  4.61260460e-17]
 [-3.51581677e-16  1.00000000e+00 -2.68452368e-17]
 [ 4.61260460e-17 -2.68452368e-17  1.00000000e+00]]
diagonals:  [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]

## get Sigma, this is the r diagonals
Sigma = tsvd[1]
print("Sigma:",Sigma)
St = np.diag(Sigma) # rxr Sigma upper left corner.
# r is rank of X.
r = len(Sigma)
print("rank of X: ",r)

Sigma: [4.82887546 4.50029029 2.38775877]
rank of X:  3

## get V, this is actually V' 
V = tsvd[2]
print("VV', should be identity:\n",V @ V.T)

VV', should be identity:
 [[1.00000000e+00 1.40814493e-16 5.08298617e-17]
 [1.40814493e-16 1.00000000e+00 4.36456746e-17]
 [5.08298617e-17 4.36456746e-17 1.00000000e+00]]

## reconstruct X from the reduced form version of the SVD
Ut = U[:,0:r] # first r columns of U
print(f'number of columns in U: " {U.shape[1]}\n')
print(f'number of columns in Ut: " {Ut.shape[1]}\n')

number of columns in U: " 20

number of columns in Ut: " 3

## this should be X
check = Ut @ St @ V

## check that check = X
cdif = check-X
print("\nshould be all 0 !!!!!\n")
print(np.abs(cdif).sum())

should be all 0 !!!!!

2.3425705819590803e-14

## check linear regression
## simulate y
beta = np.array([1,2,3]).reshape((p,1))
sigma = .5
y = X @ beta + sigma * np.random.normal(size=(n,1))
plt.scatter(X[:,2],y)

<matplotlib.collections.PathCollection at 0x760eabf1abf0>

## fit reg in sklearn
lmod = LinearRegression(fit_intercept=False)
lmod.fit(X,y)
yhat = lmod.predict(X)
plt.scatter(y,yhat)
plt.plot(y,y,c='red')

[<matplotlib.lines.Line2D at 0x760eabf9fbb0>]

## Xp = Ut * Sigma^(-1) * V, this is also the moore-penrose inverse
Xp = V.T @ np.linalg.inv(St) @ Ut.T
print('\nXp * y:\n',Xp @ y)
print('\ncheck:\n', lmod.coef_)

Xp * y:
 [[1.18160334]
 [2.43027246]
 [3.0619344 ]]

check:
 [[1.18160334 2.43027246 3.0619344 ]]

## yhat from Ut, col's of Ut are ON basis for column space of X
yhat1 = Ut @ Ut.T @ y
check = yhat - yhat1
print("check should be 0: ",np.sum(np.abs(yhat-yhat1)))

## yhat should also be XXp+ y.
yhat2 = X @ Xp @ y
print("check should be 0: ",np.sum(np.abs(yhat-yhat2)))

check should be 0:  3.697042672001771e-14
check should be 0:  1.2434497875801753e-14

##################################################
## let's see what happens in the non full-rank case
## we will append X with the sum of the columns of X
Xsum = X.sum(axis=1)
## check Xp * Xsum should be all ones
print(Xp @ Xsum)

[1. 1. 1.]

## try one of these two
X1 = np.hstack([X,Xsum.reshape((n,1))]) # add on mean of columns
#X1 = np.hstack([X,X[:,p-1].reshape(n,1)]) # repeat last column

## get svd and check basic dimensions
tsvd = np.linalg.svd(X1)
print("the length of tsvd is:",len(tsvd))
print("dim of first: ",tsvd[0].shape)
print("dim of second: ",tsvd[1].shape)
print("dim of third: ",tsvd[2].shape)

the length of tsvd is: 3
dim of first:  (20, 20)
dim of second:  (4,)
dim of third:  (4, 4)

## check out the singular values, notice that the last one is basically 0
print(tsvd[1])

## set rank (number of nonzero diagonal elements of Sigma)
r=3

## reconstruct X from the reduced form version of the SVD
Ut1 = tsvd[0][:,0:r] # first r columns of U

yhat3 = Ut1 @ Ut1.T @ y
print("check should be 0: ",np.sum(np.abs(yhat-yhat3)))

[8.38824424e+00 4.53109251e+00 2.73044217e+00 9.73053443e-16]
check should be 0:  4.884981308350689e-14

## Make generalized inverse
Vt1 = (tsvd[2].T)[:,0:r]
Xp1 = Vt1 @ np.diag(1/tsvd[1][0:r]) @ Ut1.T
bo1 = Xp1 @ y

yhat4 = X1  @ bo1
print("check should be 0: ",np.sum(np.abs(yhat-yhat4)))

check should be 0:  1.587618925213974e-14

check1 = Ut1 @ np.diag(tsvd[1][0:r]) @ Vt1.T

## check that check1 = X
cdif = check1-X1
print("should be all 0 !!!!!")
print(np.abs(cdif).sum())

should be all 0 !!!!!
7.724376693829527e-14

Simple SVD¶