比numpy更快的矩阵幂运算？

Question

我需要计算许多不同的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如此处所示）使用它是否合适？
• 在numpy中是否有已经实现了这一操作的函数？

编辑：

• 对于另一个矩阵，P不可逆。因此，我添加了一个新问题：如何更改矩阵P以使其可逆，但生成的矩阵不会过于改变？我的意思是，如果值接近实际结果，那么就可以；接近的定义是 ~0.0001。

- Maxime Chéramy

---

关于我的前两个问题，我猜测 Q 不应该有缺陷。但是我不知道我的矩阵是否有缺陷（因为我的数学背景太远了）。 - Maxime Chéramy

你可以使用平方取幂法通过 diag**10000 进一步加速。请参见我的另一个回答，其中我在numpy中实现了它。 - Claudiu

@Claudiu 哇，我天真地以为 diag**10000 会使用平方方法！然而，在我的情况下，我甚至可以使用传递浮点数的可能性。这也是我无法使用 LA.matrix_power 的东西。 - Maxime Chéramy

只要矩阵可对角化，你就可以做到这一点。如果矩阵Q是实数对称矩阵，你总是能够执行此操作。祝你好运，Maxime。 - Wok

我同意wok的观点，正如我之前所说，我不知道我的矩阵是否有缺陷（维基百科上说：“一个非对角化的方阵被称为有缺陷的。”）。矩阵Q不对称（它是三角形的）。也许我应该对矩阵的行列式做些什么，不知道，我想我需要深入研究数学定理。谢谢；)。 - Maxime Chéramy

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计算为等于P D^N P^-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值，因为当值太接近时，分解方法容易出现问题。我相信这可能会导致不可预测的错误。

- Maxime Chéramy

---

矩阵指数选项很有趣（即使用la.expm），但在我的机器上似乎比matrix_power还要慢。尽可能对角化可能是最好的选择。 - IanH

Answer 2

您已经知道了特征值是(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]])

- Jaime

---

我正在阅读您的答案，谢谢。但是有一个实时评论：我的矩阵每行的总和为1（因为它是马尔可夫链的转移矩阵）。 - Maxime Chéramy

是的，这已经被考虑在内了。a、b、c...不应该相互等同，除非它们都等于1。 - Jaime

啊啊确实是这样的 :). 但是，a、b、c等实际上是概率。这就是我应该说的，抱歉。 - Maxime Chéramy

Answer 3

问题是：

在哪种情况下无法使用此最后一种方法计算 Q^N？

关键思想是判断 Q 是否可对角化。等价地，我们应该判断 Q 是否有 n（其行/列数）个线性无关的特征向量。请注意，不同的特征值只是对角化的充分条件，但不是必要条件。

- user17412768