Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Activity 5.1

This page is an ipython notebook, meaning you can download it and run the executable cells. It directly follows Activity 5.1. Reference answers appear after the questions, challenges, predictions, and reflections. They show one reasonable solution or interpretation, but other correct approaches may exist.

Matrix Computation, Broadcasting, Vectorization, QR & SVD

Tools: Python, NumPy, Jupyter/VS Code

Learning objectives

By the end of this tutorial, you should be able to:

  • create and inspect NumPy arrays;

  • reason about array shape;

  • distinguish one-dimensional, row, and column vectors;

  • replace simple Python loops with vectorized operations;

  • apply broadcasting rules;

  • use axis correctly;

  • distinguish element-wise and matrix multiplication;

  • solve linear systems numerically;

  • compute and verify QR decompositions;

  • compute and verify SVD decompositions;

  • build low-rank approximations;

  • connect SVD to dimensionality reduction;

  • investigate conditioning and multicollinearity.

Part 0 — Set up

Create:

notebooks/week09_linear_algebra.ipynb

Import:

import numpy as np
import pandas as pd

Create a reproducible random-number generator:

rng = np.random.default_rng(556)

Part 1 — Lists vs. arrays

python_list = [1, 2, 3, 4]
python_list * 2

Now:

x = np.array([1, 2, 3, 4])
x * 2

Questions

  1. What did multiplication mean for the list?

  2. What did it mean for the array?

  3. Which interpretation is more useful for numerical computing?

Part 2 — Inspect array structure

x = np.array([10, 20, 30])

Inspect:

x.shape
x.ndim
x.dtype
x.size

Explain each output.

Reference answer

List * 2 repeats the list; NumPy array * 2 multiplies each element by 2. The element-wise interpretation is generally what numerical computing needs.

Part 3 — Vector shapes

x = np.array([1, 2, 3])

Create:

x_col = x.reshape(-1, 1)
x_row = x.reshape(1, -1)

Inspect:

x.shape
x_col.shape
x_row.shape

Explain why these shapes are different.

Part 4 — Two-dimensional arrays

X = np.array([
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9],
    [10, 11, 12]
])

Inspect:

X.shape
X.ndim
X.size

Select:

X[0, :]
X[:, 0]
X[1:3, :]

Connect this to row/column selection in pandas.

Reference answer

A vector may have shape (n,), a row matrix (1, n), and a column matrix (n, 1). These are mathematically related but have different NumPy dimensions.

Part 5 — Element-wise operations

Run:

X + 10
X * 2
X ** 2
np.sqrt(X)

What does “element-wise” mean?

Part 6 — Loop vs. vectorized computation

values = list(range(1, 100001))
array = np.array(values)

Loop-style:

squared_loop = [
    value ** 2
    for value in values
]

Vectorized:

squared_vectorized = array ** 2

Verify:

np.array_equal(
    np.array(squared_loop),
    squared_vectorized
)

Reference answer

Element-wise means the operation is applied independently to corresponding entries. NumPy vectorization performs this in compiled routines rather than explicit Python loops.

Part 7 — Benchmark

In Jupyter:

%timeit [value ** 2 for value in values]

and:

%timeit array ** 2

Reflection

Which is faster, and why?

Part 8 — Statistical vectorization

scores = rng.normal(
    80,
    10,
    size=100
)

Calculate:

mean = scores.mean()
sd = scores.std()

z = (
    scores - mean
) / sd

Verify:

z.mean()
z.std()

Where is the iteration?

Reference answer

The vectorized version is usually faster because it avoids Python’s per-iteration overhead and uses optimized compiled code.

Part 9 — Broadcasting a scalar

X = np.array([
    [1, 2, 3],
    [4, 5, 6]
])

X + 10

Write the conceptual matrix that would produce the same result.

Part 10 — Broadcasting across rows

v = np.array([10, 20, 30])

X + v

Inspect:

X.shape
v.shape

Why are the shapes compatible?

Reference answer

Broadcasting can subtract a compatible vector from each row without explicitly copying that vector into a full matrix.

Part 11 — Broadcasting failure

bad = np.array([10, 20])

