Unit 5: Matrices

I. Orientation — Matrix Computation in Python

A matrix is a rectangular arrangement of values with (m) rows and (n) columns; in Python, numerical matrices are usually represented by NumPy arrays rather than nested lists. NumPy and SciPy provide efficient storage, element-wise computation, linear algebra, random sampling, and sparse-matrix algorithms.

• Shape convention: A matrix (A\in\mathbb{R}^{m\times n}) contains (m) rows and (n) columns; an element is written (a_{ij}).
• Zero-based indexing: Python stores the first element at index 0, so mathematical (a_{ij}) corresponds to A[i-1, j-1] when (i,j) use one-based notation.
• Homogeneous storage: A NumPy array normally stores values of one data type, such as int64, float64, or complex128.
• Dimensional distinction:
  • A vector commonly has shape (n,).
  • A row matrix has shape (1, n).
  • A column matrix has shape (n, 1).
• Operation convention: A * B means element-wise multiplication, whereas A @ B means matrix multiplication.
• Performance principle: Array operations execute compiled numerical loops, making them generally faster and clearer than element-by-element Python loops.

II. Randomized Numerical Data — Controlled Sampling

A. Random number generation

Random number generation creates pseudorandom values for simulation, testing, initialization, and statistical sampling.

• Pseudorandom principle: A deterministic algorithm produces a sequence that behaves statistically like random data; the initial state is determined by a seed.
• Recommended generator: np.random.default_rng() creates an independent modern random-number generator.

import numpy as np

rng = np.random.default_rng(seed=42)
integers = rng.integers(1, 7, size=(2, 3))
uniform = rng.random((2, 3))
normal = rng.normal(loc=0.0, scale=1.0, size=5)

• Distribution meanings:
  • integers(1, 7) samples integers from (1) through (6); the upper bound is excluded.
  • random() samples uniformly from the interval ([0,1)).
  • normal(loc=μ, scale=σ) samples from a normal distribution with mean (\mu) and standard deviation (\sigma).
• Reproducibility: Reusing seed 42 reproduces the same sequence, which is useful in debugging and experiments.
• Sampling operations: rng.choice(a, size, replace=False) samples without replacement, while rng.shuffle(A) rearranges an array in place.
• Security limitation: NumPy generators are designed for numerical work, not passwords or cryptographic tokens; security-sensitive code should use Python’s secrets module.

III. Numerical Array Infrastructure — Efficient Matrix Storage

A. NumPy

NumPy is the core Python library for homogeneous multidimensional arrays and high-performance numerical operations.

• Central object: numpy.ndarray stores data in a fixed-dimensional, typed memory block.
• Construction methods:

A = np.array([[1, 2], [3, 4]], dtype=float)
Z = np.zeros((2, 3))
I = np.eye(3)
r = np.arange(0, 10, 2)
x = np.linspace(0, 1, 5)

• Concrete results: Z is a (2\times3) zero matrix, I is the (3\times3) identity matrix, and r equals [0, 2, 4, 6, 8].
• Data types: The dtype controls memory use and interpretation; for example, int32 generally uses 4 bytes per element and float64 uses 8 bytes.
• Type conversion: A.astype(np.int32) returns a converted copy; converting floating-point values to integers discards their fractional parts.
• Array advantage: Nested lists do not directly support matrix arithmetic, but NumPy arrays support broadcasting, reductions, and linear algebra.

IV. Element Selection — Accessing Matrix Regions

A. Indexing and slicing

Indexing selects individual entries, while slicing selects regular subarrays without usually copying the underlying data.

• Element access: For A = np.array([[10,20,30],[40,50,60]]), A[1,2] is 60.
• Slice syntax: start:stop:step includes start but excludes stop; therefore, A[:, 1:] selects every row and columns 1 onward.
• Negative indices: A[-1, -1] accesses the final row and final column.
• Dimensional effect:
  • A[0, :] has shape (3,).
  • A[0:1, :] preserves two dimensions and has shape (1, 3).
