Matrix Operations

A 2-D NumPy array is a matrix, and the matrix toolkit sits on top of the array primitives already covered: products, transposes, inverses, and the special-matrix constructors. The discipline that keeps matrix code correct is treating * as elementwise (it is not matrix multiplication) and reserving @ for the true matrix product. This chapter covers the constructors and the structural operations; solving and decompositions are in the linear algebra chapters.

Special Matrix Constructors

import numpy as np

I = np.eye(3)                 # 3x3 identity
Z = np.zeros((2, 3))          # zero matrix
D = np.diag([1, 2, 3])        # diagonal matrix from a vector
v = np.diag(D)                # and back: extract the diagonal
print(D.shape, v)

print(np.triu(D + np.ones((3, 3))))   # upper triangular (diagonal included)
print(np.tril(D + np.ones((3, 3))))   # lower triangular

# a matrix from a function over its indices (row, col)
M = np.fromfunction(lambda r, c: r * 10 + c, (2, 3)).astype(int)
print(M)                       # [[ 0  1  2] [10 11 12]]

eye, zeros, and full are the standard initial matrices; diag converts between vectors and diagonal matrices (useful for implementing D @ x scaling or diagonal regularization). The triangular constructors mask a matrix into its upper or lower part - the building blocks for iterative solvers and banded systems.

Products and Shape Contracts

A = np.array([[1, 2], [3, 4]])       # (2, 2)
B = np.ones((2, 3))                   # (2, 3)
v = np.array([1.0, 1.0])              # (2,)

print(A @ B)            # (2, 3) - matrix product, defined when inner dims match
print(A @ v)            # (3,) - matrix times vector
print(A.T @ A)          # (2, 2) - always defined for any shape, symmetric PSD

# the elementwise trap
print(A * A)            # elementwise - NOT the matrix product
# rule: * is elementwise, @ is the matrix product

The shape contract for A @ B: if A is (m, k) and B is (k, n), the result is (m, n) and the inner dimensions must be equal. A shape mismatch here is the most common matrix error, and the fix is almost always a missing transpose - print the shapes of both operands and the intended result, then check the contract before reading the traceback.

Transposes and Symmetry

A = np.random.default_rng(0).normal(size=(4, 3))
print(A.T.shape)              # (3, 4) - a free view, no data moved

S = A.T @ A                   # (3, 3) - the Gram matrix
print(np.allclose(S, S.T))    # True - always symmetric

# checking symmetry numerically, never with ==
diff = np.abs(S - S.T).max()
print(diff < 1e-12)           # True up to roundoff

Transposing is a view, so large matrix transposes are free. Symmetry checks must be numerical - allclose or a max-difference bound - because floating-point products introduce roundoff that an exact == would flag.

Structural Operations

A = np.arange(16, dtype=float).reshape(4, 4)

print(A[:, :2])               # first two columns
print(A[1:, 2:])              # trailing block

print(A.sum(axis=0))          # column sums (4,)
print(A.sum(axis=1))          # row sums (4,)
print(A.max(axis=1)[:, None]) # per-row max, kept as a column for broadcasting

# normalization rows: divide each row by its sum - pure broadcasting
row_sums = A.sum(axis=1, keepdims=True)
P = A / row_sums
print(P.sum(axis=1))          # [1. 1. 1. 1.] - each row normalized

Row/column selection plus axis aggregation plus broadcasting covers the day-to-day matrix processing: extracting blocks, computing per-row statistics, and normalizing. The keepdims=True detail is what makes per-row scaling one line - the result keeps its shape so it broadcasts back against the original.

Summary

Build matrices with the constructors (eye, zeros, diag, the triangular family), multiply with @ (never *), and check the (m,k)x(k,n) shape contract whenever a product fails. Transpose freely - it is a view - verify symmetry numerically with allclose, and use axis aggregation with keepdims plus broadcasting for per-row operations like normalization. With those habits, the remaining matrix machinery - solving and decompositions - slots in on the same array foundation.