Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Lecture 01: Python primer

University of Lausanne

Course: Deep Learning for Solving and Estimating Dynamic Models in Economics and Finance
Script reference: Front matter and Appendix~E — NumPy arrays, broadcasting, linear algebra
Notebook role: primer (pre-course self-study; skip if you write Python every day)
Author: Simon Scheidegger

Python Basics 9: Numpy Arrays and Linear Algebra

The ability to do basic things with numerical arrays (vectors and matrices) is extremely important in data science applications. We do numerical algebra in Python using the Numpy. Numpy is a huge Python library and it has a lot of details. We are going to cover the very basics here. When there is something specific that you want to do, the best way to figure out how to do it is to Google it.

First, let’s import numpy. We typically, do it by creating a shortcut called np:

import numpy as np

In this way you can access whatever is in numpy by going throuh np which is less characters. Let’s create a 1D numpy array from a list:

a = np.array([1.2, 2.0, 3.0, -1.0, 2.0])
print(a)
[ 1.2  2.   3.  -1.   2. ]

This behaves just like a vector in Matlab. You can multiply with a number:

2.0 * a
array([ 2.4, 4. , 6. , -2. , 4. ])
a * 3
array([ 3.6, 6. , 9. , -3. , 6. ])

You can add another vector (of the same size):

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

You can raise all elements to a power:

a ** 3
array([ 1.728, 8. , 27. , -1. , 8. ])

Or you can apply a standard function to all elements:

np.exp(a)
array([ 3.32011692, 7.3890561 , 20.08553692, 0.36787944, 7.3890561 ])
np.cos(a)
array([ 0.36235775, -0.41614684, -0.9899925 , 0.54030231, -0.41614684])
np.sin(a)
array([ 0.93203909, 0.90929743, 0.14112001, -0.84147098, 0.90929743])
np.tan(a)
array([ 2.57215162, -2.18503986, -0.14254654, -1.55740772, -2.18503986])

Notice that the functions you can use are in np. Here is how you can get how many elements you have in the array:

a.shape
(5,)

Also, to see if this is a 1D array (higher dimensional arrays, 2D, 3D, etc. are also possible), you can use this:

a.ndim
1

You access the elements of the array just like the elements of a list:

a[0]
1.2
a[1]
2.0
a[2:4]
array([ 3., -1.])
a[::2]
array([1.2, 3. , 2. ])
a[::-1]
array([ 2. , -1. , 3. , 2. , 1.2])

You can find the minimum or the maximum of the array like this:

a.min()
-1.0
a.max()
3.0

So, notice that a is an object of a class. The class is called np.ndarray:

type(a)
numpy.ndarray

and min() and max() are functions of the np.ndarray class. Another function of the class is the dot product:

a.dot(b)
20.2

There is also a submodule called numpy.linalg which has some linear algebra functions you can apply to arrays. For 1D arrays (aka vectors) the vector norm is a useful one:

np.linalg.norm(a)
4.409081537009721
help(np.linalg.norm)
Help on function norm in module numpy.linalg:

