3. Numpy.ipynb

Kernel: Python 3 (system-wide)

Table of Contents

• 1  In this notebook
• 2  More on arrays
• 3  Broadcasting
• 4  Some canonical matrices/vectors
• 5  Inverting a matrix
• 6  Solving linear systems
• 7  Exercises

(More on) Numpy

Camilo A. Garcia Trillos - 2020

---

In this notebook

We learn more about the Numpy package: we look at arrays and matrices, and some of the associated modules.

We start by importing the modules we will use in this notebook: only numpy.

import numpy as np

More on arrays

Let us look deeper at the numpy library. The numpy library contains useful tools to work with vectors and matrices, solve linear operations and generate random entries.

We have a look at data types and their implications, basic operations, some canonical matrices, matrix inversion and solution of linear systems.

a = np.array([1,3,6,1])

a

array([1, 3, 6, 1])

Two important methods for arrays are size and shape. The former give sus the total number of elements in an array. The latter tells us how these elements are organised, whcih is particularly useful to define matrices or higher order tensors.

a.size # The number of total elements in the array

4

a.shape

(4,)

This means that at this time, all elements are organised in a unique row. We can render this vector into a 'column vector' or into a 2x2 matrix with the same entries simply assigning the shape we want or by using the reshape method.

print('original')
print(a)
a.shape=(1,4)
print('column vector (inplace)')
print(a)
b = a.reshape(2,2)
print('matrix')
print(b)

original
[1 3 6 1]
column vector (inplace)
[[1 3 6 1]]
matrix
[[1 3]
 [6 1]]

When assigning the shape property, the array itself changes (this is called 'inplace'):

a

array([[1, 3, 6, 1]])

print(a[0,1]) # This is OK
print(a[0]) # this is OK
a[1] #This is an error

3
[1 3 6 1]

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-8-685c588a6827> in <module>
      1 print(a[0,1]) # This is OK
      2 print(a[0]) # this is OK
----> 3 a[1] #This is an error

IndexError: index 1 is out of bounds for axis 0 with size 1

The above error occurs because there is only one entry on the first dimension of the array.

Arrays, unlike lists, cannot in general combine different types of entries. The vector we created above was an integer vector (note the lack of dots at the end of all prints!). Hence, if we try to assign a float, we will get a roundup version of it

a[0,1] = 9.2

a

array([[1, 9, 6, 1]])

a_1 = np.array([1,2,3,4],dtype=float)
a_1

array([1., 2., 3., 4.])

As mentioned above, 9.2 becomes 9. Let su try now by creating a vector of floats from the start:

b = np.array([1.0,3,6])

b

array([1., 3., 6.])

b[1] = 9.2

b

array([1. , 9.2, 6. ])

Note the difference! The array b was initialized to have floats, thanks to writing 1.0 instead of 1. Thus, in contrast to a, it contains 9.2 (instead of 9, as the array a does)

As lists, arrays are assigned by reference. What follows should not surprise you if you followed the first notebook.

c = b
c
c[0] = 7.1
print('c:',c)
print('b:', b)

c: [7.1 9.2 6. ]
b: [7.1 9.2 6. ]

If an independent copy of a given array is needed, we can use the method copy

c=b.copy()
c
c[1]=1
print(c)
print(b)

[7.1 1.  6. ]
[7.1 9.2 6. ]

In the previous notebook, we looked at the method arange. Here are ways to use it:

d = np.arange(3, 42, 7)
d

array([ 3, 10, 17, 24, 31, 38])

d = np.arange(38, 2, -7)
d

array([38, 31, 24, 17, 10,  3])

d = np.arange(10)
d

array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

Broadcasting

We mentioned in the previous notebook that there was an exception to the rule of operating with the exact same dimensions. This exception is called broadcasting. Let us first present an example: try to figure out what is happening

Example 1

d.shape = (5,2)
d

array([[0, 1],
       [2, 3],
       [4, 5],
       [6, 7],
       [8, 9]])

e = np.array([1,3])

print(d * e)
print(d.shape, e.shape)

[[ 0  3]
 [ 2  9]
 [ 4 15]
 [ 6 21]
 [ 8 27]]
(5, 2) (2,)

e_0 = np.arange(5)
e_0.T

array([0, 1, 2, 3, 4])

Example 2

e.shape = (1,2)
print(e * d)
print(d.shape, e.shape)

