'''
Uses numpy to numerically compute the SVD and compact SVD of a given matrix.
Note that this is a good approximation, but not exact. Small numerical errors
and rounding errors may change the rank of the matrix.
'''

import numpy as np
np.set_printoptions(precision=3) #Number of decimals shown

A = (1/6)*np.array([[4, 5, 2],
                    [0, 3, 6],
                    [4, 5, 2],
                    [0, 3, 6]
                    ])

# Compute the full SVD of A
U, s, VT = np.linalg.svd(A, full_matrices=True)
Sigma=np.zeros(A.shape) #singular values is given in a vector s, put in matrix
Sigma[:min(A.shape), :min(A.shape)] = np.diag(s)

#To form the compact SVD, only save the first r columns of U and V
rank = np.linalg.matrix_rank(A)
U_compact = U[:, :rank]
Sigma_compact=Sigma[:rank,:rank]
VT_compact = VT[:rank, :]

#Display full SVD
print(" U:")
print(U)
print("\n Sigma:")
print(Sigma)
print("\n V:")
print(VT.T)

print("----------")
# Display the compact SVD results

print("\n Compact U:")
print(U_compact)
print("\n Compact Sigma:")
print(Sigma_compact)
print("\n Compact V:")
print(VT_compact.T)

#Uncomment to verify result, these should both be A:
#print(U @ Sigma @ VT)
#print(U_compact @ Sigma_compact @ VT_compact)