norm(x, ord=None, axis=None, keepdims=False)
    Matrix or vector norm.
    
    This function is able to return one of eight different matrix norms,
    or one of an infinite number of vector norms (described below), depending
    on the value of the ``ord`` parameter.
    
    Parameters
    ----------
    x : array_like
        Input array.  If `axis` is None, `x` must be 1-D or 2-D, unless `ord`
        is None. If both `axis` and `ord` are None, the 2-norm of
        ``x.ravel`` will be returned.
    ord : {non-zero int, inf, -inf, 'fro', 'nuc'}, optional
        Order of the norm (see table under ``Notes``). inf means numpy's
        `inf` object. The default is None.
    axis : {None, int, 2-tuple of ints}, optional.
        If `axis` is an integer, it specifies the axis of `x` along which to
        compute the vector norms.  If `axis` is a 2-tuple, it specifies the
        axes that hold 2-D matrices, and the matrix norms of these matrices
        are computed.  If `axis` is None then either a vector norm (when `x`
        is 1-D) or a matrix norm (when `x` is 2-D) is returned. The default
        is None.
    
        .. versionadded:: 1.8.0
    
    keepdims : bool, optional
        If this is set to True, the axes which are normed over are left in the
        result as dimensions with size one.  With this option the result will
        broadcast correctly against the original `x`.
    
        .. versionadded:: 1.10.0
    
    Returns
    -------
    n : float or ndarray
        Norm of the matrix or vector(s).
    
    See Also
    --------
    scipy.linalg.norm : Similar function in SciPy.
    
    Notes
    -----
    For values of ``ord < 1``, the result is, strictly speaking, not a
    mathematical 'norm', but it may still be useful for various numerical
    purposes.
    
    The following norms can be calculated:
    
    =====  ============================  ==========================
    ord    norm for matrices             norm for vectors
    =====  ============================  ==========================
    None   Frobenius norm                2-norm
    'fro'  Frobenius norm                --
    'nuc'  nuclear norm                  --
    inf    max(sum(abs(x), axis=1))      max(abs(x))
    -inf   min(sum(abs(x), axis=1))      min(abs(x))
    0      --                            sum(x != 0)
    1      max(sum(abs(x), axis=0))      as below
    -1     min(sum(abs(x), axis=0))      as below
    2      2-norm (largest sing. value)  as below
    -2     smallest singular value       as below
    other  --                            sum(abs(x)**ord)**(1./ord)
    =====  ============================  ==========================
    
    The Frobenius norm is given by [1]_:
    
        :math:`||A||_F = [\sum_{i,j} abs(a_{i,j})^2]^{1/2}`
    
    The nuclear norm is the sum of the singular values.
    
    Both the Frobenius and nuclear norm orders are only defined for
    matrices and raise a ValueError when ``x.ndim != 2``.
    
    References
    ----------
    .. [1] G. H. Golub and C. F. Van Loan, *Matrix Computations*,
           Baltimore, MD, Johns Hopkins University Press, 1985, pg. 15
    
    Examples
    --------
    >>> from numpy import linalg as LA
    >>> a = np.arange(9) - 4
    >>> a
    array([-4, -3, -2, ...,  2,  3,  4])
    >>> b = a.reshape((3, 3))
    >>> b
    array([[-4, -3, -2],
           [-1,  0,  1],
           [ 2,  3,  4]])
    
    >>> LA.norm(a)
    7.745966692414834
    >>> LA.norm(b)
    7.745966692414834
    >>> LA.norm(b, 'fro')
    7.745966692414834
    >>> LA.norm(a, np.inf)
    4.0
    >>> LA.norm(b, np.inf)
    9.0
    >>> LA.norm(a, -np.inf)
    0.0
    >>> LA.norm(b, -np.inf)
    2.0
    
    >>> LA.norm(a, 1)
    20.0
    >>> LA.norm(b, 1)
    7.0
    >>> LA.norm(a, -1)
    -4.6566128774142013e-010
    >>> LA.norm(b, -1)
    6.0
    >>> LA.norm(a, 2)
    7.745966692414834
    >>> LA.norm(b, 2)
    7.3484692283495345
    
    >>> LA.norm(a, -2)
    0.0
    >>> LA.norm(b, -2)
    1.8570331885190563e-016 # may vary
    >>> LA.norm(a, 3)
    5.8480354764257312 # may vary
    >>> LA.norm(a, -3)
    0.0
    
    Using the `axis` argument to compute vector norms:
    
    >>> c = np.array([[ 1, 2, 3],
    ...               [-1, 1, 4]])
    >>> LA.norm(c, axis=0)
    array([ 1.41421356,  2.23606798,  5.        ])
    >>> LA.norm(c, axis=1)
    array([ 3.74165739,  4.24264069])
    >>> LA.norm(c, ord=1, axis=1)
    array([ 6.,  6.])
    
    Using the `axis` argument to compute matrix norms:
    
    >>> m = np.arange(8).reshape(2,2,2)
    >>> LA.norm(m, axis=(1,2))
    array([  3.74165739,  11.22497216])
    >>> LA.norm(m[0, :, :]), LA.norm(m[1, :, :])
    (3.7416573867739413, 11.224972160321824)

