我需要为很多不同的 N (1 到 10000 之间) 值计算 Q^N，而 Numpy 有点太慢了。

我在 math.stackexchange.com 上问过如果我可以避免为我的特定需要计算 Q^N，并且有人回答我说使用 P D^N P^-1 方法计算 Q^N 应该非常快。

所以基本上，而不是做:

import numpy as np
from numpy import linalg as LA
...
LA.matrix_power(m, N)

我试过了:

diag, P = LA.eig(m)
DN = np.diag(diag**N)
P1 = LA.inv(P)

P*DN*P1

然后我得到与结果相同的矩阵(在单个示例上尝试)

在更复杂的矩阵上，Q:

% timeit.Timer('Q**10000', setup=setup).repeat(2, 100)
[5.87254786491394, 5.863131046295166]

% timeit.Timer('diag, P = linalg.eig(Q); DN=np.diag(diag**10000);P1=linalg.inv(P); P*DN*P1', setup=setup).repeat(2, 100)
[2.0032401084899902, 2.018735885620117]

关于我最初的问题，第二种方法允许我计算 P、diag 和 P1一次并使用它数千次。使用这种方法，速度提高了 8 倍。

我的问题是:

• 在什么情况下不能使用最后一种方法来计算 Q^N？
• 在我的情况下使用它是否合适(给定的矩阵 Q here )？
• 在 numpy 中是否有一个函数已经完成了它？

编辑:

• 对于另一个矩阵，P 似乎是不可逆的。所以我要添加一个新问题:如何更改矩阵 P 使其变得可逆但生成的矩阵不会改变太多？我的意思是，如果这些值接近真实结果就没问题，我的意思是接近 ~0.0001。

最佳答案

您已经知道您的特征值将是 (0, a, b, c, ..., 1)。让我重命名您的参数，使特征值为 (0, e1, e2, e3, ..., 1)。要找出与特征值 ej 对应的特征向量 (v0, v1, v2, ..., v(n-1))，您必须求解方程组:

v1                    = v0*ej
v1*e1 + v2*(1-e1)     = v1*ej
v2*e2 + v3*(1-e2)     = v2*ej
...
vj*ej + v(j+1)*(1-ej) = vj*ej
...
v(n-1)                = v(n-1)*ej

或多或少很清楚，如果你所有的 ei 都是不同的，并且没有一个等于 0 或 1，那么解决方案总是定义良好，并且在处理 ej 时，生成的特征向量的前 j 分量非零，其余分量为零。这保证了没有特征向量是其他特征向量的线性组合，因此特征向量矩阵是可逆的。

当您的某些 ei 为 0 或 1 或重复时，问题就来了。我还没能拿出它的证明，但是用下面的代码试验似乎你应该只担心你的 ei 中的任何两个是否相等并且不同于 1 :

>>> def make_mat(values):
...     n = len(values) + 2
...     main_diag = np.concatenate(([0], values, [1]))
...     up_diag = 1 - np.concatenate(([0], values))
...     return np.diag(main_diag) + np.diag(up_diag, k=1)
>>> make_mat([4,5,6])
array([[ 0,  1,  0,  0,  0],
       [ 0,  4, -3,  0,  0],
       [ 0,  0,  5, -4,  0],
       [ 0,  0,  0,  6, -5],
       [ 0,  0,  0,  0,  1]])
>>> a, b = np.linalg.eig(make_mat([4,5,6]))
>>> a
array([ 0.,  4.,  5.,  6.,  1.])
>>> b
array([[ 1.        ,  0.24253563, -0.18641093,  0.13608276,  0.4472136 ],
       [ 0.        ,  0.9701425 , -0.93205465,  0.81649658,  0.4472136 ],
       [ 0.        ,  0.        ,  0.31068488, -0.54433105,  0.4472136 ],
       [ 0.        ,  0.        ,  0.        ,  0.13608276,  0.4472136 ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.4472136 ]])

现在对于一些测试用例:

>>> a, b = np.linalg.eig(make_mat([1,0,3])) # having a 0 or 1 is OK
>>> b
array([[ 1.        ,  0.70710678,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.70710678,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  1.        ,  0.31622777,  0.57735027],
       [ 0.        ,  0.        ,  0.        ,  0.9486833 ,  0.57735027],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.57735027]])
>>> a, b = np.linalg.eig(make_mat([1,1,3])) # repeating 1 is OK
>>> b
array([[ 1.        ,  0.70710678,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.70710678,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  1.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  1.        ,  0.70710678],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.70710678]])
>>> a, b = np.linalg.eig(make_mat([0,0,3])) # repeating 0 is not OK
>>> np.round(b, 3)
array([[ 1.   , -1.   ,  1.   ,  0.035,  0.447],
       [ 0.   ,  0.   ,  0.   ,  0.105,  0.447],
       [ 0.   ,  0.   ,  0.   ,  0.314,  0.447],
       [ 0.   ,  0.   ,  0.   ,  0.943,  0.447],
       [ 0.   ,  0.   ,  0.   ,  0.   ,  0.447]])
>>> a, b = np.linalg.eig(make_mat([2,3,3])) # repeating other values are not OK
>>> np.round(b, 3)
array([[ 1.   ,  0.447, -0.229, -0.229,  0.447],
       [ 0.   ,  0.894, -0.688, -0.688,  0.447],
       [ 0.   ,  0.   ,  0.688,  0.688,  0.447],
       [ 0.   ,  0.   ,  0.   ,  0.   ,  0.447],
       [ 0.   ,  0.   ,  0.   ,  0.   ,  0.447]])

关于python - 比 numpy 更快的矩阵功率？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/18920125/