Linear algebra

Matrix multiplication

There are two multiplications, and confusing them is a classic mistake. The * multiplies element by element. The @ (or np.matmul) is the multiplication of matrices, rows by columns:

import numpy as np

A = np.array([[2.0, 1.0], [1.0, 3.0]])
b = np.array([3.0, 5.0])
print(A @ A)
print(A.T @ b)
print(A * A)
print(np.array_equal(A * A, A @ A))

[[ 5.  5.]
 [ 5. 10.]]
[11. 18.]
[[4. 1.]
 [1. 9.]]
False

For vectors, the @ is the dot product. With it and with the norm (the length of the vector), the cosine of the angle between two vectors is calculated, a measure of similarity used in search and recommendation:

u = np.array([3.0, 4.0])
v = np.array([1.0, 0.0])
print(u @ v, np.linalg.norm(u), (u @ v) / (np.linalg.norm(u) * np.linalg.norm(v)))

3.0 5.0 0.6

Solve a system of equations

The system 2x + y = 3 and x + 3y = 5 is written as A @ x = b. The np.linalg.solve it solves, and checking is simple: multiply back and compare:

x = np.linalg.solve(A, b)
print(x, np.allclose(A @ x, b))

[0.8 1.4] True

Do not calculate the inverse to solve systems

Many people write inv(A) @ b. Avoid: calculating the inverse is slower and less accurate than the solve. Use the solve whenever the goal is to solve a system, and leave the inverse for when you really need it.

Determinant, inverse and eigenvalues

print(round(np.linalg.det(A), 6))
print(np.round(np.linalg.inv(A), 3))
valores, vetores = np.linalg.eigh(A)
print(np.round(valores, 3))

5.0
[[ 0.6 -0.2]
 [-0.2  0.4]]
[1.382 3.618]

The eigh is the version for symmetricmatrices, faster and more stable than the general one (eig). The eigenvalues summarize "how much the matrix stretches" each direction, and are behind techniques like principal component analysis.

Fit a line: least squares

Given a set of points with noise, the np.linalg.lstsq finds the line that minimizes the squared error. Just set up the matrix with the column of x and a column of ones (for the intercept):

rng = np.random.default_rng(1)
xs = np.linspace(0, 10, 50)
ys = 3 * xs + 2 + rng.normal(0, 0.5, 50)
X = np.column_stack([xs, np.ones_like(xs)])
(coef, intercepto), *_ = np.linalg.lstsq(X, ys, rcond=None)
print(abs(coef - 3) < 0.1, abs(intercepto - 2) < 0.5)

True True

I generated the data with slope 3 and intercept 2, plus noise, and the fit retrieves the with good accuracy. It is the same calculation that chapter 40 compares with scikit-learn. works in

Batches of matrices

The @ batches: it multiplies the last two dimensions and treats the previous ones as a batch (with broadcasting). This is how a neural network processes multiple examples at once: Projecting one vector onto another, the projection of

lote_a = np.ones((10, 3, 4))
lote_b = np.ones((10, 4, 2))
print((lote_a @ lote_b).shape)

(10, 3, 2)