https://www.ibm.com/developerworks/community/blogs/jfp/entry/Elementary_Matrix_Operations_In_Python?lang=en


IT Best Kept Secret Is Optimization

Elementary Matrix Operations In Python
JeanFrancoisPuget | Dec 10 2015 | Visits (43902)
2 people like this
2

Octave andMatlab are high level languages that support vectors and matrices with a very simple syntax.  Python support for matrices is not as nice, but few little tricks should do the job. 

Let me first briefly introduce how Octave and Matlab support elementary matrices operations, then we'll look at how to achieve the same with Python.

The following is based on Octave tutorial.  The code runs fine with Matlab.

Octave and Matlab

Vector and matrices can be created from a list of elements.  Commas are used to separate elements in the same row, while semicolons are used to separate rows:

Vector creation:

>> x = [1; 3; 2]
x =

   1
   3
   2

Transposed vector (horizontal vector) can be created directly with elements separated by a comma, or obtained was the transpose of a vector using the' operator.

>> y = [1, 3, 2]
y =

   1   3   2

>> x'
ans =

   1   3   2

Matrix creation mixes the use of comma and semicolon:

>> A = [1, 1, 2; 3, 5, 8; 13, 21, 34]
A =

    1    1    2
    3    5    8
   13   21   34

There are a number of functions that return a matrix.  For instance, creating a zeroed matrix:

>> zeros(3,4)
ans =

   0   0   0   0
   0   0   0   0
   0   0   0   0

Matrix dimensions are retrieved with the size function.  Note that vectors are 2D matrices actually:

>> size(A)
ans =

   3   3

>> size(x)
ans =

   3   1

>> size(y)
ans =

   1   3

It can be convenient to create vectors as an arithmetic series.  We first create a vector, then reshape it to get a 2 row matrix:

>> B = 0:3:15
B =

    0    3    6    9   12   15

>> B = reshape(B,3,2)
B =

    0    9
    3   12
    6   15

Arithmetic operations:  +,-,*,/, can involve a scalar and a matrix.

>> A*2
ans =

    2    2    4
    6   10   16
   26   42   68

>> B/3
ans =

   0   3
   1   4
   2   5

Matrix multiplication

>> A*x
ans =

     8
    34
   144

>> y*A
ans =

   36   58   94

>> A*B
ans =

    15    51
    63   207
   267   879

Element wise operations are prefixed with a dot, while operations without prefix deal with matrices.  For instance, for exponentiation

>> A.^2
ans =

      1      1      4
      9     25     64
    169    441   1156

>> A^2
ans =

     30     48     78
    122    196    318
    518    832   1350

Let's conclude with vector products:

>> z = [3;4;5]
z =

   3
   4
   5

Inner product

>> x' * z
ans =  25

Outer product

>> x * z'
ans =

    3    4    5
    9   12   15
    6    8   10

Element wise product

>> x.*z
ans =

    3
   12
   10

Python

For technical computing, I recommend the use of Numpy arrays instead of the native Python arrays.  Indeed, Numpy is used by most scientific packages in Python, including Pandas, Scipy, and Scikit-Learn.  Numpy provides a matrix class that can be used to mimic Octave and Matlab operations.  For instance, one can create matrices using a similar syntax:

>>> import numpy as np
>>> np.matrix('1 2; 3 4')
matrix([[1, 2],
        [3, 4]]
)

However most Python scientific functions deal with 2D arrays instead of matrices.  Moving back and forth from arrays to matrices is easy, but it slows the code.  Good news is that most matrix operations can be used with 2D Numpy arrays.  Let's see how, by replicating the above Octave/Matlab examples with Numpy arrays.

Vector creation.  Vectors are one dimensional arrays in Numpy.

>>> y = np.array([1,3,2])
>>> y
array([1, 3, 2])

This creates a transposed vector (horizontal vector).  Vector, array, and matrix dimension can be retrieved using theshape attribute.

>>> y.shape
(3,)

We can turn this into a 2D matrix by adding a row  dimension to it:

>>> y = y[None,:]
array([[1, 3, 2]])

>>> y.shape
(1, 3)

Vectors (vertical vectors) can be created in a similar way, but the added dimension should be the column dimension this time

>>> x = np.array([1,3,2])[:,None]
>>> x
array([[1],
       [3],
       [2]])
>>> x.shape
(3, 1)

Matrix creation can be done in two steps, first create a 1D array, then reshape it:

>>> A = np.array([1,1,2,3,5,8,13,21,34]).reshape(3,3)
>>> A
array([[ 1,  1,  2],
       [ 3,  5,  8],
       [13, 21, 34]])

