Introduction to linear algebra

You can always make arrays by typing them in. There are many convenient ways to make special ones though. For example, you can make an array of all ones or zeros with these:

np.zeros(shape=[3, 3])

array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.]])

np.ones(shape=[3, 3])

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

You can make an identity matrix with:

np.eye(N=3)

array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

or a diagonal array:

np.diag([1, 2, 3])

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

If you need a lower triangular array:

np.tri(3)

array([[1., 0., 0.],
       [1., 1., 0.],
       [1., 1., 1.]])

It takes some getting use to how to use arrays with algebra.

To do matrix multiplication you use the @ operator (This is new in Python 3.5), or the numpy.dot function. If you are not familiar with the idea of matrix multiplication you should review it at https://en.wikipedia.org/wiki/Matrix_multiplication.

We write matrix multiplication as: \(\mathbf{A} \mathbf{B}\). We cannot multiply any two arrays; their shapes must follow some rules. We can multiply any two arrays with these shapes:

(m, c) * (c, n) = (m, n)

In other words the number of columns in the first array must equal the number of rows in the second array. This means it is not generally true that \(\mathbf{A} \mathbf{B} = \mathbf{B} \mathbf{A}\).

A @ B

array([[3., 3.],
       [7., 7.]])

This is the older way to do matrix multiplication.

np.dot(A, B)

array([[3., 3.],
       [7., 7.]])

These rules are true:

1. \((k \mathbf{A})\mathbf{B} = k(\mathbf{A} \mathbf{B}) = \mathbf{A}(k\mathbf{B})\)
2. \(\mathbf{A}(\mathbf{B}\mathbf{C}) = (\mathbf{A}\mathbf{B})\mathbf{C}\)
3. \((\mathbf{A} + \mathbf{B})\mathbf{C} = \mathbf{A}\mathbf{B} + \mathbf{A}\mathbf{C}\)
4. \(\mathbf{C}(\mathbf{A} + \mathbf{B}) = \mathbf{C}\mathbf{A} + \mathbf{C}\mathbf{A}\)

Exercise construct examples of each of these rules.

We can also multiply a matrix and vector. This is like the shapes of (m, r) * (r, 1) = (m, 1)

x = np.array([1, 2])
A @ x

array([ 5, 11])

There is a small subtle point, the x-array is 1-D:

x.shape

(2,)

Its shape is not (2, 1)! Numpy does the right thing here and figures out what you want. Not all languages allow this, however, and you have to be careful that everything has the right shape with them.

In the transpose operation you swap the rows and columns of an array. The transpose of A is denoted \(\mathbf{A}^T\).

A.T

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

There is also a function for transposing.

np.transpose(A)

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

A matrix is called symmetric if it is equal to its transpose: \(\mathbf{A} == \mathbf{A}^T\).

Q = np.array([[1, 2],
              [2, 4]])

np.allclose(Q, Q.T)

True

A matrix is called skew symmetric if \(\mathbf{A}^T = -\mathbf{A}\).

Q = np.array([[0, 1],
              [-1, 0]])

np.allclose(Q.T, -Q)

True

A matrix is called orthogonal if this equation is true: \(\mathbf{A} \mathbf{A}^T = \mathbf{I}\). Here is an example of an orthogonal matrix:

theta = 12
Q = np.array([[np.cos(theta), -np.sin(theta)],
              [np.sin(theta),  np.cos(theta)]])

Q @ Q.T

array([[1.00000000e+00, 2.19187673e-17],
       [2.19187673e-17, 1.00000000e+00]])

Here are the four rules for matrix multiplication and transposition

1. \((\mathbf{A}^T)^T = \mathbf{A}\)
2. \((\mathbf{A}+\mathbf{B})^T = \mathbf{A}^T+\mathbf{B}^T\)
3. \((\mathit{c}\mathbf{A})^T = \mathit{c}\mathbf{A}^T\)
4. \((\mathbf{AB})^T = \mathbf{B}^T\mathbf{A}^T\)

Exercise Come up with an example for each rule.

The determinant of a matrix is noted: det(A) or |A|. Many matrices are used to linearly transform vectors, and the determinant is related to the scaling magnitude.

np.linalg.det(A)

-2.0000000000000004

A matrix is invertible if and only if the determinant of the matrix is non-zero.

The inverse is defined by: \(\mathbf{A} \mathbf{A}^{-1} = \mathbf{I}\).

We compute the inverse as:

np.linalg.inv(A)

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

And here verify the definition.

A @ np.linalg.inv(A)

array([[1.00000000e+00, 1.11022302e-16],
       [0.00000000e+00, 1.00000000e+00]])

Another way to define an orthogonal matrix is \(\mathbf{A}^T = \mathbf{A}^{-1}\).

theta = 12
Q = np.array([[np.cos(theta), -np.sin(theta)],
              [np.sin(theta),  np.cos(theta)]])

np.allclose(Q.T, np.linalg.inv(Q))

True

The rank of a matrix is equal to the number of linearly independent rows in it. Rows are linearly independent if and only if they cannot be made by constants times another row or linear combinations of other rows.

np.linalg.matrix_rank(A)

2

Here is an example of a rank-deficient array. The last row is a linear combination of the first two rows.

A1 = [[1, 2, 3],
      [0, 2, 3],
      [2, 6, 9]]

np.linalg.matrix_rank(A1)

2

Here is an example of a rank-deficient array. It is deficient because the last row is just 0 times any other row.

A1 = [[1, 2, 3],
      [0, 2, 3],
      [0, 0, 0]]

np.linalg.matrix_rank(A1)

2

Note the determinant of this array is nearly zero as a result.

np.linalg.det(A1)

0.0

Also note the inverse has some enormous numbers in it. This is not a reliable inverse. It is never a good idea to have giant numbers and small numbers in the same calculations!

np.linalg.inv(A1)

LinAlgErrorTraceback (most recent call last) <ipython-input-34-ebbe6da30534> in <module> -—> 1 np.linalg.inv(A1)

<__array_function__ internals> in inv(*args, **kwargs)

~/opt/anaconda3/lib/python3.7/site-packages/numpy/linalg/linalg.py in inv(a) 545 signature = 'D->D' if isComplexType(t) else 'd->d' 546 extobj = get_linalg_error_extobj(_raise_linalgerror_singular) –> 547 ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj) 548 return wrap(ainv.astype(result_t, copy=False)) 549

~/opt/anaconda3/lib/python3.7/site-packages/numpy/linalg/linalg.py in _raise_linalgerror_singular(err, flag) 95 96 def _raise_linalgerror_singular(err, flag): —> 97 raise LinAlgError("Singular matrix") 98 99 def _raise_linalgerror_nonposdef(err, flag):

LinAlgError: Singular matrix

The condition number is a measure of the norm of an array times the inverse of the array. If it is very large, the array is said to be ill-conditioned.

np.linalg.cond(A1)

inf

What all of these mean is that we only have two independent rows in the array.