2
$\begingroup$

I was assigned to make a program that finds the largest, the N-largest, the smallest and the N-smallest eigenvalues of a symmetric matrix, using the Power Method. So far, I've been able to succesfully calculate the largest eigenvalue using the traditional Power Method, the N-largest using the Power Method with Deflation, and the smallest using the Inverse Iteration (the Inverse Iteration as described here in section 3-2: Iterative Methods).

But, right now I have no idea how to determine the N-smallest. I tried using the Inverse Iteration with Deflation to calculate the N-smallest but it is not working. When I calculate the second smallest (and so on...) I don't get the expected results, as if it's not possible to simply apply the Inverse Iteration with deflation. What am I missing? What should be the right way to calculate de N-smallest?

Your help is deeply appreciated.

3 Answers 3

1

I assume by "largest" and "smallest" you mean in absolute value. The $n$'th smallest eigenvalue of $A$ is $\pm \sqrt{t - \mu}$ where $\mu$ is the $n$'th largest eigenvalue of $t I - A^2$ if $t$ is large enough. See how this moves if you replace $A$ by $A-\epsilon I$ for small $\epsilon$ and you can tell whether it's $+$ or $-$.

  • 0
    Yes, I mean in absolute value, thanks for your reply. I was totally unaware of that way of calculating the N-smallest eigenvalues of a matrix, and definitely will keep it in mind, next time I need to do it. Now, since my homework requires me to do it using the Power Method, I wonder if you know how to determine the N-smallest eigenvalues using the Power Method, I feel like it's easy and I'm missing something obvious but I can't see it.2012-09-17
0

Is your matrix positive-definite ? (all eigenvalues positives). If yes, let $M$ be the largest eigenvalue. Find the N-largest eigenvalues of $A-MId$ by power method w/ deflation to find the N smallest eigenvalues of $A$.

0

Try this,

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.