In this lab sheet we will learn how to use Python to study linear algebra.

Building blocks

  1. TICKABLE Linear algebra is the study of systems of linear equations.

    A video describing the concept.

    A video demo.

    Many mathematical problems reduce to solving systems of linear equations. Here is one such system that we will use as a running example:

    We can represent this system of equations as a matrix equation (if you are not familiar with this you will learn more about it soon):

    Where \(A\) corresponds to the coefficients of the equations above, and \(b\) to the right hand side:

    Solving a system of linear equations can be done very efficiently when using this matrix notation. There are various Python libraries that are able to do this. We will use NumPy which is a library that specialises in efficient numerical manipulation (it does a lot more than just solve linear equations). Let us set up the above matrix equation:

    >>> import numpy as np
    >>> A = np.matrix([[5, 1, -1], [-1, 2, 4], [1, 1, 1]])
    >>> A
    matrix([[ 5,  1, -1],
            [-1,  2,  4],
            [ 1,  1,  1]])
    >>> b = np.matrix([[0], [-2], [4]])
    >>> b
    matrix([[ 0],
            [-2],
            [ 4]])
    
    

    Once we have these NumPy objects we can solve our equation using the .linalg.solve command:

    >>> np.linalg.solve(A, b)
    matrix([[-14.],
            [ 44.],
            [-26.]])
    
    

    Experiment with another linear system.

  2. TICKABLE Manipulating matrices.

    A video describing the concept.

    A video demo.

    It is possible to add two matrices together, this corresponds to adding each corresponding element:

    >>> B = np.matrix([[1, 2, 0], [-4, 2, 2], [1, 3, 1]])
    >>> A + B
    matrix([[ 6,  3, -1],
            [-5,  4,  6],
            [ 2,  4,  2]])
    
    

    Similarly we can carry out scalar multiplication of matrices:

    >>> 5 * A
    matrix([[25,  5, -5],
            [-5, 10, 20],
            [ 5,  5,  5]])
    
    

    Finally we can also carry out matrix multiplication:

    >>> A * B
    matrix([[ 0,  9,  1],
            [-5, 14,  8],
            [-2,  7,  3]])
    
    

    which in turn implies we can take the power of a matrix. Here is \(A ^ 3 = AAA\):

    >>> A ** 3
    matrix([[107,  33,  -1],
            [ -9,  24,  44],
            [ 25,  17,  15]])
    
    

    Modify the code here to try multiplying different matrices together.

  3. TICKABLE Identity matrix.

    A video describing the concept.

    A video demo.

    As we have matrix multiplication, that implies that there will be an element that does not have any effect when multiplying by it. (In normal multiplication this is the number \(1\) as \(1 \times 5 = 5\).

    This is referred to as the identity matrix and corresponds to a diagonal of \(1\)s. Here, for example is the identity matrix of size 3:

    Here is the same matrix in python:

    >>> I = np.matrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
    >>> I
    matrix([[1, 0, 0],
            [0, 1, 0],
            [0, 0, 1]])
    >>> I * A == A  # We see that all elements are equal
    matrix([[ True,  True,  True],
            [ True,  True,  True],
            [ True,  True,  True]], dtype=bool)
    >>> np.array_equal(I * A, A)  # We see that the matrices are equal
    True
    
    

    As an inverse matrix of a given size is often used Numpy has a shorthand to get one:

    >>> I = np.identity(3)
    >>> I  # This is a numpy array
    array([[ 1.,  0.,  0.],
           [ 0.,  1.,  0.],
           [ 0.,  0.,  1.]])
    >>> I = np.matrix(I)
    >>> I
    matrix([[ 1.,  0.,  0.],
            [ 0.,  1.,  0.],
            [ 0.,  0.,  1.]])
    
    

    Experiment with the above with different matrices.

  4. TICKABLE Matrix inversion.

    A video describing the concept.

    A video demo.

    Related to the idea of the identity matrix is the notion of the inverse of a matrix \(A^{-1}\) so that:

    Here is how we compute the matrix inverse in Python:

    >>> Ainv = np.linalg.inv(A)
    >>> Ainv
    matrix([[ 1. ,  1. , -3. ],
            [-2.5, -3. ,  9.5],
            [ 1.5,  2. , -5.5]])
    
    

    We see that multiplying this inverse with the original matrix gives the identity matrix:

    >>> Ainv * A
    matrix([[  1.0...
    
    

    Solving systems of linear equations corresponds to inverting a matrix:

    This can be checked here:

    >>> Ainv * b  # We obtain the same solution as before
    matrix([[-14.],
            [ 44.],
            [-26.]])
    
    

    Experiment with the above.

  5. TICKABLE

    A video describing the concept.

    A video demo.

    Related to the inverse of a matrix (it is used for its calculation) is called the determinant. We won’t go in to detail about that here but this is how it is calculated:

    >>> np.linalg.det(A)
    -2.0...
    
    

    For the various matrices you have experimented with, calculate the determinant.

  6. TICKABLE: Worked example: Large linear system.

    A video describing the concept.

    A video demo.

    Let us consider the following hypothetical problem. A set of \(N\) business people who agree together to invest in a company. The first business person and the second business person agree to contribute a combined £1, the second and the third a combined £2, the third and the fourth a combined £3, and so on and so on until the \(N-1\)th and the \(N\)th a combined £\(N-1\). Finally, the \(N\)th and the first business person contribute a combined £\(N\).

    Assume that \(N=5\), how much will each person contribute to the deal, and subsequently what is the total amount?

    If we let \(x_i\) represent the amount contributed by business person \(i\), then the above corresponds to:

    Let us build the matrix that corresponds to this linear system. First the coefficients of the equations:

    >>> N = 5
    >>> A = np.matrix(np.zeros((N, N)))  # An N by N array of zeros
    >>> for deal in range(N):  # Putting the 1s in the right position
    ...     A[deal, deal] = 1
    ...     A[deal, (deal + 1) % N] = 1
    >>> A
    matrix([[ 1.,  1.,  0.,  0.,  0.],
            [ 0.,  1.,  1.,  0.,  0.],
            [ 0.,  0.,  1.,  1.,  0.],
            [ 0.,  0.,  0.,  1.,  1.],
            [ 1.,  0.,  0.,  0.,  1.]])
    
    

    Now for the right hand side:

    >>> b = np.matrix([[i + 1] for i in range(N)])
    >>> b
    matrix([[1],
            [2],
            [3],
            [4],
            [5]])
    
    

    We can now easily find out the contributions from each business person:

    >>> contributions = np.linalg.solve(A, b)
    >>> contributions
    matrix([[ 1.5],
            [-0.5],
            [ 2.5],
            [ 0.5],
            [ 3.5]])
    >>> sum(contributions)  # Total contributions
    matrix([[ 7.5]])
    
    

    We see that the total contribution is £7.5, but that the second contributor actually gains money from the deal! (A negative contribution.)

    Further work

    These questions aim to push a bit further.

  7. Consider a similar situation to that described in question 6 except that business person \(i\) and half of the contribution of business person \(i+1\) always equal to £5.

    (For the last pairing we have that the contribution of business person \(N\) and half of that of business person \(0\) is £5).

    What is the contribution of each player and what is the total contribution for \(N=50\)?

  8. Modify the problem of question 7 with the assumption that the contributions are now made by triplets of players so that: business person \(i\) and half of the contribution of business person \(i+1\), and a third of the contributions of business person \(i+2\) is always equal to £5.

  9. Further aspects of linear algebra are called eigenvalues and eigenvectors. Read about them and see how to compute them in numpy.

  10. Scipy (another package) and Sympy can also be used for linear algebra. Explore those libraries and compare them to numpy.

Further resources

Solutions

Solutions available.