[[ 0  3]
 [ 2  9]
 [ 4 15]
 [ 6 21]
 [ 8 27]]
(5, 2) (1, 2)

Example 3

f = np.array([1,1,1,2,2]).reshape(-1,1)
print('f: ',f, 'and the shape of f:', f.shape)

f:  [[1]
 [1]
 [1]
 [2]
 [2]] and the shape of f: (5, 1)

print(f*d)
print(f.shape, d.shape)

[[ 0  1]
 [ 2  3]
 [ 4  5]
 [12 14]
 [16 18]]
(5, 1) (5, 2)

Example 4

However, recall that the following gives an error:

g = np.array([1,2,3])
print(d.shape,g.shape)
print(g*d) # This gives an error

(5, 2) (3,)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-39-276b4363f031> in <module>
      1 g = np.array([1,2,3])
      2 print(d.shape,g.shape)
----> 3 print(g*d) # This gives an error

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

If you have not figure it out, the explanation is as follows: suppose we are operating the arrays a and b. Suppose that the shape of a is $(m_1,m_2)$, while the shapes of b are $(n_1,n_2)$:

• If $n_i \neq m_i$ and both are different from one or empty for some $i=1, 2$, then the operation cannot be done.

• If $n_i \neq m_i$ and one of the two is equal to one or empty for some $i=1, 2$, then the operation is done by completing the missing dimensions via copying the remaining components.

Look again to the example to make sure you understand this process. In particular, see how the empty spaces can be accomodated to the above template.

Note: arrays can be tensors, i.e. have more than 2 dimensions. Broadcasting is generalised as you would expect.

Example 5

Observe the following example

d[0,0] = 1

d

array([[1, 1],
       [2, 3],
       [4, 5],
       [6, 7],
       [8, 9]])

1/d

array([[1.        , 1.        ],
       [0.5       , 0.33333333],
       [0.25      , 0.2       ],
       [0.16666667, 0.14285714],
       [0.125     , 0.11111111]])

The above works, because 1 is taken as an array of empty size, so it can be broadcasted by repeating itself to operate on each entry.

Note though that the following is an error:

d ** (-1) # This raises an error

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-43-d2424fe8248b> in <module>
----> 1 d ** (-1) # This raises an error

ValueError: Integers to negative integer powers are not allowed.

What's wrong? We initialized d with integers, and there is no well defined change of type for this operation.

d.dtype # Observe the type of this array

dtype('int64')

d = d.astype(float) # Set it to float

d.dtype # Observe again the type

dtype('float64')

d ** (-1) # The correct answer appears.

array([[1.        , 1.        ],
       [0.5       , 0.33333333],
       [0.25      , 0.2       ],
       [0.16666667, 0.14285714],
       [0.125     , 0.11111111]])

(d ** (-1)) == (1/d)

array([[ True,  True],
       [ True,  True],
       [ True,  True],
       [ True,  True],
       [ True,  True]])

Some canonical matrices/vectors

Let us introduce some canonical matrices and vectors.

# Zeros and ones matrix/vector

print (np.zeros(10))
a = np.ones((2,5))
print(a)

