```
import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg
```

## Matrix multiplication

All operations on an array are element-wise. Numpy used to have a "matrix" mode, where all operations were "matrix-wise"; that is, something like `*`

would do matrix multiplication instead of element-wise multiplication. This was ugly and messy, and has been replaced in Python 3.5+ with a matrix multipy operator, `@`

. (Older Python: Use `.matmul()`

or `.dot()`

.)

Let's first look at the diminsion rules for matrix multiplicaiton:

`[a, b] @ [b, c] = [a, c]`

```
(np.ones([3, 4]) @ np.ones([4, 5])).shape
```

The "inner" diminsions go away. This works for ND arrays, too:

`[a] @ [a] = scalar`

```
(np.ones([4]) @ np.ones([4])).shape
```

One of the two is allowed to have more than 2 dimensions, in which case it behaves like "stacks" of arrays:

`[a,b,c] @ [c,d] = [a,b,d]`

```
(np.ones([2, 3, 4]) @ np.ones([4, 5])).shape
```

Normal "prepend 1" broadcasting rules apply.

```
np.array([1, 2, 3]) @ np.array([[1, 2, 3]]).T
```

### Power user: Einstein summation notation

You can use Einstein summation notation for full control:

```
a = np.arange(25).reshape(5, 5)
```

```
np.trace(a)
```

```
np.einsum("ii", a)
```

```
a.T
```

```
np.einsum("ji", a)
```

```
a @ a
```

```
np.einsum("ij,jk", a, a)
```

```
np.sum(a * a)
```

```
np.einsum("ij,ij", a, a)
```

```
np.einsum("ij->", a ** 2)
```

```
b = np.array([1, 2, 3])
print(b)
```

```
A = np.array([[1, 2, 3], [22, 32, 42], [55, 66, 100]])
print(A)
```

```
np.linalg.inv(A) @ b
```

Note that for these equations, 1D vectors really should be 2D column vectors! `@`

and solve handle 1D vectors pretty well so we are safe, but be careful.

Computing the inverse is slow - there are faster algorithms when you just want to solve one case, available as `solve`

and internally using the LAPACK matrix library. We can even tell solve if we know something special about our matrix, like if we have a diagonal matrix, if we use `scipy.linalg.solve`

instead!

```
x = np.linalg.solve(A, b)
```

```
A @ x - b
```

```
x
```

```
A = np.array([[4, -2, 1], [3, 6, -4], [2, 1, 8]])
```

```
np.linalg.inv(A) - 1 / 263 * np.array([[52, 17, 2], [-32, 30, 19], [-9, -8, +30]])
```

```
x1 = np.array([12, -25, 32])
x2 = np.array([4, -10, 22])
x3 = np.array([20, -30, 40])
```

```
print(np.linalg.inv(A) @ x1)
print(np.linalg.inv(A) @ x2)
print(np.linalg.inv(A) @ x3)
```

```
print(np.linalg.solve(A, x1))
print(np.linalg.solve(A, x2))
print(np.linalg.solve(A, x3))
```

For such a tiny problem, `inv`

beats `solve`

by a hair. But if you invert only once, you can solve many problems with the same solution!

```
%%timeit
np.linalg.solve(A, x1)
np.linalg.solve(A, x2)
np.linalg.solve(A, x3)
```

```
%%timeit
np.linalg.inv(A) @ x1
np.linalg.inv(A) @ x2
np.linalg.inv(A) @ x3
```

```
%%timeit
Ai = np.linalg.inv(A)
Ai @ x1
Ai @ x2
Ai @ x3
```

```
a = scipy.linalg.hilbert(8)
```

Or, given the formula in the book we could have produced it ourselves:

```
i, j = np.ogrid[1 : len(a) + 1, 1 : len(a) + 1]
a_ours = 1 / (i + j - 1)
```

```
np.all(a == a_ours)
```

We need `b`

, which is just the first row of the matrix:

```
b = a[0]
```

Let's try solve. When you try 100x100, does this still work?

```
np.linalg.solve(a, b)
```

We can use the invhilbert function to make an inverse hilbert matrix. You can pass `exact=True`

to return integers instead of double floats. Note that this matrix will overflow 64 bit integers at 14x14, and therefore will become an inefficient python integer array.

```
scipy.linalg.invhilbert(8) @ b
```

We can also take the inverse ourselves. What happens when this becomes larger?

```
np.linalg.inv(a) @ b
```

```
I = 1 / 12 * np.array([[8, -3, -3], [-3, 8, -3], [-3, -3, 8]])
```

```
λs, ωs = np.linalg.eig(I)
```

```
λs
```

```
ωs
```

```
ω = ωs[:, 0]
print("ω:", ω)
print("Iω", I @ ω)
print("λω", λs[0] * ω)
```

```
print(I @ ωs - λs * ωs)
```

Mathematics tends to think of stacked vectors as column vectors - which is the opposite of what's natural to write in a text file, like `[[row1], [row2]]`

. This is why matrix oriented languages tend to be column major. You can stack in rows, use `.T`

when needed for matrix manipulations. You might want to start in "F" order if performance is important.