01: Matrix Derivatives

Why this matters

Throughout QTM 447, training a model means solving an optimization problem:

\min_\theta L(\theta)

where the parameters \theta are vectors or matrices, not scalars.

In principle, you could always: - expand every expression into scalar sums, - compute partial derivatives \partial L / \partial \theta_{ij} one at a time, - and then reassemble the result.

In practice, this approach is slow, error-prone, and impossible to scale.

Matrix calculus provides shortcut rules that let us differentiate entire vectors and matrices at once. Importantly:

These rules are not new mathematics — they are exactly what you would obtain by doing the scalar derivatives by hand and collecting terms.

In this course, we will use a small, focused set of matrix derivative rules that cover nearly everything needed for: - logistic regression - multinomial (softmax) regression - linear layers in neural networks - backpropagation

You should think of these rules as compression tools for derivatives you already know how to compute.

A guiding principle

Shape discipline

If L(\theta) is a scalar loss and \theta \in \mathbb{R}^d, then
\frac{\partial L}{\partial \theta} \in \mathbb{R}^d.

If \theta is a matrix, the gradient has the same shape as \theta.

If your derivative has the wrong shape, something is wrong — no exceptions.

This simple check will catch most errors in optimization and backpropagation.

Matrix Derivative Rules

We now introduce the core rules that will be reused throughout the optimization module. Each rule is shown in:

matrix form (the shortcut),
index/sum notation (what you would get by hand),
reassembled matrix form.

Rule 1: Linear Forms

Let
\mathbf{a} \in \mathbb{R}^d be constant and
\mathbf{x} \in \mathbb{R}^d be variable.

\frac{\partial}{\partial \mathbf{x}} (\mathbf{a}^T \mathbf{x}) = \mathbf{a}

This rule appears constantly: any time a loss depends linearly on a parameter vector, the gradient simply returns the coefficient vector.

Derivation via index notation

Matrix form f(\mathbf{x}) = \mathbf{a}^T \mathbf{x}
Expand into a sum f(\mathbf{x}) = \sum_{i=1}^d a_i x_i
Partial derivative \frac{\partial f}{\partial x_k} = a_k
Reassemble \nabla_{\mathbf{x}} f = \begin{bmatrix} a_1 \\ \vdots \\ a_d \end{bmatrix} = \mathbf{a}

There is no trick here — the matrix rule is just a compact way of writing the scalar result.

Rule 2: Quadratic Forms

Let
\mathbf{A} \in \mathbb{R}^{d \times d}.

\frac{\partial}{\partial \mathbf{x}} (\mathbf{x}^T \mathbf{A} \mathbf{x}) = (\mathbf{A} + \mathbf{A}^T)\mathbf{x}

If \mathbf{A} is symmetric: \frac{\partial}{\partial \mathbf{x}} (\mathbf{x}^T \mathbf{A} \mathbf{x}) = 2\mathbf{A}\mathbf{x}

Quadratic forms appear whenever we:

study curvature,
compute Hessians,
analyze Newton’s method.

Derivation via index notation

Matrix form f(\mathbf{x}) = \mathbf{x}^T \mathbf{A} \mathbf{x}
Expand f(\mathbf{x}) = \sum_i \sum_j x_i A_{ij} x_j
Differentiate w.r.t. x_k
Terms contribute when i = k or j = k: \frac{\partial f}{\partial x_k} = \sum_j A_{kj} x_j + \sum_i A_{ik} x_i
Recognize matrix products = (\mathbf{A}\mathbf{x})_k + (\mathbf{A}^T\mathbf{x})_k
Reassemble \nabla_{\mathbf{x}} f = (\mathbf{A} + \mathbf{A}^T)\mathbf{x}

Rule 3: Least Squares / Norms

Let
\mathbf{A} \in \mathbb{R}^{n \times d},
\mathbf{b} \in \mathbb{R}^n.

