8  Python and NumPy Primer

Condensed reference for the toolset used throughout these notes.

8.1 Arrays and Basic Operations

NoteDefinition: NumPy Arrays

The ndarray is the fundamental structure.

  • Creation: np.array([1, 2]), np.zeros((m, n)), np.random.randn(m, n).

  • Indexing: 0-indexed. A[i, j] for entries; A[i, :] for rows; A[:, j] for columns.

TipRemark

(Shapes and Rank-1 Arrays) A.shape returns (m, n). Note that a 1D array has shape (n,), which is not equivalent to a column vector (n, 1). These `rank-1 arrays'' are a common source of broadcasting bugs. Usex[:, None]orx.reshape(-1, 1)` to convert to 2D.

WarningExercise

  1. Create \(4 \times 4\) identity matrix. Verify shape and dtype.

  2. Extract row 2 and column 3 from a random \(3 \times 3\) \(\bA\).

  3. Confirm that np.zeros((3,)) fails where a (3, 1) column vector is required in a matrix-vector product.

8.2 Linear Algebra Operations

Core routines in numpy.linalg:

TipRemark

(Performance Law: Vectorization) Explicit for loops in Python are slow. Push heavy lifting into NumPy/BLAS routines by vectorizing operations. For example, use A @ B instead of triple-nested loops.

TipRemark

(Stability: solve vs inv) np.linalg.solve(A, b) is faster and more stable than np.linalg.inv(A) @ b. Explicit inversion is almost never desirable in scientific computation.

TipRemark

(Computed equality) Use np.allclose(x, y) or np.isclose(x, y) for computed floating-point quantities. Exact equality is appropriate for integers, booleans, and deliberately exact structural checks.

TipRemark

(Sparse matrices) Large discretized PDE and graph problems should usually use scipy.sparse, not dense NumPy arrays. Use scipy.sparse.linalg for sparse solves, eigenvalue computations, CG, and GMRES.

WarningExercise

  1. Solve a \(2 \times 2\) system and verify the residual \(\|\bA\bx - \bb\|\).

  2. Profile solve vs inv @ b for \(n=1000\).

  3. Compare np.linalg.eig and np.linalg.eigh (optimized for symmetric matrices).

8.3 Key Idioms

  • Element-wise: * and ** are entry-wise. Matrix powers use np.linalg.matrix\_power.

  • Broadcasting: Automatic dimension expansion. A * v scales columns of \(\bA\) by entries of \(v\).

  • Aliases, Views, Copies: B = A makes a second name for the same array. Slices are usually views. Use B = A.copy() for independent data.

  • Sparse construction: use scipy.sparse.diags for banded matrices and scipy.sparse.csr\_matrix for efficient arithmetic.

WarningExercise

  1. Let \(\bA = \begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix}\). Compare A * A and A @ A.

  2. Demonstrate the ``view pitfall’’ by modifying a slice of a matrix.

  3. Use broadcasting to perform diagonal scaling \(\bA\bD\) without constructing \(\bD = \text{diag}(d)\).