Try:

X + bad

Read the error and explain why the shapes are incompatible.

Reference answer

Shapes are broadcast-compatible when aligned dimensions are equal or one of them is 1.

Part 12 — Fix the broadcasting problem

bad_col = bad.reshape(-1, 1)

X + bad_col

Inspect:

bad_col.shape

How did reshaping change the meaning of the operation?

Reference answer

The shapes fail because an aligned dimension is neither equal nor 1. Reshape the vector or choose the intended axis so dimensions align.

Part 13 — Understanding axis

X = np.array([
    [10, 20],
    [12, 25],
    [14, 30]
])

Predict the shape, then calculate:

X.mean(axis=0)
X.mean(axis=1)

Reference answer

axis=0 reduces down rows and returns one value per column; axis=1 reduces across columns and returns one value per row.

Part 14 — Column standardization

means = X.mean(axis=0)
sds = X.std(axis=0)

Z = (
    X - means
) / sds

Verify:

Z.mean(axis=0)
Z.std(axis=0)

Which part depends on broadcasting?

Part 15 — keepdims

Compare:

means_a = X.mean(axis=0)

with:

means_b = X.mean(
    axis=0,
    keepdims=True
)

Inspect shapes.

Reference answer

keepdims=True leaves reduced axes with length 1, which helps later broadcasting against the original array.

Part 16 — Element-wise vs. matrix multiplication

A = np.array([
    [1, 2],
    [3, 4]
])

B = np.array([
    [5, 6],
    [7, 8]
])

Compare:

A * B
A @ B

Write out the first element of each result manually.

Reference answer

Element-wise multiplication multiplies corresponding entries. Matrix multiplication forms dot products of rows and columns.

Part 17 — Shape rules for matrix multiplication

A = rng.normal(size=(5, 3))
B = rng.normal(size=(3, 2))

Predict:

(A @ B).shape

Then create an incompatible matrix and explain the resulting error.

Reference answer

For (m,n) @ (n,p), the inner dimensions must match and the output shape is (m,p). The first output element is the dot product of the first row and first column.

Part 18 — Dot products

x = np.array([1, 2, 3])
y = np.array([4, 5, 6])

x @ y

Verify manually.

Reference answer

A dot product multiplies corresponding vector entries and sums the products, yielding a scalar.

Part 19 — Model predictions

X = np.array([
    [1.0, 2.0],
    [2.0, 1.0],
    [3.0, 4.0],
    [4.0, 3.0]
])

beta = np.array([2.0, -1.0])
beta_0 = 5.0

y_hat = beta_0 + X @ beta

If the leading space causes an indentation error, fix it.

Identify where broadcasting occurs.

Part 20 — Solve a linear system

A = np.array([
    [3.0, 1.0],
    [1.0, 2.0]
])

b = np.array([9.0, 8.0])

x = np.linalg.solve(A, b)

Verify:

A @ x

Part 21 — Avoid explicit inversion

x_inverse = np.linalg.inv(A) @ b

Compare with x.

Why might solve() still be preferable?

Part 22 — Build a regression design matrix

n = 100

x1 = rng.normal(size=n)
x2 = rng.normal(size=n)

X = np.column_stack([
    np.ones(n),
    x1,
    x2
])

Inspect:

X.shape
X[:5]

Explain the first column.

Reference answer

The first design-matrix column is all ones and represents the regression intercept; later columns contain predictors.

Part 23 — Simulate a response

beta_true = np.array([
    5.0,
    2.0,
    -1.5
])

y = (
    X @ beta_true
    + rng.normal(0, 1, size=n)
)

Fix any accidental indentation if necessary.

Which operation creates the systematic component of the model?

Reference answer

The simulated response equals the linear predictor plus random noise, so fitted coefficients should be close to but not exactly the simulation parameters.

Part 24 — Least squares

beta_hat, residuals, rank, singular_values = (
    np.linalg.lstsq(
        X,
        y,
        rcond=None
    )
)

Compare:

beta_hat
beta_true

Why are they close but not identical?

Reference answer