• Boolean indexing: A[A > 30] returns [40, 50, 60]; this is useful for filtering and conditional assignment.
• Fancy indexing: A[[0, 1], [2, 0]] selects coordinates (0,2) and (1,0), producing [30, 40].
• Views and copies: A basic slice such as B = A[:, :2] is usually a view, so changing B can change A; use .copy() when independent storage is required.

V. Array Metadata — Understanding Structure

A. Attributes of a NumPy array

Array attributes describe an array’s dimensions, element type, memory use, and layout.

• ndim: Gives the number of axes; a conventional matrix has A.ndim == 2.
• shape: Returns axis lengths; a (3\times4) matrix has A.shape == (3, 4).
• size: Counts all elements and equals the product of the shape dimensions: (3\times4=12).
• dtype: Identifies the stored type, such as dtype('float64').
• itemsize: Gives bytes per element; a float64 normally has itemsize == 8.
• nbytes: Gives element-buffer size and satisfies:

nbytes = size × itemsize

• Transpose attribute: A.T reverses the axes of a two-dimensional array, changing shape (m, n) to (n, m).
• Memory flags: A.flags reports properties such as C-contiguous row-major storage; these can influence interoperability and performance.

VI. Matrix Arithmetic — Computing with Arrays

A. Basic mathematical operations

NumPy distinguishes element-wise arithmetic, matrix products, reductions, and linear-algebra operations.

• Element-wise operators: For equally shaped arrays, A+B, A-B, A*B, A/B, and A**2 operate on corresponding entries.
• Matrix product: If (A) has shape ((m,n)) and (B) has shape ((n,p)), then (C=A@B) has shape ((m,p)), where:

cᵢⱼ = Σₖ aᵢₖbₖⱼ

Here (i) identifies a row, (j) a column, and (k) runs across the shared dimension.

A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
C = A @ B                 # [[19, 22], [43, 50]]

• Broadcasting: Arrays with compatible trailing dimensions combine without explicit repetition; adding shape (3,) to shape (2,3) adds the vector to each row.
• Reductions: A.sum(), A.mean(), A.min(), and A.max() reduce all entries; axis=0 processes columns and axis=1 processes rows.
• Linear algebra: np.linalg.solve(A, b) is preferable to np.linalg.inv(A) @ b for solving (Ax=b), because it is typically faster and numerically more stable.
• Numerical comparison: Use np.allclose(x, y) rather than exact equality for floating-point results affected by rounding.

VII. Reshaping and Combining — Structural Transformations

A. Array manipulation functions

Array manipulation functions change shape, orientation, ordering, or composition without changing the intended numerical information.

• Reshaping: A.reshape(2, 3) reorganizes six elements into two rows and three columns; the total element count must remain unchanged.
• Automatic dimension: In A.reshape(2, -1), NumPy infers the missing dimension.
• Flattening: A.ravel() returns a flattened view when possible, whereas A.flatten() always returns a copy.
• Joining arrays:
  • np.concatenate((A, B), axis=0) joins existing rows.
  • np.vstack((A, B)) stacks vertically.
  • np.hstack((A, B)) stacks horizontally when row counts match.
• Splitting: np.split(A, 2, axis=0) divides an axis into equal sections; unequal divisions require np.array_split.
• Axis adjustment: np.expand_dims(x, axis=1) can convert shape (n,) into (n,1), while np.squeeze() removes axes of length one.
• Reordering: np.transpose(A), np.flip(A, axis=0), and np.rot90(A) transpose, reverse, and rotate array data respectively.

VIII. Factorization — Revealing Matrix Structure

A. Matrix decomposition

Matrix decomposition expresses a matrix as a product of simpler matrices for solving systems, reducing dimensions, and analyzing numerical properties.

• LU decomposition: Represents (A=PLU), where (P) is a permutation matrix, (L) is lower triangular, and (U) is upper triangular; pivoting improves stability.
• QR decomposition: Represents (A=QR), where (Q^TQ=I) and (R) is upper triangular; it is central to least-squares calculations.

