如何在涉及 matrix_power 的 matplotlib 中绘制矩阵的标量函数?

How to plot scalar functions of matrices in matplotlib involving matrix_power?

我有以下代码来绘制标量 x 与标量 f(x),其中函数内部有一些矩阵乘法:

import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import matrix_power
P=np.array([\
 [0,0,0.5,0,0.5],\
 [0,0,1,0,0], \
 [.25,.25,0,.25,.25], \
 [0,0,.5,0,.5], \
 [0,0,0,0,1], \
 ])
t=np.array([0,1,0,0,0])
ones=np.array([1,1,1,1,0])

def f(x):
    return t.dot(matrix_power(P,x)).dot(ones)
x=np.arange(1,20)
plt.plot(x, f(x))

现在,这个函数本身就可以正常工作了。

>>> f(1)
1.0
>>> f(2)
0.75 

但是绘图引发了错误 exponent must be an integer

换句话说,我如何在数组上计算这个函数?例如

f(np.array([1,2]))

我试着用 plt.plot(x, map(f,x))

但这并没有帮助。

我该如何解决这个问题?

In [1]: P=np.array([\ 
   ...:  [0,0,0.5,0,0.5],\ 
   ...:  [0,0,1,0,0], \ 
   ...:  [.25,.25,0,.25,.25], \ 
   ...:  [0,0,.5,0,.5], \ 
   ...:  [0,0,0,0,1], \ 
   ...:  ])                                                                     
In [2]:                                                                         
In [2]: P                                                                       
Out[2]: 
array([[0.  , 0.  , 0.5 , 0.  , 0.5 ],
       [0.  , 0.  , 1.  , 0.  , 0.  ],
       [0.25, 0.25, 0.  , 0.25, 0.25],
       [0.  , 0.  , 0.5 , 0.  , 0.5 ],
       [0.  , 0.  , 0.  , 0.  , 1.  ]])

In [4]: np.linalg.matrix_power(P,3)                                             
Out[4]: 
array([[0.   , 0.   , 0.25 , 0.   , 0.75 ],
       [0.   , 0.   , 0.5  , 0.   , 0.5  ],
       [0.125, 0.125, 0.   , 0.125, 0.625],
       [0.   , 0.   , 0.25 , 0.   , 0.75 ],
       [0.   , 0.   , 0.   , 0.   , 1.   ]])

In [5]: np.linalg.matrix_power(P,np.arange(0,4))                                
---------------------------------------------------------------------------    
TypeError: exponent must be an integer

所以只要给它它想要的整数:

In [10]: [f(i) for i in range(4)]                                               
Out[10]: [1.0, 1.0, 0.75, 0.5]

pylab.plot(np.arange(25), [f(i) for i in np.arange(25)]) 

来自matrix_power代码:

a = asanyarray(a)
_assertRankAtLeast2(a)
_assertNdSquareness(a)

try:
    n = operator.index(n)
except TypeError:
    raise TypeError("exponent must be an integer")
....

这是它为 n=3 所做的:

In [5]: x = np.arange(9).reshape(3,3)                                           
In [6]: np.linalg.matrix_power(x,3)                                             
Out[6]: 
array([[ 180,  234,  288],
       [ 558,  720,  882],
       [ 936, 1206, 1476]])
In [7]: x@x@x                                                                   
Out[7]: 
array([[ 180,  234,  288],
       [ 558,  720,  882],
       [ 936, 1206, 1476]])

您可以定义一个 matrix_power 接受幂数组的函数:

def matrix_power(P,x):
    return np.array([np.linalg.matrix_power(P,i) for i in x])

用这个 matrix_power(P,np.arange(25)) 会产生一个 (25,5,5) 数组。而你的 f(x) 实际上确实适用于它,返回一个 (25,) 形状数组。但我想知道,这是偶然的,还是故意的?您在编写 f 时是否考虑了 3d 幂数组?

t.dot(matrix_power(P,x)).dot(ones)