\frac{\partial}{\partial \mathbf{x}} \|\mathbf{A}\mathbf{x} - \mathbf{b}\|^2 = 2\mathbf{A}^T(\mathbf{A}\mathbf{x} - \mathbf{b})

This is one of the most important gradients in machine learning.

It already has the structure we will see repeatedly: \text{gradient} = (\text{inputs})^T \times (\text{residuals})

Later:

in logistic regression, the residual becomes p - y
in multinomial regression, it becomes P - Y
in neural networks, this is exactly the linear-layer backprop rule

Rule 4: Matrix–Matrix Derivatives (Trace Trick)

For matrices \mathbf{X} and \mathbf{A} of the same shape:

\frac{\partial}{\partial \mathbf{X}} \operatorname{tr}(\mathbf{A}^T \mathbf{X}) = \mathbf{A}

This is useful because many matrix expressions can be rewritten using the trace, for example:

\|\mathbf{X}\|_F^2 = \operatorname{tr}(\mathbf{X}^T \mathbf{X})

The trace trick lets us extend vector calculus ideas cleanly to matrix-valued parameters.

Numerical Verification

Matrix calculus rules should always be checked numerically when implementing new models.

Below, we verify the quadratic-form gradient using a finite-difference approximation.

Code

import numpy as np

def numerical_gradient(f, x, epsilon=1e-5):
    """
    Compute numerical gradient of f at x using central differences.
    """
    grad = np.zeros_like(x)
    it = np.nditer(x, flags=['multi_index'], op_flags=['readwrite'])
    while not it.finished:
        idx = it.multi_index
        orig_val = x[idx]

        x[idx] = orig_val + epsilon
        f_plus = f(x)

        x[idx] = orig_val - epsilon
        f_minus = f(x)

        grad[idx] = (f_plus - f_minus) / (2 * epsilon)
        x[idx] = orig_val
        it.iternext()
    return grad

np.random.seed(42)
d = 5
A = np.random.randn(d, d)
x = np.random.randn(d)

def f_quad(x):
    return x.T @ A @ x

analytic_grad = (A + A.T) @ x
numeric_grad = numerical_gradient(f_quad, x)

print("Quadratic Form Gradient Check:")
print(f"Analytic: {analytic_grad[:3]}...")
print(f"Numeric:  {numeric_grad[:3]}...")
print(f"Difference: {np.linalg.norm(analytic_grad - numeric_grad):.2e}")

Quadratic Form Gradient Check:
Analytic: [-0.32822323 -2.76535917  1.29888003]...
Numeric:  [-0.32822323 -2.76535917  1.29888003]...
Difference: 3.28e-11

Common mistake: missing transposes

If \mathbf{A} is not symmetric, the gradient is

(\mathbf{A} + \mathbf{A}^T)\mathbf{x},

not 2\mathbf{A}\mathbf{x}.

This mistake is easy to make when memorizing formulas, and it silently breaks:

Hessian calculations,
Newton’s method,
curvature-based reasoning later in the course.

Key takeaways

Matrix derivative rules are compressed versions of scalar derivatives, not new mathematics.
Most gradients used in machine learning reduce to the pattern \text{gradient} = (\text{inputs})^T \times (\text{residuals}).
Shape checks are your first line of defense against implementation bugs.
Always verify new analytic gradients with a numerical check before trusting them.

Self-check questions

What is
\frac{\partial}{\partial \mathbf{x}} (\mathbf{x}^T \mathbf{x})?
If \mathbf{A} is not symmetric, why does
\mathbf{x}^T \mathbf{A} \mathbf{x} depend on both \mathbf{A} and \mathbf{A}^T?
Using the rules above, derive the gradient of
f(\mathbf{w}) = \|\mathbf{X}\mathbf{w} - \mathbf{y}\|^2, and identify which term plays the role of the inputs and which term plays the residuals.

