import numpy as np from numpy.linalg import norm
from random import normalvariate from math import sqrt import
matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D
def randomUnitVector(n):
unnormalized = [normalvariate(0, 1) for _ in range(n)]
theNorm = sqrt(sum(x * x for x in unnormalized))
return [x / theNorm for x in unnormalized]
def svd_1d(A, epsilon=1e-10,orth=[]):
''' The one-dimensional SVD, find smallest eigen '''
n, m = A.shape
x = randomUnitVector(min(n,m))
lastV = None
currentV = x
if n > m:
B = np.linalg.inv(np.dot(A.T, A))
else:
B = np.linalg.inv(np.dot(A, A.T))
iterations = 0
while True:
iterations += 1
lastV = currentV
currentV = np.dot(B, lastV)
for v in orth:
currentV -= currentV.dot(v) * v #ensure orthogonality with vectors in orth
currentV /= np.linalg.norm(currentV)
if abs(np.dot(currentV, lastV)) > 1 - epsilon:
# print("converged in {} iterations!".format(iterations))
return currentV
def svd(A, k=None, epsilon=1e-10):
'''
Compute the singular value decomposition of a matrix A
using the power method. A is the input matrix, and k
is the number of singular values you wish to compute.
If k is None, this computes the full-rank decomposition.
'''
rng = np.random.RandomState(42)
A = np.array(A, dtype=float)
n, m = A.shape
svdSoFar = []
if k is None:
k = min(n, m)
for i in range(k):
matrixFor1D = A.copy()
orth = [np.array(v) for singularValue, u, v in svdSoFar[:i]]
if n > m:
v = svd_1d(matrixFor1D, epsilon=epsilon,orth=orth) # next singular vector
u_unnormalized = np.dot(A, v)
sigma = norm(u_unnormalized) # next singular value
u = u_unnormalized / sigma
else:
u = svd_1d(matrixFor1D, epsilon=epsilon,orth=orth) # next singular vector
v_unnormalized = np.dot(A.T, u)
sigma = norm(v_unnormalized) # next singular value
v = v_unnormalized / sigma
svdSoFar.append((sigma, u, v))
singularValues, us, vs = [np.array(x) for x in zip(*svdSoFar)]
return singularValues, us.T, vs
def plot_all(A,v,ax):
ax.scatter(A[:,0],A[:,1],A[:,2])
x=[0,v[0]]
y=[0,v[1]]
z=[0,v[2]]
ax.set_xlim(-1, 1)
ax.set_ylim(-1, 1)
ax.set_zlim(-1, 1)
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
ax.plot(x,y,z)
if name == "main":
def funC(t,rng,r):
return t +r*2*(rng.rand()-0.5)
rng = np.random.RandomState(42)
xaux=[-1+float(i)/float(50) for i in range(100)]
yaux=[funC(i,rng,0.5) for i in xaux]
zaux=[funC(i,rng,0.01) for i in xaux]
movieRatings= np.array([[xaux[i],yaux[i],zaux[i]] for i in range(100)])
# v1 = svd_1d(movieRatings,orth=[np.array([1,0,0])])
S,_,V=(svd(movieRatings,k=2))
# fig = plt.figure()
# ax = fig.add_subplot(111, projection='3d')
# plot_all(movieRatings,v1,ax)
fig2 = plt.figure()
ax2 = fig2.add_subplot(111, projection='3d')
print (S)
print (V)
for v in V:
plot_all(movieRatings,v,ax2)
plt.show()
# print(v1)
#print(svd(movieRatings,k=1))
apologies for the messy code.