Q, R = np.linalg.qr(A)

• Eigenvalue decomposition: For suitable square matrices, (A=V\Lambda V^{-1}); columns of (V) are eigenvectors and diagonal entries of (\Lambda) are eigenvalues satisfying (Av=\lambda v).
• Singular value decomposition: Every (m\times n) matrix has:

A = UΣVᴴ

Here (U) and (V) are unitary or orthogonal matrices, (\Sigma) contains nonnegative singular values, and (V^H) is the conjugate transpose.

• Cholesky decomposition: A symmetric positive-definite matrix can be written (A=LL^T); this is efficient for relevant linear systems.
• Practical rule: Use decomposition-based solvers rather than explicitly computing inverses, especially for large or ill-conditioned matrices.

IX. Efficient Storage of Mostly Zero Data

A. Sparse matrices and associated data structures

Sparse matrices store only significant nonzero entries, reducing memory and computation when zeros dominate.

• COO format: Coordinate storage keeps parallel arrays of row indices, column indices, and values; it is convenient for constructing a matrix from triplets ((i,j,a_{ij})).
• CSR format: Compressed Sparse Row uses data, indices, and indptr; it supports efficient row slicing and matrix-vector multiplication.
• CSC format: Compressed Sparse Column is analogous to CSR but favors column operations.
• DOK and LIL formats: Dictionary of Keys and List of Lists support incremental insertion; they are commonly converted to CSR or CSC before computation.
• Construction example:

from scipy.sparse import csr_matrix

S = csr_matrix(([5, 8], ([0, 2], [1, 2])), shape=(3, 3))

This stores nonzero values (S{0,1}=5) and (S{2,2}=8).

• Limitation: Sparse storage is inefficient for dense matrices, and calling S.toarray() may consume prohibitive memory for large dimensions.

X. Scientific Algorithms — Extending NumPy

A. SciPy

SciPy builds on NumPy by providing specialized algorithms for scientific computing.

• Linear algebra: scipy.linalg supplies decompositions, matrix functions, and structured solvers beyond NumPy’s core facilities.
• Sparse computation: scipy.sparse defines sparse formats, while scipy.sparse.linalg provides routines such as spsolve, cg, and sparse eigenvalue solvers.
• Optimization: scipy.optimize supports root finding, minimization, curve fitting, and constrained optimization.
• Integration: scipy.integrate.quad approximates definite integrals, and solve_ivp handles initial-value differential equations.
• Statistics and signals: scipy.stats implements probability distributions and tests; scipy.signal supports filtering and spectral operations.
• Interoperation: SciPy functions usually accept NumPy arrays and return arrays or related scientific objects, allowing both libraries to form one numerical workflow.

XI. Eliminating Python-Level Loops — Array-Wide Computation

A. Vectorization of code

Vectorization expresses computations as whole-array operations so that optimized compiled routines perform the repeated work.

• Loop replacement: Squaring one million values is written as y = x**2, not as an explicit Python loop with repeated append.
• Concrete transformation:

# Loop-based
result = np.empty_like(x)
for i in range(x.size):
    result[i] = 3 * x[i] + 2

# Vectorized
result = 3 * x + 2

• Conditional computation: np.where(A >= 0, A, 0) replaces negative entries with zero without manually traversing indices.
• Broadcasted outer operation: If x has shape (m,1) and y has shape (1,n), then x + y produces an (m,n) matrix of pairwise sums.
• Universal functions: Operations such as np.sin, np.exp, and np.sqrt apply element-wise and support reductions through methods such as np.add.reduce.
• Benefits: Vectorized code is usually shorter, less error-prone, and faster because it minimizes interpreted Python operations.
• Trade-off: Broadcasting can create large temporary arrays; in memory-sensitive work, use in-place operations, chunking, np.einsum, or specialized matrix routines.