Python中的矩阵求幂 [英] Matrix exponentiation in Python
问题描述
我正在尝试对Python中的复杂矩阵求幂,并且遇到了一些麻烦.我正在使用scipy.linalg.expm
函数,并且在尝试以下代码时出现了一个非常奇怪的错误消息:
I'm trying to exponentiate a complex matrix in Python and am running into some trouble. I'm using the scipy.linalg.expm
function, and am having a rather strange error message when I try the following code:
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]
运行第二个实验时的错误是:
The error when the second experiment is run is:
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
的机制,因此将不胜感激.
This seems really strange since all I changed was the range of t
I was using. Is it because the Hamiltonian is diagonal? In general, the Hamiltonians won't be, but I also want it to work for diagonal ones. I don't really know the mechanics of expm
, so any help would be greatly appreciated.
推荐答案
这很有趣.我可以说的一件事是,该问题特定于np.matrix
子类.例如,以下方法可以正常工作:
That is interesting. One thing I can say is that the problem is specific to the np.matrix
subclass. For example, the following works fine:
h = np.array(hamiltonian)
unitary = [linalg.expm(-(1j)*t*h) for t in t_list]
深入研究追溯,在scipy.sparse.linalg.matfuncs.py
的_fragment_2_1
中引发了异常,特别是
Digging a little deeper into the traceback, the exception is being raised in _fragment_2_1
in scipy.sparse.linalg.matfuncs.py
, specifically 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]
).
suggests to me that exp_diag[k]
ought to be a scalar, but is instead returning a vector (and you can't assign a vector to X[k, k]
, which is a scalar).
设置断点并检查这些变量的形状可以确认这一点:
Setting a breakpoint and examining the shapes of these variables confirms this:
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
,因此在大多数情况下,最好使用后者.
The underlying problem is that exp_diag
is assumed to be either 1D or a column vector, but the diagonal of an np.matrix
object is a row vector. This highlights a more general point that np.matrix
is generally less well-supported than np.ndarray
, so in most cases it's better to use the latter.
一种可能的解决方案是使用np.ravel()
将diag_T
展平为一维np.ndarray
:
One possible solution would be to use np.ravel()
to flatten diag_T
into a 1D np.ndarray
:
diag_T = np.ravel(T.diagonal().copy())
这似乎可以解决您遇到的问题,尽管我可能还没有发现与np.matrix
有关的其他问题.
This seems to fix the problem you're encountering, although there may be other issues relating to np.matrix
that I haven't spotted yet.
我已经在此处打开了拉取请求.
I've opened a pull request here.
这篇关于Python中的矩阵求幂的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!