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

$\endgroup$
1
$\begingroup$

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

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

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

$\endgroup$
-1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.