--- title: "01: Matrix Derivatives" format: html: code-fold: true jupyter: python3 --- # Why this matters Throughout QTM 447, training a model means **solving an optimization problem**: $$ \min_\theta L(\theta) $$ where the parameters $\theta$ are **vectors or matrices**, not scalars. In principle, you could always: - expand every expression into scalar sums, - compute partial derivatives $\partial L / \partial \theta_{ij}$ one at a time, - and then reassemble the result. In practice, this approach is slow, error-prone, and impossible to scale. **Matrix calculus provides shortcut rules** that let us differentiate *entire vectors and matrices at once*. Importantly: > These rules are not new mathematics — they are exactly what you would obtain by doing the scalar derivatives by hand and collecting terms. In this course, we will use a **small, focused set of matrix derivative rules** that cover nearly everything needed for: - logistic regression - multinomial (softmax) regression - linear layers in neural networks - backpropagation You should think of these rules as **compression tools** for derivatives you already know how to compute. --- # A guiding principle ::: {.callout-tip} ## Shape discipline If $L(\theta)$ is a **scalar loss** and $\theta \in \mathbb{R}^d$, then $$ \frac{\partial L}{\partial \theta} \in \mathbb{R}^d. $$ If $\theta$ is a **matrix**, the gradient has the **same shape** as $\theta$. If your derivative has the wrong shape, something is wrong — no exceptions. ::: This simple check will catch most errors in optimization and backpropagation. --- # Matrix Derivative Rules We now introduce the core rules that will be reused throughout the optimization module. Each rule is shown in: 1. matrix form (the shortcut), 2. index/sum notation (what you would get by hand), 3. reassembled matrix form. --- ## Rule 1: Linear Forms Let $\mathbf{a} \in \mathbb{R}^d$ be constant and $\mathbf{x} \in \mathbb{R}^d$ be variable. $$ \frac{\partial}{\partial \mathbf{x}} (\mathbf{a}^T \mathbf{x}) = \mathbf{a} $$ This rule appears constantly: any time a loss depends *linearly* on a parameter vector, the gradient simply returns the coefficient vector. ### Derivation via index notation 1. **Matrix form** $$ f(\mathbf{x}) = \mathbf{a}^T \mathbf{x} $$ 2. **Expand into a sum** $$ f(\mathbf{x}) = \sum_{i=1}^d a_i x_i $$ 3. **Partial derivative** $$ \frac{\partial f}{\partial x_k} = a_k $$ 4. **Reassemble** $$ \nabla_{\mathbf{x}} f = \begin{bmatrix} a_1 \\ \vdots \\ a_d \end{bmatrix} = \mathbf{a} $$ There is no trick here — the matrix rule is just a compact way of writing the scalar result. --- ## Rule 2: Quadratic Forms Let $\mathbf{A} \in \mathbb{R}^{d \times d}$. $$ \frac{\partial}{\partial \mathbf{x}} (\mathbf{x}^T \mathbf{A} \mathbf{x}) = (\mathbf{A} + \mathbf{A}^T)\mathbf{x} $$ If $\mathbf{A}$ is symmetric: $$ \frac{\partial}{\partial \mathbf{x}} (\mathbf{x}^T \mathbf{A} \mathbf{x}) = 2\mathbf{A}\mathbf{x} $$ Quadratic forms appear whenever we: - study curvature, - compute Hessians, - analyze Newton’s method. ### Derivation via index notation 1. **Matrix form** $$ f(\mathbf{x}) = \mathbf{x}^T \mathbf{A} \mathbf{x} $$ 2. **Expand** $$ f(\mathbf{x}) = \sum_i \sum_j x_i A_{ij} x_j $$ 3. **Differentiate w.r.t. $x_k$** Terms contribute when $i = k$ or $j = k$: $$ \frac{\partial f}{\partial x_k} = \sum_j A_{kj} x_j + \sum_i A_{ik} x_i $$ 4. **Recognize matrix products** $$ = (\mathbf{A}\mathbf{x})_k + (\mathbf{A}^T\mathbf{x})_k $$ 5. **Reassemble** $$ \nabla_{\mathbf{x}} f = (\mathbf{A} + \mathbf{A}^T)\mathbf{x} $$ --- ## Rule 3: Least Squares / Norms Let $\mathbf{A} \in \mathbb{R}^{n \times d}$, $\mathbf{b} \in \mathbb{R}^n$. $$ \frac{\partial}{\partial \mathbf{x}} \|\mathbf{A}\mathbf{x} - \mathbf{b}\|^2 = 2\mathbf{A}^T(\mathbf{A}\mathbf{x} - \mathbf{b}) $$ This is one of the **most important gradients in machine learning**. It already has the structure we will see repeatedly: $$ \text{gradient} = (\text{inputs})^T \times (\text{residuals}) $$ Later: - in logistic regression, the residual becomes $p - y$ - in multinomial regression, it becomes $P - Y$ - in neural networks, this is exactly the linear-layer backprop rule --- ## Rule 4: Matrix–Matrix Derivatives (Trace Trick) For matrices $\mathbf{X}$ and $\mathbf{A}$ of the same shape: $$ \frac{\partial}{\partial \mathbf{X}} \operatorname{tr}(\mathbf{A}^T \mathbf{X}) = \mathbf{A} $$ This is useful because many matrix expressions can be rewritten using the trace, for example: $$ \|\mathbf{X}\|_F^2 = \operatorname{tr}(\mathbf{X}^T \mathbf{X}) $$ The trace trick lets us extend vector calculus ideas cleanly to matrix-valued parameters. --- # Numerical Verification Matrix calculus rules should **always** be checked numerically when implementing new models. Below, we verify the quadratic-form gradient using a finite-difference approximation. ```{python} import numpy as np def numerical_gradient(f, x, epsilon=1e-5): """ Compute numerical gradient of f at x using central differences. """ grad = np.zeros_like(x) it = np.nditer(x, flags=['multi_index'], op_flags=['readwrite']) while not it.finished: idx = it.multi_index orig_val = x[idx] x[idx] = orig_val + epsilon f_plus = f(x) x[idx] = orig_val - epsilon f_minus = f(x) grad[idx] = (f_plus - f_minus) / (2 * epsilon) x[idx] = orig_val it.iternext() return grad np.random.seed(42) d = 5 A = np.random.randn(d, d) x = np.random.randn(d) def f_quad(x): return x.T @ A @ x analytic_grad = (A + A.T) @ x numeric_grad = numerical_gradient(f_quad, x) print("Quadratic Form Gradient Check:") print(f"Analytic: {analytic_grad[:3]}...") print(f"Numeric: {numeric_grad[:3]}...") print(f"Difference: {np.linalg.norm(analytic_grad - numeric_grad):.2e}") ``` ::: {.callout-note} ## Common mistake: missing transposes If $\mathbf{A}$ is **not symmetric**, the gradient is $$ (\mathbf{A} + \mathbf{A}^T)\mathbf{x}, $$ **not** $2\mathbf{A}\mathbf{x}$. This mistake is easy to make when memorizing formulas, and it silently breaks: - Hessian calculations, - Newton’s method, - curvature-based reasoning later in the course. ::: --- # Key takeaways - Matrix derivative rules are **compressed versions of scalar derivatives**, not new mathematics. - Most gradients used in machine learning reduce to the pattern $$ \text{gradient} = (\text{inputs})^T \times (\text{residuals}). $$ - Shape checks are your first line of defense against implementation bugs. - Always verify new analytic gradients with a numerical check before trusting them. --- # Self-check questions 1. What is $$ \frac{\partial}{\partial \mathbf{x}} (\mathbf{x}^T \mathbf{x})? $$ 2. If $\mathbf{A}$ is not symmetric, why does $\mathbf{x}^T \mathbf{A} \mathbf{x}$ depend on both $\mathbf{A}$ and $\mathbf{A}^T$? 3. Using the rules above, derive the gradient of $$ f(\mathbf{w}) = \|\mathbf{X}\mathbf{w} - \mathbf{y}\|^2, $$ and identify which term plays the role of the **inputs** and which term plays the **residuals**.