# How to determine the N-smallest eigenvalues of a symmetric matrix using the Power Method?

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?

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 $-$.

• 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. – Robst Sep 17 '12 at 1:39

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$.

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()
# plot_all(movieRatings,v1,ax)

fig2 = plt.figure()

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.