np.outer 输出在存储在数组中时会发生变化



我有一个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)

如果有人解释原因,我会很高兴,但这至少是一个解决方案。

最新更新