r - 在 R 中实现幂方法时的无限递归

Question

我正在尝试实现用于对矩阵的特征值进行数值评估的幂方法，代码如下:

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)

这个函数获取第一个特征值没有问题，但是在获取第二个特征值时失败了(见代码的最后一行)。你能帮我找出原因吗？谢谢。

我也用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)

我得到了以下错误:

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

Answer 1

我终于让它工作了:

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)

看来python的深度递归有问题，所以我稍微改了一下代码:

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)

我还提出了上述 R 代码的非递归版本:

# 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
}