Questions

Consider the two vectors:

r1 = np.array([1.0, -2.0, 3.0])
r2 = np.array([2.0, 2.0, -3.0])

Use numpy functionality to:

  • Find a unit vector in the direction of r1:

# Your code here
  • Find the angle between r1 and r2 in radians:

# Your code here
  • The projection of r2 on r1:

# Your code here
  • Use numpy.cross to find the cross product between r2 and r1. Then, verify numerically (using numpy functionality) that r1(r1×r2)=0\vec{r}_1\cdot\left(\vec{r}_1\times\vec{r}_2\right) = 0.

# Your code here

Matrices

Now, let’s go to two dimensional arrays. Of course, you can think of the vectors as two dimensional arrays. In particular, you can think of them as a matrix with a single column. To do this in numpy you need to reshape the vectors:

# The original vector
a
array([ 1.2, 2. , 3. , -1. , 2. ])
# The vector as a single column matrix
acm = a.reshape((5, 1))
acm
array([[ 1.2], [ 2. ], [ 3. ], [-1. ], [ 2. ]])
# There is one more way you can achieve the same thing
# (and I am going to be using it because it requires less typing)
acm = a[:, None]
acm
array([[ 1.2], [ 2. ], [ 3. ], [-1. ], [ 2. ]])

The code a[:, None] means “pretend that the array has one more index that is always zero.”

Now, since acm is a 2D array:

acm.ndim
2

Look also at the shape:

acm.shape
(5, 1)

So, the shape is (5, 1). This means that there are 5 rows and one colum. You can access elements in this 2D array with indices like this:

acm[0, 0]
1.2
acm[1, 0]
2.0

And so on. Of course, the second index can only be zero because you only have one column.

Let’s now make a 2D matrix to play with:

A = np.array([[2.0, -1.0, 0.0, 0.0, 0.0],
              [-1.0, 2.0, -1.0, 0.0, 0.0],
              [0.0, -1.0, 2.0, -1.0, 0.0],
              [0.0, 0.0, -1.0, 2.0, -1.0],
              [0.0, 0.0, 0.0, -1.0, 2.0]
             ])
A
array([[ 2., -1., 0., 0., 0.], [-1., 2., -1., 0., 0.], [ 0., -1., 2., -1., 0.], [ 0., 0., -1., 2., -1.], [ 0., 0., 0., -1., 2.]])

Notice that we used an list of rows and each row was a list of numbers. Here are the dimensions of A:

A.ndim
2
A.shape
(5, 5)

Let’s access some elements of A:

A[0, 0]
2.0
A[2, 2]
2.0

You can access the 3rd row like this:

A[2]
array([ 0., -1., 2., -1., 0.])

You can get a submatrix, say the first 3x3 submatrix like this:

A[:3, :3]
array([[ 2., -1., 0.], [-1., 2., -1.], [ 0., -1., 2.]])

And so on.

You can multiply matrices with numbers:

2.0 * A
array([[ 4., -2., 0., 0., 0.], [-2., 4., -2., 0., 0.], [ 0., -2., 4., -2., 0.], [ 0., 0., -2., 4., -2.], [ 0., 0., 0., -2., 4.]])
A * 3.0
array([[ 6., -3., 0., 0., 0.], [-3., 6., -3., 0., 0.], [ 0., -3., 6., -3., 0.], [ 0., 0., -3., 6., -3.], [ 0., 0., 0., -3., 6.]])

You can raise them to a power:

A ** 3
array([[ 8., -1., 0., 0., 0.], [-1., 8., -1., 0., 0.], [ 0., -1., 8., -1., 0.], [ 0., 0., -1., 8., -1.], [ 0., 0., 0., -1., 8.]])

