我有一个NxN矩阵,我对它的特征值分解感兴趣,以获得P = exp(tA(。为了评估P_ij(t(函数,我想准备一个E=NxNxN矩阵,其中E是v[i]和np.linalg.inv(v)[i]
的外积,v[i]
是特征向量之一。
我的 3x3 玩具矩阵的示例代码如下
def eig_trial(A,n):
v,w=np.linalg.eig(A)
w_inv=np.linalg.inv(
eigs=np.array([[[0 for k in xrange(len(v))] for j in xrange(n)] for i in xrange(n)])
for i in range(len(v)):
q=np.outer(np.array(w[:,i]),np.array(w_inv[i]))
eigs[i]=q
pprint.pprint(eigs)
A=np.asarray([[6,3,-2],[-4, -1, 2],[13,9,-3]])
eig_trial(A,3)
问题是,正确的值q=np.outer(...)
不等于eigs[i]
。
有什么原因不能将其存储到矩阵中吗?如果有帮助,我正在使用Canopy版本:2.1.9.3717(64位(和Python 2.7.13。
谢谢!
问题是您初始化的eigs
数组具有 dtypeint
In [96]: alist = ([[0 for k in range(5)] for j in range(3)])
In [97]: arr = np.array(alist)
In [98]: alist
Out[98]: [[0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0]]
In [99]: arr
Out[99]:
array([[0, 0, 0, 0, 0],
[0, 0, 0, 0, 0],
[0, 0, 0, 0, 0]])
并做一些作业:
In [100]: arr[0] = np.arange(5)
In [101]: alist[0] = np.arange(5)
In [102]: arr[1] = np.arange(5)*.1
In [103]: alist[1] = np.arange(5)*.1
arr
仍然是整数,即使最后一个赋值包含浮点数
In [104]: arr
Out[104]:
array([[0, 1, 2, 3, 4],
[0, 0, 0, 0, 0],
[0, 0, 0, 0, 0]])
In [105]: alist
Out[105]: [array([0, 1, 2, 3, 4]), array([0. , 0.1, 0.2, 0.3, 0.4]), [0, 0, 0, 0, 0]]
您可以通过以下方式更清楚地看到截断:
In [108]: arr[1] = np.arange(5)*1.3
In [109]: arr
Out[109]:
array([[0, 1, 2, 3, 4],
[0, 1, 2, 3, 5],
[0, 0, 0, 0, 0]])
如果数组以浮点型 dtype 开头:
In [110]: arr = np.zeros((3,5)) # float dtype default
In [111]: arr[0] = np.arange(5)
In [112]: arr[1] = np.arange(5)*1.3
In [113]: arr
Out[113]:
array([[0. , 1. , 2. , 3. , 4. ],
[0. , 1.3, 2.6, 3.9, 5.2],
[0. , 0. , 0. , 0. , 0. ]])
您也可以使用0.
而不是0
来初始化eigs
In [114]: np.array([0. for _ in range(5)])
Out[114]: array([0., 0., 0., 0., 0.])
但我上面使用zeros
更快(也更简单(
您无需分配列表的完整嵌套深度。只要第一层就足够了。 列表可以分配任何内容:
In [115]: alist = [None for _ in range(5)]
In [116]: alist
Out[116]: [None, None, None, None, None]
In [117]: alist[0] = np.arange(5)
In [118]: alist[1] = 32
In [119]: alist[2] = np.ones((3,2))
In [120]: alist
Out[120]:
[array([0, 1, 2, 3, 4]), 32, array([[1., 1.],
[1., 1.],
[1., 1.]]), None, None]
数组分配arr[0] = np.arange(5)
不同。 有一个隐含的arr[0,:,:] = ...
. 如果arr
中的选定插槽与 RHS 阵列不匹配,则会出现错误。match
表示形状完全匹配,或广播术语匹配。
当您遇到此类问题时,请尝试我演示的代码的小部分。 使用小型交互式示例跟踪问题要容易得多。
我真的不知道,但我发现了一个奇怪的解决方法。问题是由于 eig 是一个数组。如果它是列表列表,则它有效。
def eig_trial(A,n):
v,w=np.linalg.eig(A)
w_inv=np.linalg.inv(
eigs=[[[0 for k in xrange(len(v))] for j in xrange(n)] for i in xrange(n)] #np.array dropped
for i in range(len(v)):
q=np.outer(np.array(w[:,i]),np.array(w_inv[i]))
eigs[i]=q
eigs=np.asarray(eigs) #now convert it to an np array
pprint.pprint(eigs)
A=np.asarray([[6,3,-2],[-4, -1, 2],[13,9,-3]])
eig_trial(A,3)
如果有人解释原因,我会很高兴,但这至少是一个解决方案。