Matrix Multiplication and np.linalg

So far every operation worked element by element. Linear algebra adds operations that work on whole matrices: matrix multiplication, inverses, solving equations and more. They are the maths behind graphics, machine learning, statistics and engineering. NumPy puts them in the @ operator and the np.linalg module.

How to run the examples

Create a new file named linear_algebra.py in your NumpyTutorials folder (see the setup in lesson 01). Any name works; this one just matches the topic. Start the file with these lines:

import numpy as np

Then copy the examples of this page into the same file, one after another, in the order shown, and run it with python linear_algebra.py (with your virtual environment active). Some examples use variables created in earlier ones, so keep everything in the one file. Each Output box shows what that example prints.

1. * vs @

A = np.array([[1, 2],
              [3, 4]])
B = np.array([[5, 6],
              [7, 8]])

print("A * B (element by element):\n", A * B)
print("A @ B (matrix multiplication):\n", A @ B)
# A @ B, row 0 col 0 = 1*5 + 2*7 = 19

Output

A * B (element by element):
 [[ 5 12]
 [21 32]]
A @ B (matrix multiplication):
 [[19 22]
 [43 50]]

A * B multiplies matching elements: 1x5, 2x6, 3x7, 4x8.
A @ B is matrix multiplication: each result is a row of A times a column of B, added up. Top-left: 1x5 + 2x7 = 19. Top-right: 1x6 + 2x8 = 22.

Mixing up * and @ is the most common linear-algebra bug. Both work on two square matrices without any error, but they give completely different answers.

2. @, matmul and dot

print("np.matmul(A, B):\n", np.matmul(A, B))
print("np.dot(A, B):\n", np.dot(A, B))

v = np.array([1, 2, 3])
w = np.array([4, 5, 6])
print("dot of two 1D arrays v @ w:", v @ w, "(1*4 + 2*5 + 3*6)")

Output

np.matmul(A, B):
 [[19 22]
 [43 50]]
np.dot(A, B):
 [[19 22]
 [43 50]]
dot of two 1D arrays v @ w: 32 (1*4 + 2*5 + 3*6)

A @ B, np.matmul(A, B) and np.dot(A, B) are the same for 2D matrices. For two 1D arrays, @ gives the dot product, a single number: 1x4 + 2x5 + 3x6 = 32. Prefer @ in new code; it is the shortest and clearest.

3. Shapes: (m, n) @ (n, p) gives (m, p)

C = np.arange(6).reshape(2, 3)      # 2x3
D = np.arange(6).reshape(3, 2)      # 3x2
print("(2, 3) @ (3, 2) ->", (C @ D).shape, "\n", C @ D)
print("(3, 2) @ (2, 3) ->", (D @ C).shape, "(order matters!)")

try:
    C @ C                               # (2, 3) @ (2, 3): 3 != 2
except ValueError as error:
    print("Error:", error)

Output

(2, 3) @ (3, 2) -> (2, 2) 
 [[10 13]
 [28 40]]
(3, 2) @ (2, 3) -> (3, 3) (order matters!)
Error: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 2 is different from 3)

The inner sizes must match: (2, 3) @ (3, 2) works because 3 = 3.
The result takes the outer sizes: (2, 2).
Order matters: (3, 2) @ (2, 3) gives a (3, 3) matrix. In general A @ B is not equal to B @ A.
(2, 3) @ (2, 3) fails because 3 and 2 don't match.

4. Real example: bill of 3 customers

# rows = customers, columns = products (pen, notebook, bag)
quantities = np.array([[10, 2, 0],
                       [0, 5, 1],
                       [3, 3, 3]])
prices = np.array([10, 60, 450])            # price of each product

print("bill per customer (quantities @ prices):", quantities @ prices)

# two price lists at once: columns = normal price, sale price
price_lists = np.array([[10, 8],
                        [60, 50],
                        [450, 400]])
print("bill per customer for each price list:\n", quantities @ price_lists)

Output

bill per customer (quantities @ prices): [ 220  750 1560]
bill per customer for each price list:
 [[ 220  180]
 [ 750  650]
 [1560 1374]]

