Python中的矩阵求幂
我试图在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 ，所以任何帮助将不胜感激。
这太有趣了。 我可以说的一件事是，问题是特定于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引发了异常，具体如下：
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是1D或列向量，但np.matrix对象的对角线是行向量。 这凸显了一个更普遍的观点，即np.matrix通常不如np.ndarray支持np.ndarray ，所以在大多数情况下最好使用后者。
一种可能的解决方案是使用np.ravel()将diag_T成一维np.ndarray ：
diag_T = np.ravel(T.diagonal().copy())
这似乎解决了你遇到的问题，尽管可能还有其他与np.matrix有关的问题，我还没有发现。
我在这里打开了拉取请求。
链接地址: http://www.djcxy.com/p/32173.html