>>> A.shape
(3, 3)

There are a number of array creations functions that return a 2D array directly.  For instance, creating a zeroed 2D array of integers:

>>> np.zeros((3,4),dtype=int)
array([[0, 0, 0, 0],
       [0, 0, 0, 0],
       [0, 0, 0, 0]])

It can be convenient to create vectors as an arithmetic series.  We first create a vector, then reshape it to get a 2 row matrix:

>>> B = np.arange(0,16,3)
>>> B
array([ 0,  3,  6,  9, 12, 15])
>>> B = B/3B.reshape(3,2)
array([[ 0,  3],
       [ 6,  9],
       [12, 15]])

Note that the arguments of the arange function are different from the rang expression in Octave/Matlab.  The stop isn't included in Python and it comes in second position rather than third.

The order of elements in B is also different from the Octave/Matlab case. Python uses C ordering by default, while Octave/Matlab use Fortran ordering.  If we want the same order as in Octave/Matlab, then we have to specify a Fortran ordering.

>>> np.arange(0,16,3).reshape((3,2), order='F')
array([[ 0,  9],
       [ 3, 12],
       [ 6, 15]])

Arithmetic operations:  +,-,*,/, can involve a scalar and a matrix.

>>> A*2
array([[ 2,  2,  4],
       [ 6, 10, 16],
       [26, 42, 68]])

>>> B/3
array([[0, 1],
       [2, 3],
       [4, 5]])

Matrix multiplication is the dot() function.  That's one of the case where Numpy matrices differ from Numpy arrays.  With Numpy matrices, * is the matrix multiplication.

>>> A*x
array([[ 1,  1,  2],
       [ 9, 15, 24],
       [26, 42, 68]])
>>> y
array([1, 3, 2])
>>> np.dot(A,x)
array([[  8],
       [ 34],
       [144]])
>>> np.dot(y,A)
array([36, 58, 94])
>>> np.dot(A,B)
array([[ 30,  42],
       [126, 174],
       [534, 738]])

Note that * is the element wise operation in Python.  Matrix multiplication

In general, Numpy operations are element wise operations.  Example for exponentiation:

>>> A**2
array([[   1,    1,    4],
       [   9,   25,   64],
       [ 169,  441, 1156]])

Squaring a matrix must use the dot function:

>>> np.dot(A,A)
array([[  30,   48,   78],
       [ 122,  196,  318],
       [ 518,  832, 1350]])

Let's conclude with vector products:

>>> z = np.array([3,4,5])[:,None]
>>> z
array([[3],
       [4],
       [5]])

Inner product can be obtained via different ways, simplest is probably using sum:np.sum(x*z)

>>> np.sum(x*z)
25

or

>>> (x*z).sum()
25

Outer product

>>> np.dot(x,z.T)
array([[ 3,  4,  5],
       [ 9, 12, 15],
       [ 6,  8, 10]])

Element wise product

>>> x*z
array([[ 3],
       [12],
       [10]])

Numpy Broadasting

The above assumed we systematically reshape 1D arrays into 2D arrays.

This is not required in general thanks to Numpy broadcasting rules.  Indeed, when an operation is applied between two arrays of differing dimensions, Numpy will automatically expand the smallest one by adding dimensions in front of it.  A 1D array of length n will be automatically expanded into a1xn 2D array if need be.  Here is an example.

>>> t = np.array([3,4,5])
>>> t
array([3, 4, 5])
>>> B
array([[ 0,  3],
       [ 6,  9],
       [12, 15]])
>>> np.dot(t,B)
array([ 84, 120])

We did not have to transform t into a 2D matrix via t[None,:] .

Note however that Numpy broadcasting adds dimensions in front of the array, not at the rear.  Therefore, we need to explicitly add the dimension when we want to use a 1D array of lengthm as a mx1 2D array.  For instance, here is the multiplication of A by T as a vertical vector:

>>> np.dot(A,t[:,None])
array([[ 17],
       [ 69],
       [293]])

 

Careful use of Numpy broadcasting and reshaping is enough to implement usual matrix operations on arrays.  The syntax isn't nearly as elegant as with Octave/Matlab, but it is usable, and more efficient than using the matrix class.

 

There is an advantage of Numpy however: arrays can have more than 2 dimensions.  The above tricks can be used with higher dimension operations.  The only thing to really care about is to explicitly add dimensions where they are needed when the default broadcasting rule does not fit the need. 



Tags:  pythonmatlaboctavenumpy


Logo

开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!

更多推荐