Rows are customers, columns are products. quantities @ prices multiplies each customer's quantities by the prices and adds them up: customer 0 pays 10x10 + 2x60 + 0x450 = 220. With a matrix of two price lists, one @ calculates every customer's bill under both lists at once.

5. np.linalg.inv: the inverse

The inverse of a matrix M is the matrix that "undoes" it: M @ inverse = identity (the "1" of matrices, lesson 02).

M = np.array([[4.0, 7.0],
              [2.0, 6.0]])
M_inv = np.linalg.inv(M)
print("inverse:\n", M_inv)
print("M @ inverse (should be identity):\n", (M @ M_inv).round(10))
print("allclose to eye(2):", np.allclose(M @ M_inv, np.eye(2)))

Output

inverse:
 [[ 0.6 -0.7]
 [-0.2  0.4]]
M @ inverse (should be identity):
 [[ 1. -0.]
 [ 0.  1.]]
allclose to eye(2): True

The -0. in the result is a tiny rounding leftover such as -0.0000000000000001, rounded. Floats are rarely exact, so compare matrices with np.allclose, not == (lesson 14).

6. np.linalg.det: the determinant

print("det(M)        :", round(np.linalg.det(M), 4), "(4*6 - 7*2 = 10)")
singular = np.array([[1.0, 2.0],
                     [2.0, 4.0]])           # row 2 = 2 x row 1
print("det(singular) :", np.linalg.det(singular))

try:
    np.linalg.inv(singular)
except np.linalg.LinAlgError as error:
    print("Error:", error, "-> a matrix with det 0 has no inverse")

Output

det(M)        : 10.0 (4*6 - 7*2 = 10)
det(singular) : 0.0
Error: Singular matrix -> a matrix with det 0 has no inverse

For a 2x2 matrix [[a, b], [c, d]] the determinant is a x d - b x c, so 4x6 - 7x2 = 10. When the determinant is 0, the matrix is singular: one row is just a multiple of another, it has no inverse, and NumPy raises LinAlgError.

7. np.linalg.solve: solving equations

A classic school problem: 2 pens and 3 notebooks cost 190; 4 pens and 1 notebook cost 130. What does each cost? Write the numbers in front of the unknowns as a matrix and the totals as a vector, and solve does the rest.

# 2 pens + 3 notebooks = 190
# 4 pens + 1 notebook  = 130      What does one pen / notebook cost?
coefficients = np.array([[2, 3],
                         [4, 1]])
totals = np.array([190, 130])

pen, notebook = np.linalg.solve(coefficients, totals)
print(f"pen = {pen:.0f}, notebook = {notebook:.0f}")
print("check:", coefficients @ np.array([pen, notebook]))

Output

pen = 20, notebook = 50
check: [190. 130.]

A pen costs 20 and a notebook 50. Multiplying back with @ gives the original totals, so the answer is right. solve works for any number of equations, and it is faster and more accurate than calculating inv(A) @ b yourself.

8. np.linalg.eig and eigh

An eigenvector of a matrix is a direction that the matrix only stretches, without turning it. The eigenvalue is how much it is stretched. They are used in physics, Google's PageRank and PCA in machine learning.

S = np.array([[2.0, 1.0],
              [1.0, 2.0]])                   # symmetric: S equals S.T

print("eig  values:", np.linalg.eig(S).eigenvalues, "(eig works for any square matrix, may give complex numbers)")

values, vectors = np.linalg.eigh(S)          # eigh: for symmetric matrices, real results
print("eigh values:", values)
print("eigenvectors (columns):\n", vectors.round(4))

last = vectors[:, -1]                        # eigenvector of the biggest value (3)
print("S @ v       :", (S @ last).round(4))
print("lambda * v  :", (values[-1] * last).round(4), "<- the same: S only stretches v")

Output

eig  values: [3.+0.j 1.+0.j] (eig works for any square matrix, may give complex numbers)
eigh values: [1. 3.]
eigenvectors (columns):
 [[-0.7071  0.7071]
 [ 0.7071  0.7071]]