np.linalg.lstsq() finds coefficients minimizing squared residuals. Sampling noise prevents exact recovery of the simulation coefficients.

Part 25 — QR decomposition

Q, R = np.linalg.qr(
    X,
    mode="reduced"
)

Inspect shapes and verify:

np.allclose(
    X,
    Q @ R
)

Reference answer

QR factorization writes X = QR; the orthonormal/triangular structure provides a numerically stable route to least squares.

Part 26 — Verify orthogonality

Q.T @ Q

Then:

np.allclose(
    Q.T @ Q,
    np.eye(Q.shape[1])
)

What mathematical property are you checking?

Part 27 — Least squares using QR

beta_qr = np.linalg.solve(
    R,
    Q.T @ y
)

Compare:

beta_qr
beta_hat

Check:

np.allclose(
    beta_qr,
    beta_hat
)

Reference answer

The QR solution should closely match np.linalg.lstsq(). Compare with a tolerance because floating-point paths can differ slightly.

Part 28 — SVD

M = rng.normal(size=(6, 4))

U, s, Vt = np.linalg.svd(
    M,
    full_matrices=False
)

Inspect:

U.shape
s.shape
Vt.shape
s

Part 29 — Reconstruct from SVD

S = np.diag(s)

M_reconstructed = (
    U
    @ S
    @ Vt
)

Check:

np.allclose(
    M,
    M_reconstructed
)

Why use allclose() instead of ==?

Part 30 — Low-rank approximation

k = 2

M_k = (
    U[:, :k]
    @ np.diag(s[:k])
    @ Vt[:k, :]
)

Calculate error:

error = np.linalg.norm(
    M - M_k
)

error

Reference answer

Use np.allclose() rather than exact == because floating-point calculations can differ by tiny rounding amounts.

Part 31 — Approximation error vs. rank

Repeat for:

k = 1, 2, 3, 4

Store the errors.

What happens as more singular components are retained?

Part 32 — SVD on correlated data

x1 = rng.normal(size=200)

x2 = (
    2 * x1
    + rng.normal(
        0,
        0.1,
        size=200
    )
)

x3 = rng.normal(size=200)

X_corr = np.column_stack([
    x1,
    x2,
    x3
])

X_centered = (
    X_corr
    - X_corr.mean(axis=0)
)

Compute:

U, s, Vt = np.linalg.svd(
    X_centered,
    full_matrices=False
)

s

Does one direction appear weaker than the others?

Reference answer

With correlated columns, SVD typically shows one or more dominant singular values and smaller values corresponding to redundant directions.

Part 33 — Condition number

np.linalg.cond(X_centered)

Compare with:

X_independent = rng.normal(
    size=(200, 3)
)

np.linalg.cond(X_independent)

How does this connect to multicollinearity?

Reference answer

The condition number measures numerical sensitivity. A large value indicates an ill-conditioned problem where small perturbations can greatly affect a solution.

Part 34 — Matrix rank

np.linalg.matrix_rank(X_centered)

Now:

X_exact = np.column_stack([
    x1,
    2 * x1,
    x3
])

np.linalg.matrix_rank(X_exact)

Explain the result.

Part 35 — Vectorized distance challenge

points = rng.normal(
    size=(1000, 2)
)

target = np.array([
    0.5,
    -0.5
])

Compute all Euclidean distances without a Python loop:

distances = np.sqrt(
    (
        (points - target) ** 2
    ).sum(axis=1)
)

Identify where broadcasting, element-wise computation, and reduction occur.

Reference answer

Vectorized pairwise distances replace nested Python loops with broadcasting/array operations. Check that the distance matrix is symmetric and has zeros on the diagonal.

Part 36 — Loop version

Implement the same distance calculation with a Python loop.

Verify:

np.allclose(
    distances_loop,
    distances
)

Benchmark both versions.

Part 37 — Mini statistical-computing challenge

Simulate:

n = 500

X = rng.normal(
    size=(n, 5)
)

beta = np.array([
    2.0,
    -1.0,
    0.5,
    3.0,
    -2.0
])

y = (
    4
    + X @ beta
    + rng.normal(
        0,
        2,
        size=n
    )
)

