 # Intro to Linear Algebra Using NumPy

Last updated: November 19th, 2019  # Intro to Linear Algebra using NumPy¶

Numpy arrays, when created with the correct set of dimensions, will resemble algebraic structures like vectors or matrices. As expected, NumPy will be able to perform the expected primitive operations associated with them. It has a very comprehensive linear algebra sub-system. ## Hands on!¶

In :
import sys
import numpy as np

$$\large A = \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \end{bmatrix}$$
In :
A = np.array([[1, 2, 3], [4, 5, 6]])
A

Out:
array([[1, 2, 3],
[4, 5, 6]])
$$\large B = \begin{bmatrix} 5 & 5 & 5 \\ 5 & 5 & 5 \end{bmatrix}$$
In :
B = np.full((2, 3), 5)
B

Out:
array([[5, 5, 5],
[5, 5, 5]])
$$\normalsize v = \begin{bmatrix} 6 \\ 7 \\ 8 \end{bmatrix}$$
In :
v = np.array([6, 7, 8])


### Operations with Matrices and Vectors¶

Sum of matrices:

\begin{equation*} \large A + \large B = \begin{bmatrix} 1 & 2 & 3\\ 4 & 5 & 6 \end{bmatrix} + \begin{bmatrix} 5 & 5 & 5 \\ 5 & 5 & 5 \end{bmatrix} = \begin{bmatrix} 6 & 7 & 8 \\ 9 & 10 & 11 \end{bmatrix} \end{equation*}
In :
A + B

Out:
array([[ 6,  7,  8],
[ 9, 10, 11]])

Subtraction of matrices:

\begin{equation*} \large A - \large B = \begin{bmatrix} 1 & 2 & 3\\ 4 & 5 & 6 \end{bmatrix} - \begin{bmatrix} 5 & 5 & 5 \\ 5 & 5 & 5 \end{bmatrix} = \begin{bmatrix} -4 & -3 & -2 \\ -1 & 0 & 1 \end{bmatrix} \end{equation*}
In :
A - B

Out:
array([[-4, -3, -2],
[-1,  0,  1]])

(Not to be confused with matrix multiplication or dot product):

\begin{equation*} \large A \circ \large B = \begin{bmatrix} 1 & 2 & 3\\ 4 & 5 & 6 \end{bmatrix} \circ \begin{bmatrix} 5 & 5 & 5 \\ 5 & 5 & 5 \end{bmatrix} = \begin{bmatrix} 5 & 10 & 15 \\ 20 & 25 & 30 \end{bmatrix} \end{equation*}
In :
A * B

Out:
array([[ 5, 10, 15],
[20, 25, 30]])

Division also works "per-entry":

In :
A / B

Out:
array([[0.2, 0.4, 0.6],
[0.8, 1. , 1.2]])

Matrix multiplication:

\begin{equation*} \large A \cdot \normalsize v = \begin{bmatrix} 1 & 2 & 3\\ 4 & 5 & 6 \end{bmatrix} \cdot \begin{bmatrix} 6 \\ 7 \\ 8 \\ \end{bmatrix} = \begin{bmatrix} 44 \\ 107 \end{bmatrix} \end{equation*}
In :
np.dot(A, v)

Out:
array([ 44, 107])

In Python 3.5, the @ operator was introduced for matrix multiplication:

In :
A @ v

Out:
array([ 44, 107])

Matrix transpose:

In :
A.T

Out:
array([[1, 4],
[2, 5],
[3, 6]])

### The identity matrix¶

An identity matrix is a square matrix, of size n, with ones on the main diagonal and zeros elsewhere:

\begin{equation*} \normalsize I_{n} = \begin{bmatrix} 1&0&0&\cdots &0\\ 0&1&0&\cdots &0\\ 0&0&1&\cdots &0\\ \vdots &\vdots &\vdots &\ddots &\vdots \\ 0&0&0&\cdots &1 \end{bmatrix} \end{equation*}

So, for example, an identity matrix of size 3:

\begin{equation*} \normalsize I_{3} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \end{equation*}
In :
I = np.identity(3)
I

Out:
array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])
\begin{equation*} \normalsize C = \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{bmatrix} \end{equation*}
In :
np.arange(1, 10)

Out:
array([1, 2, 3, 4, 5, 6, 7, 8, 9])
In :
C = np.arange(1, 10).reshape((3, 3))
C

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

Properties of the identity matrix:

\begin{equation*} \large A_{n} \cdot I_{n} = A_{n} \end{equation*}\begin{equation*} \large I_{n} \cdot A_{n} = A_{n} \end{equation*}
In :
I @ C

Out:
array([[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]])
In :
C @ I

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

### Inverse of a matrix:¶

The formula to find the inverse of a 2x2 matrix is:

\begin{equation*} \normalsize A ^ {-1} = \begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix}^{-1} = \frac{1}{\lvert A\rvert} \begin{bmatrix} a_{22} & -a_{12} \\ -a_{21} & a_{11}\end{bmatrix} \end{equation*}

|A| is the determinant of a matrix, which can be computed with np.linalg.det. So, as an example, given the following matrix:

\begin{equation*} \normalsize D = \begin{bmatrix} 1 & 2 \\ 3 & 4 \\ \end{bmatrix} \end{equation*}
In :
D = np.array([
[1, 2],
[3, 4]
])

