在python中计算矩阵的指数

Question

我想计算一个 200x200 矩阵 (expm(B)) 的指数并得到以下问题。非常感谢您的帮助。

exp_matrix2 = expm(B)

文件“.../python2.7/site-packages/scipy/linalg/matfuncs.py”，行 261，在 expm 返回 scipy.sparse.linalg.expm(A)

文件“.../python2.7/site-packages/scipy/sparse/linalg/matfuncs.py”， 第 582 行，在 expm 中 return _expm(A, use_exact_onenorm='auto')

文件“.../python2.7/site-packages/scipy/sparse/linalg/matfuncs.py”， 第 618 行，在 _expm 中 eta_1 = max(h.d4_loose, h.d6_loose)

文件“.../python2.7/site-packages/scipy/sparse/linalg/matfuncs.py”， 第 457 行，在 d4_loose 结构=self.structure)**(1/4.)

文件“.../python2.7/site-packages/scipy/sparse/linalg/matfuncs.py”， 第 301 行，在 _onenormest_matrix_power MatrixPowerOperator(A, p, structure=structure))

文件“.../python2.7/site-packages/scipy/sparse/linalg/_onenormest.py”， 第 95 行，最正常的 est, v, w, nmults, nresamples = _onenormest_core(A, A.H, t, itmax)

文件“/python2.7/site-packages/scipy/sparse/linalg/_onenormest.py”， 第 424 行，在 _onenormest_core Z = np.asarray(AT_linear_operator.matmat(S))

文件“.../python2.7/site-packages/scipy/sparse/linalg/interface.py”， 第 326 行，matmat Y = self._matmat(X)

文件“.../python2.7/site-packages/scipy/sparse/linalg/interface.py”， 第 468 行，在 _matmat 中 return super(_CustomLinearOperator, self)._matmat(X)

文件“.../python2.7/site-packages/scipy/sparse/linalg/interface.py”， 第 174 行，在 _matmat return np.hstack([self.matvec(col.reshape(-1,1)) for col in X.T])

文件 “/home/dk2518/anaconda2/lib/python2.7/site-packages/scipy/sparse/linalg/interface.py”， 第 219 行，在 matvec 中 y = self._matvec(x)

文件“.../python2.7/site-packages/scipy/sparse/linalg/interface.py”， 第 471 行，在 _matvec 中 返回 self.__matvec_impl(x)

文件“.../python2.7/site-packages/scipy/sparse/linalg/interface.py”， 第 266 行，在 rmatvec 中 y = self._rmatvec(x)

文件“.../python2.7/site-packages/scipy/sparse/linalg/matfuncs.py”， 第 203 行，在 _rmatvec x = A_T.dot(x)

ValueError：形状 (207,207) 和 (1,207) 未对齐：207 (dim 1) != 1 (暗淡 0)

Answer 1

正如@TomNash 的链接中所讨论的，一个大的np.matrix 是问题所在。

ndarray 和 sparse 矩阵工作正常：

In [309]: slg.expm(np.ones((200,200)));                                         
In [310]: slg.expm(sparse.csc_matrix(np.ones((200,200))));                      
In [311]: slg.expm(np.matrix(np.ones((200,200))));  
ValueError: shapes (200,200) and (1,200) not aligned: 200 (dim 1) != 1 (dim 0)

不是每个np.matrix 都会出问题：

In [313]: slg.expm(np.matrix(np.eye(200)));

将np.matrix 转回ndarray 有效：

In [315]: slg.expm(np.matrix(np.ones((200,200))).A);

这使用slg.matfuncs._expm(A, use_exact_onenorm='auto')

其中有一个测试，早期的，用于：

if use_exact_onenorm == "auto":
    # Hardcode a matrix order threshold for exact vs. estimated one-norms.
    use_exact_onenorm = A.shape[0] < 200

这部分解释了为什么我们会遇到 (200,200) 矩阵而不是 (199,199) 的问题。

这行得通：

slg.matfuncs._expm(M, use_exact_onenorm=True);

它以False 失败。但是从那里我迷失了它如何设置_ExpmPadeHelper并尝试Pade order 3的细节。

简而言之 - 避免使用np.matrix，尤其是在 (200,200) 或更大的情况下。