Fix any accidental indentation before execution.

Complete:

  1. Standardize all predictor columns using broadcasting.

  2. Create a design matrix with an intercept.

  3. Estimate coefficients using np.linalg.lstsq().

  4. Compute fitted values.

  5. Compute residuals.

  6. Compute residual sum of squares.

  7. Compute a QR decomposition.

  8. Verify QR reconstruction.

  9. Compute singular values.

  10. Compute the condition number.

Reference answer

For the challenge, express the calculation with arrays, predict the shapes first, and use vectorization/broadcasting only when it preserves the intended mathematics.

Part 38 — Compare implementations

For one computation, implement:

Version A

Explicit Python loop.

Version B

Vectorized NumPy expression.

Compare:

  • clarity;

  • number of lines;

  • speed;

  • correspondence with mathematical notation.

Reference answer

Equivalent implementations should match numerically within tolerance. Verify correctness first, then compare clarity, memory use, and speed.

Part 39 — Numerical validation

For every decomposition, ask:

Shape

Do the matrix dimensions make sense?

Reconstruction

Does Q @ R recover X?

Does U @ S @ Vt recover the matrix?

Orthogonality

Does Q.T @ Q approximate the identity matrix?

Tolerance

Are numerical comparisons performed with:

np.allclose()

rather than exact equality?

Part 40 — Git checkpoint

git status
git add .
git commit -m "Complete Week 9 matrix computation exercises"
git push

Part 41 — Final reflection

1. Arrays

Why are NumPy arrays better suited than Python lists for many numerical computations?

2. Shape

Why is (n,) different from (n, 1)?

3. Vectorization

What does vectorization mean?

4. Broadcasting

Explain broadcasting in your own words.

5. Axis

What does axis=0 mean when calculating a matrix mean?

6. Multiplication

What is the difference between A * B and A @ B?

7. Linear systems

Why prefer solve() to explicitly calculating a matrix inverse?

8. QR

What are the main properties of Q and R?

9. SVD

What do singular values tell us?

10. Statistics

How can SVD or conditioning reveal problems such as multicollinearity?

Completion checklist

  • Created Week 9 notebook

  • Created NumPy arrays

  • Inspected shape, dimension, size, and dtype

  • Compared (n,), (n,1), and (1,n)

  • Used element-wise operations

  • Compared loops and vectorization

  • Benchmarked vectorized operations

  • Standardized data using vectorization

  • Applied scalar broadcasting

  • Applied row/column broadcasting

  • Diagnosed a broadcasting error

  • Used axis=0 and axis=1

  • Used keepdims=True

  • Distinguished element-wise from matrix multiplication

  • Used matrix multiplication shape rules

  • Computed dot products

  • Generated model predictions with matrix algebra

  • Solved a linear system

  • Compared solving with explicit inversion

  • Built a regression design matrix

  • Estimated least squares with lstsq

  • Computed QR decomposition

  • Verified QR reconstruction

  • Verified orthogonality

  • Solved least squares using QR

  • Computed SVD

  • Reconstructed a matrix from SVD

  • Built a low-rank approximation

  • Measured approximation error

  • Explored correlated predictors

  • Calculated a condition number

  • Calculated matrix rank

  • Implemented vectorized distances

  • Compared loop and vectorized implementations

  • Completed the statistical-computing challenge

  • Used np.allclose() for numerical checks

  • Committed Week 9 work to Git

  • Pushed work to GitHub

What you should now understand

Statistical problem
       ↓
Array representation
       ↓
Inspect shapes
       ↓
Vectorize
       ↓
Broadcast
       ↓
Matrix multiplication
       ↓
Numerical linear algebra
       ↓
QR / SVD
       ↓
Stable and efficient computation

Next week we will use this numerical foundation for random-number generation, inverse-transform sampling, and Monte Carlo simulation.

Reference answer

Core ideas are shape reasoning, broadcasting, axis semantics, vectorization, matrix multiplication, stable linear algebra, and tolerance-based numerical comparisons.

9. Reference Activity Notebooks
Activity 4.2
9. Reference Activity Notebooks
Activity 5.2