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? 

Your help is deeply appreciated.
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 $-$.
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 &gt; 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)) &gt; 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 &gt; 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.