You can pass them through functions:

np.exp(A)
array([[7.3890561 , 0.36787944, 1. , 1. , 1. ], [0.36787944, 7.3890561 , 0.36787944, 1. , 1. ], [1. , 0.36787944, 7.3890561 , 0.36787944, 1. ], [1. , 1. , 0.36787944, 7.3890561 , 0.36787944], [1. , 1. , 1. , 0.36787944, 7.3890561 ]])

You can add them together:

A + A
array([[ 4., -2., 0., 0., 0.], [-2., 4., -2., 0., 0.], [ 0., -2., 4., -2., 0.], [ 0., 0., -2., 4., -2.], [ 0., 0., 0., -2., 4.]])

You can multiply them together (assuming that the dimensions match):

A.dot(A)
array([[ 5., -4., 1., 0., 0.], [-4., 6., -4., 1., 0.], [ 1., -4., 6., -4., 1.], [ 0., 1., -4., 6., -4.], [ 0., 0., 1., -4., 5.]])

You can multiply a matrix with a vector:

A.dot(a)
array([ 0.4, -0.2, 5. , -7. , 5. ])

Remember the column matrix acm:

acm
array([[ 1.2], [ 2. ], [ 3. ], [-1. ], [ 2. ]])

You can also multiply the matrix with it since the dimensions match:

A.dot(acm)
array([[ 0.4], [-0.2], [ 5. ], [-7. ], [ 5. ]])

You can also multiply acm (5x1) with its transpose (1x5). The transpose is:

acm.T
array([[ 1.2, 2. , 3. , -1. , 2. ]])

And here is the result of the multiplication:

acm.dot(acm.T)
array([[ 1.44, 2.4 , 3.6 , -1.2 , 2.4 ], [ 2.4 , 4. , 6. , -2. , 4. ], [ 3.6 , 6. , 9. , -3. , 6. ], [-1.2 , -2. , -3. , 1. , -2. ], [ 2.4 , 4. , 6. , -2. , 4. ]])

Here is another example of the transpose operator:

B = np.arange(25).reshape(5, 5)
B
array([[ 0, 1, 2, 3, 4], [ 5, 6, 7, 8, 9], [10, 11, 12, 13, 14], [15, 16, 17, 18, 19], [20, 21, 22, 23, 24]])
B.T
array([[ 0, 5, 10, 15, 20], [ 1, 6, 11, 16, 21], [ 2, 7, 12, 17, 22], [ 3, 8, 13, 18, 23], [ 4, 9, 14, 19, 24]])

There are some special numpy functions for creating arrays. Here is how to make an array of zeros:

np.zeros((4, 5))
array([[0., 0., 0., 0., 0.], [0., 0., 0., 0., 0.], [0., 0., 0., 0., 0.], [0., 0., 0., 0., 0.]])

An array of ones:

np.ones((5, 2))
array([[1., 1.], [1., 1.], [1., 1.], [1., 1.], [1., 1.]])

A unit matrix:

np.eye(5)
array([[1., 0., 0., 0., 0.], [0., 1., 0., 0., 0.], [0., 0., 1., 0., 0.], [0., 0., 0., 1., 0.], [0., 0., 0., 0., 1.]])

Now, we typically want to do some more linear algebra with matrices. For example, we may want to get the determinant of a matrix:

np.linalg.det(A)
6.0

Or we may want to solve a linear system of the form:

Ax=b,A x = b,

where bb is a given vector and xx is the uknown vector we are looking for. Let’s create and solve such a linear system.

# I already have 5x5 matrix A, just define the b:
b = np.array([2.0, -1.0, 3.0, 4.0, 2.0])
# Here is how to solve the linear system:
x = np.linalg.solve(A, b)
x
array([4.16666667, 6.33333333, 9.5 , 9.66666667, 5.83333333])

Let’s see if the solution works as expected by calculating AxAx:

