Vectorization and matrix multiplication by scalars

Question

I am new to python/numpy. I need to do the following calculation: for an array of discrete times t, calculate $e^{At}$ for a $2\times 2$ matrix $A$

What I did:

def calculate(t_,x_0,v_0,omega_0,c):

    # define A
    a_11,a_12, a_21, a_22=0,1,-omega_0^2,-c
    A =np.matrix([[a_11,a_12], [a_21, a_22]]) 
    print A
    # use vectorization 
    temps = np.array(t_)
    A_ = np.array([A  for k in  range (1,n+1,1)])
    temps*A_
    x_=scipy.linalg.expm(temps*A)
    v_=A*scipy.linalg.expm(temps*A)
    return x_,v_

n=10
omega_0=1
c=1
x_0=1
v_0=1
t_ = [float(5*k*np.pi/n)  for k in  range (1,n+1,1)]
x_, v_ = calculate(t_,x_0,v_0,omega_0,c)

However, I get this error when multiplying A_ (array containing n times A ) and temps (containg the times for which I want to calculate exp(At) :

ValueError: operands could not be broadcast together with shapes (10,) (10,2,2)

As I understand vectorization, each element in A_ would be multiplied by element at the same index from temps; but I think i don't get it right. Any help/ comments much appreciated

Aside as I just noticed it; is there any particular reason to use `np.matrix`? The Numpy creators basically wish it didn't exist, so it's generally better to use a multidimensional `np.array` except for a very select number of operations. — roganjosh, Feb 05 '18 at 22:14
@Azevedo I added the exponential. I didn't put it in the question because I wanted to isolate the problem — Conjecture, Feb 06 '18 at 07:54
@roganjosh, the purpose of this line is to do vectorizarion. I thought I had to provide 2 numpy arrays of the same number of elements to be able to do vectorized multiplication. There was no particular reason for using np.matrix, I just used it because It was the first solution I found as I am new to numpy — Conjecture, Feb 06 '18 at 07:55

hpaulj · Accepted Answer · 2018-02-06T07:25:29.183

A pure numpy calculation of t_ is (creates an array instead of a list):

In [254]: t = 5*np.arange(1,n+1)*np.pi/n
In [255]: t
Out[255]: 
array([ 1.57079633,  3.14159265,  4.71238898,  6.28318531,  7.85398163,
        9.42477796, 10.99557429, 12.56637061, 14.13716694, 15.70796327])

In [256]: a_11,a_12, a_21, a_22=0,1,-omega_0^2,-c
In [257]: a_11
Out[257]: 0
In [258]: A = np.array([[a_11,a_12], [a_21, a_22]]) 
In [259]: A
Out[259]: 
array([[ 0,  1],
       [-3, -1]])
In [260]: t.shape
Out[260]: (10,)
In [261]: A.shape
Out[261]: (2, 2)
In [262]: A_ = np.array([A  for k in  range (1,n+1,1)])
In [263]: A_.shape
Out[263]: (10, 2, 2)

A_ is np.ndarray. I made A a np.ndarray as well; yours is np.matrix, but your A_ will still be np.ndarray. np.matrix can only be 2d, where as A_ is 3d.

So t * A will be array elementwise multiplication, hence the broadcasting error, (10,) (10,2,2).

To do that elementwise multiplication right you need something like

In [264]: result = t[:,None,None]*A[None,:,:]
In [265]: result.shape
Out[265]: (10, 2, 2)

But if you want matrix multiplication of the (10,) with (10,2,2), then einsum does it easily:

In [266]: result1 = np.einsum('i,ijk', t, A_)
In [267]: result1
Out[267]: 
array([[   0.        ,   86.39379797],
       [-259.18139392,  -86.39379797]])

np.dot can't do it because its rule is 'last with 2nd to last'. tensordot can, but I'm more comfortable with einsum.

But that einsum expression makes it obvious (to me) that I can get the same thing from the elementwise *, by summing on the 1st axis:

In [268]: (t[:,None,None]*A[None,:,:]).sum(axis=0)
Out[268]: 
array([[   0.        ,   86.39379797],
       [-259.18139392,  -86.39379797]])

Or (t[:,None,None]*A[None,:,:]).cumsum(axis=0) to get a 2x2 for each time.

You could use dot with rollaxis. Also, minor nitpick, but I'm pretty sure the class is called np.ndarray, not np.array. +1 — Mad Physicist, Feb 06 '18 at 05:51

score 0 · Answer 2 · answered Feb 05 '18 at 23:22

This is what I would do.

import numpy as np
from scipy.linalg import expm
A = np.array([[1, 2], [3, 4]])
for t in np.linspace(0, 5*np.pi, 20):
    print(expm(t*A))

No attempt at vectorization here. The expm function applies to one matrix at a time, and it surely takes the bulk of computation time. No need to worry about the cost of multiplying A by a scalar.

Vectorization and matrix multiplication by scalars

2 Answers2

Linked