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.