print('Ax = ', A.dot(x))
print('b  = ', b)

You can also find the eigenvalues of a matrix:

lam, v = np.linalg.eig(A)

The eigenvalues are:

lam

The eigenvectors are put in a matrix:

v

So, each column of that matrix is an eigenvector. Let’s test that

Avi=λivi,A v_i = \lambda_i v_i,

for the ii-th eigenvalue:

# Which eigenvalue do you want to test:
i = 0
print('Av_i        = ', A.dot(v[:, i]))
print('lamda_i v_i = ', lam[i] * v[:, i])

By the way, v[:, i] gives you the ii-th column of the matrix v.

Questions

  • Consider the quadratic function:

    f(x1,x2)=x12x1x2+x22+3x1+4x2.f(x_1, x_2) = x_1^2 -x_1x_2 + x_2^2 + 3x_1 + 4x_2.

    We are interested in finding its minimum. A necessary condition for the minimum is that the gradient of f(x1,x2)f(x_1, x_2) is zero. The gradient is:

    f(x1,x2)x1=2x1x2+3=0,\frac{\partial f(x_1,x_2)}{\partial x_1} = 2x_1 -x_2 + 3 = 0,

    and

    f(x1,x2)x2=x1+2x2+4=0.\frac{\partial f(x_1, x_2)}{\partial x_2} = -x_1 + 2x_2 + 4 = 0.

    Rearranging, we get the two equations:

    2x1x2=3,x1+2x2=4.\begin{array}{ccc} 2x_1 - x_2 &=& -3,\\ -x_1 + 2x_2 &=& -4.\\ \end{array}

    These can also be written in matrix-vector form: $$

[2112]\begin{bmatrix} 2 &-1\\ -1 & 2 \end{bmatrix}
[x1x2]\begin{bmatrix} x_1\\ x_2 \end{bmatrix}

= \begin{bmatrix} -3\ -4 \end{bmatrix}. $Usethenumpyfunctionalitytosolvethislinearsystemandthusfindtheminimumof Use the numpy functionality to solve this linear system and thus find the minimum of f(x_1,x_2)$.

# your code here
  • For the point you found above to be a true minimum of the f(x1,x2)f(x_1,x_2) (and not a maximum for example) is that the matrix of second derivatives (the so-called Hessian matrix) must be positive definite, i.e., that the matrix has only positive eigenvalues. The hessian matrix in this problem is: $$ H = \begin{bmatrix} \frac{\partial^2 f(x_1,x_2)}{\partial x_1^2}& \frac{\partial^2 f(x_1,x_2)}{\partial x_1\partial x_2}\ \frac{\partial^2 f(x_1,x_2)}{\partial x_1\partial x_2} & \frac{\partial^2 f(x_1,x_2)}{\partial x_1^2} \end{bmatrix} =

[2112]\begin{bmatrix} 2 &-1\\ -1 & 2 \end{bmatrix}

$Usenumpyfunctionalitytoverifythatthematrix Use numpy functionality to verify that the matrix H$ is indeed positive definite.

# your code here
  • If HH is a positive definite matrix, then the following property is true for any vector x0x\not = 0:

    xTHx0.x^THx \ge 0.

    In words, if you sandwich HH around the same vector produce a scalar, that scalar is a positive number. Test this hypothesis numerically using numpy functioanlity. Hint: Use the dot function and twich to calculate xTHxx^THx.

# here is a random vector for you to play with:
x = np.random.randn(2)
print('x = ', x)
# Your code here: Compute x^T * H * x:

Multi-dimensional numpy arrays

You can have 3D numpy arrays:

C = np.random.randn(10, 3, 4)
C
C.ndim
C.shape
C[0, 2, 3]
C[0]
C[0, 1]

Or even higher dimensional:

D = np.random.randn(3, 4, 5, 2, 5)
D.ndim
D.shape

Of course, most of the time we will be working with 1D (Vectors) and 2D (Matrices).