S @ v       : [2.1213 2.1213]
lambda * v  : [2.1213 2.1213] <- the same: S only stretches v

np.linalg.eig works for any square matrix. Its results can be complex numbers, which is why it prints 3.+0.j here.
When the matrix is symmetric (it equals its own transpose), use np.linalg.eigh: it is faster and gives real numbers, sorted from small to big.
The eigenvectors are the columns of the second result.
The check: S @ v equals 3 * v. The matrix only stretched v by its eigenvalue, 3.

9. np.linalg.norm: length and distance

v = np.array([3, 4])
print("length of [3, 4]          :", np.linalg.norm(v))
print("distance (1,1) to (4,5)   :", np.linalg.norm(np.array([4, 5]) - np.array([1, 1])))
print("L1 norm |3| + |4|         :", np.linalg.norm(v, ord=1))
print("matrix (Frobenius) norm   :", np.linalg.norm(A).round(4))
print("unit vector v / norm(v)   :", v / np.linalg.norm(v))

Output

length of [3, 4]          : 5.0
distance (1,1) to (4,5)   : 5.0
L1 norm |3| + |4|         : 7.0
matrix (Frobenius) norm   : 5.4772
unit vector v / norm(v)   : [0.6 0.8]

The norm of a vector is its length: √(3² + 4²) = 5 (Pythagoras).
The norm of the difference of two points is the distance between them.
ord=1 gives the "city block" length |3| + |4| = 7.
Dividing a vector by its norm gives a unit vector (length 1) pointing the same way.

10. np.linalg.svd

SVD (singular value decomposition) splits any matrix into three parts, U, s and Vt, that multiply back to the original. The values in s are sorted by importance. Keeping only the biggest ones gives a simpler matrix that still captures the main pattern. This is the idea behind recommendation systems and image compression.

# ratings: 4 users x 3 movies
ratings = np.array([[5.0, 4.0, 1.0],
                    [4.0, 5.0, 1.0],
                    [1.0, 1.0, 5.0],
                    [2.0, 1.0, 4.0]])
U, s, Vt = np.linalg.svd(ratings, full_matrices=False)
print("U shape", U.shape, " s", s.round(3), " Vt shape", Vt.shape)
print("rebuilt equals original:", np.allclose(U @ np.diag(s) @ Vt, ratings))

# Keep only the biggest value -> a simpler "summary" of the table
rank1 = s[0] * np.outer(U[:, 0], Vt[0])
print("approximation with 1 value:\n", rank1.round(1))

Output

U shape (4, 3)  s [10.065  5.421  1.143]  Vt shape (3, 3)
rebuilt equals original: True
approximation with 1 value:
 [[4.  3.8 2.7]
 [4.  3.8 2.7]
 [2.3 2.1 1.5]
 [2.4 2.3 1.6]]

The ratings table has two clear groups: users 0 and 1 like movies 0 and 1, and users 2 and 3 like movie 2. U @ np.diag(s) @ Vt rebuilds the original exactly. Using only the first and biggest value (10.065) gives an approximation that keeps the main "taste" pattern in a single number per user and per movie.

Exercises

Try these exercises on your own before moving on:

Multiply a 3x2 and a 2x4 matrix. What is the shape of the result?
Solve: x + y = 10 and x - y = 2.
Find the inverse of [[2, 1], [5, 3]] and check that M @ inv = identity.
Find the distance between the points (2, 3, 6) and (0, 0, 0).
Why does np.linalg.inv fail for [[1, 2], [3, 6]]?

Summary

* multiplies element by element; @ (or np.matmul/np.dot) is matrix multiplication.
(m, n) @ (n, p) gives (m, p); the inner sizes must match and order matters.
inv gives the inverse and det the determinant; a matrix with det 0 has no inverse.
np.linalg.solve(A, b) solves equations; prefer it to inv(A) @ b.
eig/eigh give eigenvalues and eigenvectors; use eigh for symmetric matrices.
norm measures length and distance; svd splits a matrix into its most important parts.

Next: Lesson 14: Missing Data, Rounding, clip, isclose.