Infinite recursion when implementing power method in R

Question

I am trying to implemente the power method for numerically evaluating the eigenvalues of a matrix, here is the code:

A = matrix(c(1, 1, 2, 0), 2, 2, byrow=TRUE)
x0 = c(1, 0)

powerm = function(A, x0, thresh)
{
    m0 = max(x0)
    x1 = A %*% (x0 / m0)
    m1 = max(x1)
    if(abs(m1 - m0) < thresh)
    {
        return(m1)
    }
    else
    {
        powerm(A, x1, thresh)
    }
}

ev1 = powerm(A, x0, 1e-4)
ev2 = powerm(A - diag(2)*ev1, x0, 1e-4)

This function gets the first eigenvalue without problem, but fails when getting the second one (see the last line of the code). Could you please help me find out why? Thanks.

I have also rewritten it in python:

import numpy as np

A = np.matrix([[1, 1], [2, 0]])
x0 = np.matrix([1, 0]).reshape(2, 1)

def powerm(A, x0, thresh):
    m0 = x0.max()
    x1 = A * (x0 / m0)
    m1 = x1.max()
    if abs(m1 - m0) < thresh:
        return m1
    else:
        return powerm(A, x1, thresh)

ev1 = powerm(A, x0, 1e-12)
A1 = A - np.identity(2) * ev1
ev2 = powerm(A1, x0, 1e-12)

And I got the following error:

AttributeError: 'numpy.ndarray' object has no attribute '_collapse'

Answer 1

I finally got it working:

A = matrix(c(1, 1, 2, 0), 2, 2, byrow=TRUE)
x0 = rnorm(2)
thresh = 1e-22

powerm = function(A, x0, thresh)
{
    m0 = x0[which.max(abs(x0))]
    x1 = A %*% (x0 / m0)
    m1 = x1[which.max(abs(x1))]
    cat(m1, '\n')
    if(abs(m1 - m0) < thresh)
    {
        return(m1)
    }
    else
    {
        powerm(A, x1, thresh)
    }
}

ev1 = powerm(A, x0, thresh)
A1 = A - diag(2)*ev1
ev2 = ev1 + powerm(A1, x0, thresh)

It looks like python has a problem with deep recursion, so I changed the code a bit:

import numpy as np

A = np.matrix([[1, 1], [2, 0]])
x0 = np.matrix([1, 0]).reshape(2, 1)
thresh = 1e-33

def powerm(A, x0, thresh):
    m0 = x0.flat[abs(x0).argmax()]
    x1 = A * (x0 / m0)
    m1 = x1.flat[abs(x1).argmax()]
    while abs(m1 - m0) > thresh:
        m0 = m1
        x1 = A * (x1 / m1)
        m1 = x1.flat[abs(x1).argmax()]
    return m1;

ev1 = powerm(A, x0, thresh)
A1 = A - np.identity(2) * ev1
ev2 = ev1 + powerm(A1, x0, thresh)

I also came up with a non-recursive version of the R code above:

# power method without recursion
powerm_nr = function(A, x0, thresh)
{
    m0 = x0[which.max(abs(x0))]
    x1 = A %*% (x0 / m0)
    m1 = x1[which.max(abs(x1))]
    cat(m1, '\n')
    while(abs(m1 - m0) > thresh)
    {
        m0 = m1
        x1 = A %*% (x1 / m1)
        m1 = x1[which.max(abs(x1))]
    }
    m1
}