D

Out:
array([[1, 2],
[3, 4]])

The determinant can be calculated using:

In :
np.linalg.det(D)

Out:
-2.0000000000000004

And the inverse with np.linalg.inv:

In :
np.linalg.inv(D)

Out:
array([[-2. ,  1. ],
[ 1.5, -0.5]])

A square matrix that is not invertible is called singular or degenerate. A square matrix is singular if and only if its determinant is 0. Singular matrices are rare in the sense that a square matrix randomly selected from a continuous uniform distribution on its entries will almost never be singular.

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

np.linalg.det(singular)

Out:
0.0

Calculating its inverse matrix will raise an error:

In :
np.linalg.inv(singular)

---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
<ipython-input-21-2285bd57e24f> in <module>
----> 1 np.linalg.inv(singular)

/usr/local/lib/python3.6/site-packages/numpy/linalg/linalg.py in inv(a)
526     signature = 'D->D' if isComplexType(t) else 'd->d'
527     extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> 528     ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
529     return wrap(ainv.astype(result_t, copy=False))
530

/usr/local/lib/python3.6/site-packages/numpy/linalg/linalg.py in _raise_linalgerror_singular(err, flag)
87
88 def _raise_linalgerror_singular(err, flag):
---> 89     raise LinAlgError("Singular matrix")
90
91 def _raise_linalgerror_nonposdef(err, flag):

LinAlgError: Singular matrix

The main characteristic of the inverse matrix is that, it acts as a reciprocal of the matrix:

\begin{equation*} \large A \cdot A^{-1} = I \end{equation*}\begin{equation*} \large A^{-1} \cdot A = I \end{equation*}
In :
D @ np.linalg.inv(D)

Out:
array([[1.0000000e+00, 0.0000000e+00],
[8.8817842e-16, 1.0000000e+00]])
In :
np.round(D @ np.linalg.inv(D))

Out:
array([[1., 0.],
[0., 1.]]) ### Useful Numpy functions¶

We'll now take the chance to see in more detail the functions we've been using, the most common and useful ones in the numpy world:

arange

Return evenly spaced values within a given interval.

In :
range(5)

Out:
range(0, 5)
In :
np.arange(10)

Out:
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In :
np.arange(5, 10)

Out:
array([5, 6, 7, 8, 9])
In :
np.arange(0, 1, .1)

Out:
array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])

linspace

Return evenly spaced numbers over a specified interval.

In :
np.linspace(0, 1, 5)

Out:
array([0.  , 0.25, 0.5 , 0.75, 1.  ])
In :
np.linspace(0, 1, 20)

Out:
array([0.        , 0.05263158, 0.10526316, 0.15789474, 0.21052632,
0.26315789, 0.31578947, 0.36842105, 0.42105263, 0.47368421,
0.52631579, 0.57894737, 0.63157895, 0.68421053, 0.73684211,
0.78947368, 0.84210526, 0.89473684, 0.94736842, 1.        ])
In :
np.linspace(0, 1, 20, False)

Out:
array([0.  , 0.05, 0.1 , 0.15, 0.2 , 0.25, 0.3 , 0.35, 0.4 , 0.45, 0.5 ,
0.55, 0.6 , 0.65, 0.7 , 0.75, 0.8 , 0.85, 0.9 , 0.95])

zeros, ones, empty

In :
np.zeros(5)

Out:
array([0., 0., 0., 0., 0.])
In :
np.zeros((3, 3))

Out:
array([[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.]])
In :
np.zeros((3, 3), dtype=np.int)

Out:
array([[0, 0, 0],
[0, 0, 0],
[0, 0, 0]])
In :
np.ones(5)

Out:
array([1., 1., 1., 1., 1.])
In :
np.ones((3, 3))

Out:
array([[1., 1., 1.],
[1., 1., 1.],
[1., 1., 1.]])
In :
np.empty(5)

Out:
array([1., 1., 1., 1., 1.])
In :
np.empty((2, 2))

Out:
array([[0.25, 0.5 ],
[0.75, 1.  ]])

Full filled

In :
np.full((2,2), 7)

Out:
array([[7, 7],
[7, 7]])
In :
np.full((3, 3), 1)

Out:
array([[1, 1, 1],
[1, 1, 1],
[1, 1, 1]])

Random filled

In :
np.random.random((2,2))

Out:
array([[0.95994826, 0.11944492],
[0.67259477, 0.26004165]])

Identity and eye

In :
np.identity(3)

Out:
array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])
In :
np.eye(3, 3)

Out:
array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])
In :
np.eye(8, 4)

Out:
array([[1., 0., 0., 0.],
[0., 1., 0., 0.],
[0., 0., 1., 0.],
[0., 0., 0., 1.],
[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.]])
In :
np.eye(8, 4, k=1)

Out:
array([[0., 1., 0., 0.],
[0., 0., 1., 0.],
[0., 0., 0., 1.],
[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.]])
In :
np.eye(8, 4, k=-3)

Out:
array([[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 0.],
[1., 0., 0., 0.],
[0., 1., 0., 0.],
[0., 0., 1., 0.],
[0., 0., 0., 1.],
[0., 0., 0., 0.]])
In :
np.eye(8, 4, k=-3, dtype=np.int)

Out:
array([[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0],
[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1],
[0, 0, 0, 0]]) 