Programming with Matrices

Programming with Matrices (Numpy and Eigen)

Numpy and Eigen are libraries for matrix calculations in Python and C++, respectively. We will be doing a lot of matrix calculations to calculate probabilities shortly, but before that, we do a session on Numpy and Eigen.

Numpy

We start by importing numpy

import numpy.linalg
import numpy as np

We create some containers containing matrices and vectors. A matrix is a list of rows.

a = np.array([1,2,3])
b = np.array([[1,2,3],[4,5,6]])
c = [[1,2,3],[4,5,6]] d = np.array(c)

We can specify the type of the entries by using the parameter dtype. For this course, float will be our default value.

c = np.array([1,2,3],dtype=float) d = np.array([[1,2,3],[4,5,6]],dtype=float)

Let's see how to query the shape of a matrix, and change it.

print(b.shape)
b.reshape(3,2)
print(b)

You can't reshape if the total number of elements changes.

We will be working with large matrices, which we will fill with entries from files, or programs that we write. So it can be useful to create a matrix filled with zeros. We will also need the identity matrix.

a = np.zeros(10,dtype=float)
b = np.zeros((10,10),dtype=float)
unit = np.identity(5, dtype=float)

You can access and change elements in a matrix

a[1,2] = 10.0
print(a[1,2])

The format is a[i,j] where i is the row-number and j is the column number. Note that, unlike in mathematics, column and row indices start at 0.

Numpy arrays can be accessed and sliced just like Python containers. These operations are for instance useful if you want to pick out a submatrix of a larger matrix.

a = np.array([[1,2,3],[4,5,6],[7,8,9]])

# slicing rows
print(a[1:])
print(a[:2])
print(a[1:3])

# slicing columns
print(a[:,1:])
print(a[:,:2])
print(a[:,1:3])

# slicing both
print(a[1:,1:])
print(a[:2,:2])
print(a[:3,1:3])

Slicing often creates what is called a view of the larger matrix. That means that changing an entry can have effect on both the original and the slice. To avoid dealing with views (which is most relevant for large matrices), you can add

a = np.copy(b)

to ensure that you are creating a copy of your matrix

Let's look at some matrix operations which will be useful for this module. We start with the transpose, we can be expressed in two ways

a = np.array([[1,2,3],[4,5,6],[7,8,9]], dtype=float)
print(a)
# transpose of a
b = a.T
c = a.transpose()

More advanced matrix operations can be found in the module numpy.linalg. We can for instance do addition, subtraction, scalar multiplication, matrix-inverse and matrix multiplication by

# some matrix computations
a = np.array([[1,0,2],[0,5,0],[1,2,3]], dtype=float)
b = np.identity(3)

# compute the inverse of a
inva = np.linalg.inv(a)

# compute the matrix product og a and inva
print(np.dot(a,inva))

# compute the matrix product of a and b
print(np.dot(a,b))

# add a and b
print(np.add(a,b))

# subtract b from a
print(np.subtract(a,b))
# scalar multiplication
print(np.multiply(2.0,a))

The matrix product of a and inva is supposed to compute to the identity matrix, but doesn't. Can you come up with an explanation why?

Eigen

Unlike SFML we do not need to download object files when we install Eigen. Nor do we need to precompile those. The simplest installation procedure would be to install the whole of Eigen where you keep your other personal .hpp files. As long as you make sure that your project is set to search in that directory, you should be good to go.

We now do a short Eigen tutorial similar to the numpy tutorial above. First

MatrixXd M1(4, 4);
M1.setZero();
std::cout << M1 << std::endl << std::endl;

creates a 4x4 matrix called M1 which is set to be the zero-matrix. The type of the entries of this matrix is double, that is all we need to know for now. I refer you to the documentation for more information about various other matrix-types implemented in Eigen. We see that std::cout works fine with matrices, which is good news.

There are several methods of initialising matrices with values. For example we can use random values, or we can fill in all entries with the same value.

MatrixXd M3 = MatrixXd::Random(2, 3);
MatrixXd M4 = MatrixXd::Constant(3, 3, 1.2);

Eigen implements an interface similar to iostream and string for working with matrix enties. For example, the 2x2 matrix M2 can get new values by

M2 << 1, 2, 3, 4;

This creates the matrix $$M2=\begin{pmatrix}1 & 2 \\ 3 & 4\end{pmatrix}.$$ Entries are filled in row by row.

Now a 4x4 matrix like M1 can be thought of as a 2x2 matrix where each entry is also a 2x2 matrix. So the following works fine

M1 << M2, M2, M2, M2;

We can fill in individual rows and columns

M2.row(1) << 0, 0;
M2.col(0) << 10, 10;

Note that row and column indexes start at 0 and that this is not what we do in mathematics. We can of course access individual elements, for example by

M3(1, 1) = 12.0;

There are several ways to look at or change a whole submatrix (similar to slicing in numpy). Here are some examples.

// We pick out the 2x2 submatrix with
// indices (0,0) -- (2,2) and set it to be zero
M1.topLeftCorner(2, 2).setZero();
std::cout << M1 << std::endl << std::endl;

// We pick out the 2x2 submatrix with indices
// (2,2) -- (4,4) and set it to the identity matrix
M1.bottomRightCorner(2, 2).setIdentity();
std::cout << M1 << std::endl << std::endl;

// We pick out a 2x2 block with top-left corner at (1,1)
M1.block<2, 2> (1, 1) << 10, 11, 12, 13;
std::cout << M1 << std::endl << std::endl;

// We pick out a block with top-left corner (1,1)
//and bottom right corner (2,2)
M1.block(1, 1, 2, 2) << 20, 21, 22, 23;
std::cout << M1 << std::endl << std::endl;

// We fill in the first 1 rows
M1.topRows(1) << 1, 2, 3, 4;
std::cout << M1 << std::endl << std::endl;

// We fill in two leftmost columns. We need 8 numbers
// for two columns in M1
M1.leftCols(2) << 1, 2, 3, 4, 5, 6, 7, 8;
std::cout << M1 << std::endl << std::endl;

That should be more than enough information about accessing matrices. We will now look at some matrix operations used in this module. We create three random matrices and perform some self-explanatory operations on them.

M1 = MatrixXd::Random(2, 2);
M2 = MatrixXd::Random(2, 3);
M3 = MatrixXd::Random(2, 2);
std::cout << M1 * M2 ;
std::cout << M1 + M3 ;
std::cout << M1 * 2 + 2 * M3;
// std::cout << M2 * M1; // FAILS - do you see why?

Two more operations which could be of use in this module is the transpose and the inverse.

std::cout << M1.inverse() << std::endl;
std::cout << M1.transpose();

We will look at more matrix operations when we need to later in the course.