17

我需要为很多不同的 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。
4

3 回答 3

3

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

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分量非零，其余分量为零. 这保证了没有特征向量是其他特征向量的线性组合，因此特征向量矩阵是可逆的。

当您的一些eiis 0， or 1， or 重复时，问题就来了。我无法提出证明，但尝试使用以下代码似乎你应该只担心你的任何两个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]])
于 2013-09-20T17:35:57.423 回答
3

我部分回答了我的问题：

根据源代码，我认为 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 值，因为当值太接近时，分解方法是明智的。我很确定这可能会导致不可预测的错误。

于 2013-09-23T16:04:21.607 回答
0

有问题：

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

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

于 2022-02-13T08:30:37.777 回答