Frobenius 范数不是算子范数，它是线性算子/矩阵向量空间的范数，这不是一回事。只需将其更改为任何其他预设规范，它应该可以工作。

您计算矩阵幂的方法也不稳定。Numpy 中使用的算法是基本的重复平方，它没有归一化或正则化步骤，并且对于较大的可能会爆炸。计算高矩阵幂的一种更稳定的方法是使用特征分解，但您可以直接从中计算光谱半径，因此看起来没有实际意义。你能更准确地描述你的目标和局限性吗？你这样做是为了学校吗？你的矩阵是不是太大了？等等。$k$

事实证明，这个问题非常简单（与 numpy 的矩阵幂函数本身的限制无关）。我最初认为传播了一些数值浮点错误 - 但我正在测试的示例矩阵仅包含整数。问题在于，在时，矩阵条目的值超过了 numpy 的最大可能 int64 值——它们在那里变成了的数量级。$k=20$$10^{20}$

如果我改变

a = np.array([[9, -1, 2], [-2, 8, 4], [1, 1, 8]])

到

a = np.array([[9.0, -1, 2], [-2, 8, 4], [1, 1, 8]])

在我的代码中，当我们增加的值时，Gelfand 的公式会按预期收敛。$k$

一些情节来证明这个问题。第一个（不正确的）对应于上面的整数类型的 numpy 数组，而第二个（正确的）对应于浮点类型的数组。

这个解决方案并不令人惊讶，但也许在这里发布它会帮助遇到类似问题的其他人。随着的增加，您仍然会遇到数值问题（我开始处获取inf值），因此我认为如果您想做得更好，则需要更复杂的方法。$k$$k=154$

代码：

import matplotlib.pyplot as plt
import numpy as np
from numpy import linalg as la

def kth_term(matrix, k=20):
    """Gelfand's formula"""
    matrix = la.matrix_power(matrix, k)
    f_norm = la.norm(matrix, 'fro')
    term = f_norm**(1.0/k)

    return term


def test_convergence_of_gelfands_formula():
    a = np.array([[9.0, -1, 2], [-2, 8, 4], [1, 1, 8]])

    for k in range(1, 100):
        term = kth_term(a, k)
        print("k={}: {}".format(k, term))
        plt.plot(k, term, 'rx')

    plt.plot(k, term, 'rx', label="kth term")
    plt.title("Convergence of Gelfand's formula for the spectral radius")
    plt.xlabel("k")
    plt.axhline(y=10, color='b', linestyle='-', label='Spectral radius')
    plt.legend()
    plt.show()


def main():
    test_convergence_of_gelfands_formula()


if __name__ == '__main__':
    main()