[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[[1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]]

By default these are float. If an integer is desired, this has to be stated explicitly.

b=np.zeros(10, dtype=int)
print (b)

[0 0 0 0 0 0 0 0 0 0]

Some other useful functions are zeros_like and ones_like: they create zeros and ones of the lentgh and type of the parameter

print(np.zeros_like(a))
print(np.ones_like(b))

[[0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]]
[1 1 1 1 1 1 1 1 1 1]

We have also an identity matrix and the function diag. Let us look at what they do:

c=np.eye(5)

c[2][2]=10
print(c)
print('=====')
print(np.diag(c))

[[ 1.  0.  0.  0.  0.]
 [ 0.  1.  0.  0.  0.]
 [ 0.  0. 10.  0.  0.]
 [ 0.  0.  0.  1.  0.]
 [ 0.  0.  0.  0.  1.]]
=====
[ 1.  1. 10.  1.  1.]

Diag can also be used to build a diagonal matrix, as shown below.

np.diag(np.arange(8))

array([[0, 0, 0, 0, 0, 0, 0, 0],
       [0, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 2, 0, 0, 0, 0, 0],
       [0, 0, 0, 3, 0, 0, 0, 0],
       [0, 0, 0, 0, 4, 0, 0, 0],
       [0, 0, 0, 0, 0, 5, 0, 0],
       [0, 0, 0, 0, 0, 0, 6, 0],
       [0, 0, 0, 0, 0, 0, 0, 7]])

Inverting a matrix

Recall how to generate a random Gaussian matrix: here we generate a matrix of size 3x3, where each entry is i.i.d. with distribution $\mathcal N(0,1)$

rng = np.random.default_rng()
g = rng.normal(0,1,(3,3))
g

array([[-0.03925264,  1.26893031,  0.72882569],
       [-1.55485771,  0.99537217,  0.41630284],
       [ 0.23072996, -1.35665215, -2.23083365]])

We can check that the matrix is invertible, by looking at the rank of the matrix

np.linalg.matrix_rank(g)

3

If the rank is 3 then we can invert the matrix g. We already showed that g ** (-1) does not give the matrix inverse. Instead it computes the reciprocals of each component:

g ** (-1)

array([[-25.47599223,   0.78806534,   1.37207019],
       [ -0.64314567,   1.00464935,   2.40209746],
       [  4.33407094,  -0.73710862,  -0.44826292]])

Instead the package numpy has a specific command to obtain the matrix inverse:

g_inv = np.linalg.inv(g)

g_inv

array([[ 0.58206759, -0.64755304,  0.0693228 ],
       [ 1.185619  ,  0.02833322,  0.39263563],
       [-0.6608167 , -0.08420539, -0.67986922]])

Let's check:

g@g_inv

array([[ 1.00000000e+00, -5.59342334e-18, -4.63920999e-17],
       [ 3.17230549e-17,  1.00000000e+00,  2.69221949e-17],
       [-3.53374221e-16, -2.41117475e-17,  1.00000000e+00]])

Apart from rounding errors this is basically the identity matrix.

We need to be careful though:

• Inverting matrices is expensive

• Inversion is done numerically, and this algorithm is not stable when applied to certain matrices.

Let us illustrate this point with an example. We are going to use another canonical matrix, the Hilbert matrix. This matrix is famous for being ill conditioned (that is, there is a huge gap between its eigenvalues), which makes it a bad example for inversion algorithms.

from scipy.linalg import hilbert

dim = 13
ill_conditioned = hilbert(dim)
#ill_conditioned

np.linalg.cond(ill_conditioned)

3.344143497338461e+18

ill_inverse = np.linalg.inv(ill_conditioned)
np.diag(ill_conditioned @ ill_inverse )   # If inverse is perfect, all entries here would be one (can you see why?)

array([1.        , 0.99999934, 0.99995757, 1.00088501, 1.01743457,
       1.15358779, 0.99473281, 1.15498136, 0.526914  , 0.66719027,
       0.78461637, 1.04838723, 0.98567695])

b = np.ones(dim)
ill_conditioned @ (ill_inverse @ b) # If inverse is perfect, all entries here would be one (can you see why?)

array([-0.98183828,  0.23957478,  0.45739563,  0.5548842 ,  0.61760592,
        0.66364201,  0.69945406,  0.72825585,  0.75195987,  0.77181838,
        0.7886987 ,  0.80322435,  0.81585616])

Note that the above tests show that we are not really getting a good approximation of the inverse.

There are several ways to account for this problem. One can add a small amount to the diagonal

ill_min = np.min(ill_conditioned)
ill_max = np.max(ill_conditioned)
print(ill_min,ill_max)

0.04 1.0

better_cond = ill_conditioned + np.diag (ill_min*0.05*np.ones(dim) )
better_inv = np.linalg.inv(better_cond)

ill_conditioned @ (better_inv @ b) # If inverse is perfect, all entries here would be one (can you see why?)

array([1.00160478, 0.99015151, 0.99379173, 1.00807935, 1.01921242,
       1.02373631, 1.02185175, 1.01473545, 1.00365078, 0.98968786,
       0.97371675, 0.95640769, 0.93826756])

Note that this solution is much closer to the one vector, and so the inverse is close to what we would like.

But this solution has a bias, that is, we have a systematic error. A more expensive alternative is to add an error that is neglected in probability: we add a small random amount to the matrix, and then average out a couple of solutions

noise_factor = np.diag(rng.random(dim) *0.02)

np.min(noise_factor), np.mean(noise_factor), np.max(noise_factor)

(0.0, 0.000710371431005939, 0.018187671656453738)

better_cond2 = ill_conditioned + noise_factor
better_inv2 = np.linalg.inv(better_cond2)

ill_conditioned @ (better_inv2 @ b)

array([1.00673411, 0.99159814, 1.00147924, 1.01555965, 1.02350878,
       1.02394675, 1.01798683, 1.00716335, 0.99285003, 0.97614086,
       0.95786698, 0.93864577, 0.91892998])

This works specially well when the matrix is a calculated covariance matrix: in this case one can randomly perturb the data before calculating the covariance.

Solving linear systems

Sometimes we are not interested in finding the inverse matrix but rather in solving a liner system like

$$A {\bf x} ={\bf b}$$

Doing this directly is usually more efficient.

There is a whole theory on how to solve linear equations(this is at the core of numerical analysis), but we will simply use the predifined functions.

x = np.linalg.solve(ill_conditioned, b)
x

array([ 1.92939810e+00, -5.07018480e+02,  2.71906204e+04, -5.77279548e+05,
        6.28427763e+06, -3.99817112e+07,  1.59846383e+08, -4.17481359e+08,
        7.22593695e+08, -8.21390449e+08,  5.89180874e+08, -2.41720847e+08,
        4.32198922e+07])

We test it:

ill_conditioned @ x

array([0.99999998, 1.00000003, 1.00000001, 1.        , 1.        ,
       1.        , 1.00000001, 1.00000001, 0.99999999, 1.        ,
       1.        , 1.        , 0.99999998])

Very close to $b$ (a vector of ones). Compare with the result we obain using the inverse:

x2 = ill_inverse@b
ill_conditioned@x2

array([-0.98183828,  0.23957478,  0.45739563,  0.5548842 ,  0.61760592,
        0.66364201,  0.69945406,  0.72825585,  0.75195987,  0.77181838,
        0.7886987 ,  0.80322435,  0.81585616])

And compare with what we obtain using the modified inverse (with the small perturbation):

x3 = better_inv@b
ill_conditioned@x3

array([1.00160478, 0.99015151, 0.99379173, 1.00807935, 1.01921242,
       1.02373631, 1.02185175, 1.01473545, 1.00365078, 0.98968786,
       0.97371675, 0.95640769, 0.93826756])

Better... but solving the linear system is superior.

Remark: Finding an inverse might be advantageous if several linear systems with the smae matrices have to be solved one after the other.

Exercises

1. Create a function 'special_gaussian(n)' that generates a sample of $n$ independent Gaussians, where the variance of the $i$th component is $1/i$. (Note: You must use np.random.randn). Do not forget to avoid erros and test your function.

def special_gaussian(n):
    rng = np.random.default_rng()
    aux = np.arange(n,dtype=float)
    for i in range(1,n+1):
        aux[i-1] = rng.normal(0,np.sqrt(1/i),1)
    return aux

special_gaussian(10)

array([-1.40600606, -0.3442405 ,  0.29321358,  0.1680397 ,  0.12872512,
       -0.33458345,  0.53077877,  0.12625882, -0.41157362, -0.37847876])

2. We consider a one-period market model in a finite probability space.

Create a function 'complete_market_1p(S0,S1)' that receives a vector S0 with the initial prices of the assets in a market, and a matrix S1 that gives the prices of each asset in the market under each scenario, and returns a boolean that states if the market is complete or not.

import numpy as np
def complete_market_1p(S_0,S_1):
    disc_S_0 = S_0.reshape(-1,1)
    risk_free_asset = S_1[0,:]
    disc_S_1 = (1/risk_free_asset)*S_1
    m = np.hstack((disc_S_1,disc_S_0))
    r_1 = np.linalg.matrix_rank(m)
    r_2 = np.linalg.matrix_rank(disc_S_1)
    Q = np.linalg.pinv(disc_S_1)@S0
    if (r_1 == r_2)&(Q[(Q>0)&(Q<1)]==Q).all()&(r_1==min(S_1.shape[0],S_1.shape[1])):
        return True
    else:
        return False

S0=np.array([1,3,4])
S1=np.array([
    [1,2,4,4],
    [1,3,6,5],
    [1,3,7,5]
])
print(complete_market_1p(S0,S1))

False