比numpy更快的矩阵功率？

Question

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

我已经在math.stackexchange.com上询问是否可以避免为我的特定需求计算Q ^ N而且有人回答我计算Q ^ N应该使用该P D^N P^-1方法非常快.

所以基本上,而不是做:

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 and P1只计算一次并使用它数千次.使用此方法的速度提高了8倍.

我的问题是:

• 在这种情况下,不可能使用这最后一种方法来计算Q ^ N？
• 在我的情况下使用它是否可以(这里给出的矩阵Q )？
• numpy中是否有一个已经完成的功能？

编辑:

• 看起来对于另一个矩阵,P是不可逆的.所以我要添加一个新问题:如何更改矩阵P使其变为可逆,但结果矩阵不会改变太多？我的意思是,如果值接近真实结果就可以了,接近我的意思是~0.0001.

Answer 1

我部分回答我的问题：

根据源代码，我认为 Numpy 使用的是平方求幂：

# binary decomposition to reduce the number of Matrix
# multiplications for n > 3.
beta = binary_repr(n)
Z, q, t = M, 0, len(beta)
while beta[t-q-1] == '0':
    Z = N.dot(Z, Z)
    q += 1
result = Z
for k in range(q+1, t):
    Z = N.dot(Z, Z)
    if beta[t-k-1] == '1':
        result = N.dot(result, Z)
return result

n在我的情况下，当较大时，这比计算特征值和特征向量以及计算 M^N 等于 PD^NP^-1慢。

现在，关于我的问题：

在哪种情况下不可能使用最后一个方法来计算 Q^N？

当某些特征值相等时，将无法反转 P。有人建议在问题跟踪器上的 Numpy 中进行反转。答案是：“你的方法只对无缺陷的稠密矩阵有效。”

在我的情况下使用它可以吗（此处给出的矩阵 Q）？

并非总是如此，我可能有几个相等的特征值。

numpy 中是否有一个函数已经做到了这一点？

我认为它在SciPy中：https://github.com/scipy/scipy/blob/v0.12.0/scipy/linalg/matfuncs.py#L57

所以我们也可以这样做：

LA.expm(n*LA.logm(m))

计算 m^n。

如何更改矩阵 P 使其变得可逆，但生成的矩阵不会发生太大变化？我的意思是，如果这些值接近实际结果就可以了，我所说的接近是指~0.0001。

我不能简单地添加 epsilon 值，因为当值太接近时，分解方法是明智的。我很确定这可能会导致不可预测的错误。