7.3 Matrices in Python

Source: Ideas from Computational Physics (Landau, Páez, Bordeianu), §7.3–7.3.3. Rewritten in original wording. Continues 6.8b Nonlinear least-squares labs. Labs: 7.3b Matrices in Python labs.

Lists as arrays

Python's built-in sequence type is the list — an ordered, mutable container:

L = [1, 2, 3]
L[0]        # 1 (0-indexed)
L[0] = 5    # L is now [5, 2, 3]
len(L)      # 3

A tuple looks similar but is immutable — T = (1, 2, 3, 4); attempting T[0] = 5 raises an error. Lists are dynamic (no need to declare a size up front), and a "list of lists" gives compound structure: L = [[1,2],[3,4],[5,6]] with L[0] giving [1, 2].

Handy list operations: L1 + L2 concatenates, x in L tests membership, L[i:j] slices, L.append(x), L.sort(), L.reverse(), L.count(x), L.index(x), L.remove(x).

Why NumPy instead of raw lists

Plain lists are flexible but not built for math: [1,2,3] + [1,2,3] concatenates to a 6-element list rather than adding elementwise, and there is no real notion of a mathematical matrix. NumPy arrays fix this — a uniform (single dtype), fixed-shape container optimized for vectorized numerical operations:

from numpy import *
vector1 = array([1, 2, 3, 4, 5])
vector2 = vector1 + vector1     # [2, 4, 6, 8, 10] — real elementwise addition
vector2 = 3 * vector1           # [3, 6, 9, 12, 15] — scalar multiply

A NumPy array can have up to 32 dimensions (indices), but every element must share the same type — mixing int and float input upcasts the whole array to float64. The number of indices is the array's ndim; what mathematicians call a matrix's "size" is the array's shape; size gives the total element count.

import numpy as np
a = np.arange(12).reshape((3, 4))  # 0..11 reshaped into 3 rows, 4 cols
a.shape   # (3, 4)
a.ndim    # 2
a.size    # 12
a.T       # transpose: 4 rows, 3 cols

Slicing extracts sub-arrays with start:stop:step, and — critically — indices count from 0, so 1:3 means indices 1 and 2 (not 3):

a[:2, :]     # first two rows
a[:, 1:3]    # columns 1 and 2

Slicing avoids copying large chunks of memory and is one reason vectorized NumPy code runs much faster than an equivalent Python loop.

Elementwise * vs. the real matrix product

This is the single most common NumPy pitfall: * between two 2D arrays multiplies elementwise, not as a mathematical matrix product. For real matrix multiplication, use dot (or the @ operator):

matrix1 = array([[0, 1], [1, 3]])
dot(matrix1, matrix1)     # [[1, 3], [3, 10]]  — true matrix product
matrix1 * matrix1         # [[0, 1], [1, 9]]   — elementwise product

NumPy also supports broadcasting — assigning or operating on many elements with a single statement:

w = zeros(100, float)
w = 23.7          # every element of w becomes 23.7

stuff = zeros(10, float)
t = arange(4)
stuff[3:7] = sqrt(t + 1)   # sqrt() applied elementwise, result written into a slice

Broadcasting lets Python push the loop into fast C code instead of iterating element-by-element in the interpreter — this is where most of NumPy's speed comes from.

Solving linear systems with NumPy's linear algebra library

NumPy arrays by themselves are not mathematical matrices — for that, use numpy.linalg, a thin, robust interface to the battle-tested LAPACK library. The standard problem is

$$A\mathbf x=\mathbf b,$$

with $A$ a known $N\times N$ matrix, $\mathbf b$ a known vector, and $\mathbf x$ unknown:

from numpy import array
from numpy.linalg import solve

A = array([[1, 2, 3], [22, 32, 42], [55, 66, 100]])
b = array([1, 2, 3])
x = solve(A, b)              # solves Ax=b directly (Gaussian elimination / LU)
residual = dot(A, x) - b     # ~ machine precision, essentially zero

An alternative is computing the explicit inverse and multiplying: x = dot(inv(A), b). This works, but solve is both faster and more numerically stable — computing a full inverse does more work than necessary and can amplify round-off error. Prefer solve unless you specifically need the inverse matrix itself for something else.

Eigenvalue problems $A\boldsymbol\omega=\lambda\boldsymbol\omega$ (e.g. finding the principal axes of an inertia tensor) are solved with numpy.linalg.eig:

from numpy.linalg import eig

I = array([[2/3, -1/4], [-1/4, 2/3]])
eigenvalues, eigenvectors = eig(I)

Each column of eigenvectors is an eigenvector paired with the corresponding entry of eigenvalues; the pairing is verified by checking that $I\boldsymbol\omega-\lambda\boldsymbol\omega$ comes out to (near) zero.

Why use a library at all? Hand-written matrix routines are rarely competitive: library implementations (LAPACK and friends) are typically an order of magnitude faster, are engineered to minimize round-off error, and transparently adapt to the local hardware (desktop vs. cluster). Writing your own Gaussian elimination is valuable for understanding the algorithm (as in the spline and quadratic-fit exercises earlier in this chapter) — but production numerical work should reach for numpy.linalg.

Quiz

Matrices in Python

1 / 4

A NumPy array, unlike a Python list, requires that…

Pick an answer to see if you are right.

Next: 7.3b Matrices in Python labs.