3

我正在尝试实现用于数值评估矩阵特征值的幂方法，这是代码：

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'
4

1 回答 1

2

我终于让它工作了：

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
}
于 2013-07-28T14:05:39.323 回答