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

Question

I have the following code to plot scalar x vs scalar f(x) where there is some matrix multiplication inside the function:

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))

Now, the function by itself works fine.

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

But the plotting raises the error exponent must be an integer.

To put it another way, how do I evaluate this function upon an array? e.g.

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

I tried replacing the plot line with plt.plot(x, map(f,x))

But this didn't help.

How can I fix this?

Don't "omit something for simplicity"! – ImportanceOfBeingErnest Mar 04 '19 at 01:47 — ImportanceOfBeingErnest, Mar 04 '19 at 01:47

hpaulj · Accepted Answer · 2019-03-04T19:10:52.407

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

So just give it the integer that it wants:

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)])

From the matrix_power code:

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

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

Here's what it does for 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]])

You could define a matrix_power function that accepts an array of powers:

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

With this matrix_power(P,np.arange(25)) would produce a (25,5,5) array. And your f(x) actually does work with that, returning a (25,) shape array. But I wonder, was that just fortuitous, or was it intentional? Did you write f with a 3d power array in mind?

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

can you show how to use this in the `plt.plot` function, as a function, rather than a list of values? I think that is what I really need. — Tim, Mar 04 '19 at 03:19
Maybe I'm not clear: Can you use a function, **not** a comprehension list of values? In other words, I want the paradigm of `x=np.arange(25)\ plt.plot(x,np.sin(x))`? ie. `sin(x)` seems to work without requiring me to make a list comprehension.My function doesn't. I want to make my function work as simply as the sin(x) paradigm. Your answer seems 'over-specified' since it requires the range to be specified twice — Tim, Mar 04 '19 at 06:07
And the real problem "exponent must be an integer" remains unexplained. — Tim, Mar 04 '19 at 06:13
@Tim ` np.array` returns an array, not an int (seemingly `matrix_power` can't handle arrays). Also look [here](https://stackoverflow.com/questions/35215161/most-efficient-way-to-map-function-over-numpy-array) — user8408080, Mar 04 '19 at 11:05
`np.linalg.matrix_power(x,3)` does `x@x@x`, that is, repeated matrix products. That only makes sense for an integer value. Right at the start it checks that `n` is a scalar integer. It is not coded to work with a list or array of powers. You get to do that, as I illustrated. — hpaulj, Mar 04 '19 at 18:16
'seemingly `matrix_power` can't handle arrays' : @user8408080 is the crux for me. I expected that numpy's functions always know how to handle arrays. But only some, apparently. — Tim, Mar 04 '19 at 21:06
`matrix_power` does know how to handle an array - the first argument. But the docs are quite explicit - the second argument is an integer. — hpaulj, Mar 04 '19 at 21:16

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

1 Answers1