Python中的矩阵求幂

Posted 2023-03-12

【中文标题】Python中的矩阵求幂

【英文标题】：Matrix exponentiation in Python

【发布时间】：2016-05-30 21:37:07

【问题描述】：

我正在尝试在 Python 中对一个复杂的矩阵求幂，但遇到了一些麻烦。我正在使用scipy.linalg.expm 函数，当我尝试以下代码时出现了一个相当奇怪的错误消息：

import numpy as np
from scipy import linalg

hamiltonian = np.mat('[1,0,0,0;0,-1,0,0;0,0,-1,0;0,0,0,1]')

# This works
t_list = np.linspace(0,1,10)
unitary = [linalg.expm(-(1j)*t*hamiltonian) for t in t_list]

# This doesn't
t_list = np.linspace(0,10,100)
unitary = [linalg.expm(-(1j)*t*hamiltonian) for t in t_list]

第二个实验运行时的错误是：

This works!
Traceback (most recent call last):
  File "matrix_exp.py", line 11, in <module>
    unitary_t = [linalg.expm(-1*t*(1j)*hamiltonian) for t in t_list]
  File "/usr/lib/python2.7/dist-packages/scipy/linalg/matfuncs.py",     line 105, in expm
    return scipy.sparse.linalg.expm(A)
  File "/usr/lib/python2.7/dist- packages/scipy/sparse/linalg/matfuncs.py", line 344, in expm
    X = _fragment_2_1(X, A, s)
  File "/usr/lib/python2.7/dist-  packages/scipy/sparse/linalg/matfuncs.py", line 462, in _fragment_2_1
    X[k, k] = exp_diag[k]
TypeError: only length-1 arrays can be converted to Python scalars

这看起来很奇怪，因为我改变的只是我使用的t 的范围。是因为哈密顿量是对角线的吗？一般来说，哈密顿量不会，但我也希望它适用于对角线。我真的不知道expm 的机制，所以任何帮助将不胜感激。

【问题讨论】：

您可以尝试将计算移至 for 循环而不是列表理解。然后你至少可以弄清楚它失败的 t 值。

程序失败的第一个数字是t=2.12121212121。这似乎完全是任意的......该程序不适用于t=2.ax，其中a > 0。而且它根本不适用于t=3.x ...

【参考方案1】：

这很有趣。我可以说的一件事是问题是特定于np.matrix 子类的。例如，以下工作正常：

h = np.array(hamiltonian)
unitary = [linalg.expm(-(1j)*t*h) for t in t_list]

---

深入挖掘回溯，在_fragment_2_1 中的scipy.sparse.linalg.matfuncs.py 中引发了异常，特别是these lines：

n = X.shape[0]
diag_T = T.diagonal().copy()

# Replace diag(X) by exp(2^-s diag(T)).
scale = 2 ** -s
exp_diag = np.exp(scale * diag_T)
for k in range(n):
    X[k, k] = exp_diag[k]

错误信息

    X[k, k] = exp_diag[k]
TypeError: only length-1 arrays can be converted to Python scalars

向我建议 exp_diag[k] 应该是一个标量，而是返回一个向量（并且您不能将向量分配给 X[k, k]，它是一个标量）。

设置断点并检查这些变量的形状证实了这一点：

ipdb> l
    751     # Replace diag(X) by exp(2^-s diag(T)).
    752     scale = 2 ** -s
    753     exp_diag = np.exp(scale * diag_T)
    754     for k in range(n):
    755         import ipdb; ipdb.set_trace()  # breakpoint e86ebbd4 //
--> 756         X[k, k] = exp_diag[k]
    757 
    758     for i in range(s-1, -1, -1):
    759         X = X.dot(X)
    760 
    761         # Replace diag(X) by exp(2^-i diag(T)).

ipdb> exp_diag.shape
(1, 4)
ipdb> exp_diag[k].shape
(1, 4)
ipdb> X[k, k].shape
()

根本问题是exp_diag 被假定为一维向量或列向量，但np.matrix 对象的对角线是行向量。这突出了一个更普遍的观点，即np.matrix 通常不如np.ndarray 得到很好的支持，因此在大多数情况下最好使用后者。

一种可能的解决方案是使用np.ravel() 将diag_T 扁平化为一维np.ndarray：

diag_T = np.ravel(T.diagonal().copy())

这似乎解决了您遇到的问题，尽管我还没有发现与np.matrix 相关的其他问题。

---

我已经打开了一个拉取请求here。

【讨论】：

非常感谢您！我想我可能只使用ndarray 而不是mat 并在需要时将其转换回矩阵。我想我对这两个类的内部运作了解不够。不相关，但是如何在 Python 中设置断点？您使用的是特定的 IDE 还是可以在 emacs 等中使用？

您不需要 IDE 来设置断点。我使用的是ipdb，但您也可以使用标准的Python 调试器pdb。我使用 SublimeText3 作为编辑器，它有一个扩展来方便地设置断点，但肯定有一个 Emacs 等价物......