05: Multinomial Logistic Derivatives

Why this matters

Multinomial logistic regression (softmax regression) is the first place where parameter tensors appear naturally.

The parameters are no longer vectors, but matrices. The gradients are no longer vectors, but matrices built from outer products.

This is not a new idea — it is exactly the same backpropagation logic we used in the binary case, now applied to a higher-dimensional computational graph.

This section is the bridge from classical statistics to neural network layers.

Model setup (single observation)

For a single observation \mathbf x \in \mathbb R^d with K possible classes (note that the intercept will be handled by adding a 1 to the front of the feature vector):

Parameters
\mathbf W \in \mathbb R^{d \times K} \quad\text{(one column per class)}
Logits
\mathbf z = \mathbf W^\top \mathbf x \in \mathbb R^K
Probabilities (softmax)
p_k = \frac{e^{z_k}}{\sum_{j=1}^K e^{z_j}}, \qquad \mathbf p \in \mathbb R^K
Targets
\mathbf y \in \mathbb R^K \quad\text{(one-hot encoded)}
Loss
L_i = -\sum_{k=1}^K y_k \log p_k

The computation pipeline is:

\mathbf x \;\longrightarrow\; \mathbf z = \mathbf W^\top \mathbf x \;\longrightarrow\; \mathbf p = \text{softmax}(\mathbf z) \;\longrightarrow\; L_i(\mathbf p, \mathbf y).

Computational graph perspective

This graph mirrors the binary case:

a linear node producing logits
a nonlinear normalization (softmax)
a loss node

The only difference is dimensionality.

Backpropagation proceeds exactly the same way: we compute local derivatives at each node and multiply them using the chain rule.

Gradient derivation via the chain rule

We now walk backward through the graph for a single observation.

Step 1: Loss derivative with respect to logits

Loss derivative (output node)

For a single observation, the cross-entropy loss is

L_i = -\sum_{k=1}^K y_k \log p_k,

where \mathbf y is one-hot encoded and \mathbf p = \text{softmax}(\mathbf z).

Differentiating with respect to the probability vector \mathbf p gives

\frac{\partial L_i}{\partial p_k} = -\frac{y_k}{p_k}, \qquad k = 1,\dots,K.

In vector form:

\frac{\partial L_i}{\partial \mathbf p} = -\frac{\mathbf y}{\mathbf p},

where the division is understood elementwise.

Shape check

\frac{\partial L_i}{\partial \mathbf p} \in \mathbb R^{K}.

Softmax derivative (activation node)

The softmax function is

p_k = \frac{e^{z_k}}{\sum_{j=1}^K e^{z_j}}.

Its Jacobian is a K \times K matrix with entries

\frac{\partial p_k}{\partial z_j} = p_k(\delta_{kj} - p_j),

where \delta_{kj} is the Kronecker delta.

In matrix form, the Jacobian is

\frac{\partial \mathbf p}{\partial \mathbf z} = \mathbf J_{\text{softmax}} = \operatorname{diag}(\mathbf p) - \mathbf p \mathbf p^\top.

Chain rule at the logits

Applying the chain rule:

\frac{\partial L_i}{\partial \mathbf z} = \left( \frac{\partial \mathbf p}{\partial \mathbf z} \right)^\top \frac{\partial L_i}{\partial \mathbf p}.

Substituting the expressions above:

\frac{\partial L_i}{\partial \mathbf z} = \left(\operatorname{diag}(\mathbf p) - \mathbf p \mathbf p^\top\right) \left(-\frac{\mathbf y}{\mathbf p}\right).

Carrying out the multiplication (using the fact that \mathbf y is one-hot encoded) yields the remarkable simplification:

\frac{\partial L_i}{\partial \mathbf z} = \mathbf p - \mathbf y.

Shape check

\mathbf p - \mathbf y \in \mathbb R^{K}.

Interpretation

As in binary logistic regression, the final gradient depends only on a residual:

\text{prediction} - \text{truth}.

The complexity of the softmax Jacobian collapses once it is paired with the cross-entropy loss. This is why softmax and cross-entropy are almost always used together in practice.

Step 2: Logits with respect to the weight matrix

The logits are defined by

\mathbf z = \mathbf W^\top \mathbf x.

Each logit is

z_k = \mathbf w_k^\top \mathbf x,

where \mathbf w_k is the k-th column of \mathbf W.

The derivative of \mathbf z with respect to \mathbf W produces an outer product:

\frac{\partial L_i}{\partial \mathbf W} = \mathbf x (\mathbf p - \mathbf y)^\top.

Why an outer product shows up here

At first glance, the gradient \frac{\partial L_i}{\partial \mathbf W} = \mathbf x(\mathbf p - \mathbf y)^\top can feel like it came out of nowhere. It didn’t. It is the most natural consequence of:

each class k having its own weight vector \mathbf w_k, and
the linear node producing logits by dot products with the same input \mathbf x.

Step A: Think column-by-column

Recall that \mathbf W is a matrix with one column per class: \mathbf W = [\mathbf w_1 \;\; \mathbf w_2 \;\; \cdots \;\; \mathbf w_K], \qquad \mathbf w_k \in \mathbb R^d.

The logits are: \mathbf z = \mathbf W^\top \mathbf x \quad\Longrightarrow\quad z_k = \mathbf w_k^\top \mathbf x.

So the loss depends on \mathbf W through the K separate dot products \mathbf w_k^\top \mathbf x.

Step B: Apply the chain rule to one column

Backpropagation tells us: \frac{\partial L_i}{\partial \mathbf w_k} = \frac{\partial L_i}{\partial z_k} \cdot \frac{\partial z_k}{\partial \mathbf w_k}.

But z_k = \mathbf w_k^\top \mathbf x \quad\Rightarrow\quad \frac{\partial z_k}{\partial \mathbf w_k} = \mathbf x.

So: \frac{\partial L_i}{\partial \mathbf w_k} = \left(\frac{\partial L_i}{\partial z_k}\right)\mathbf x.

Now use the standard softmax + cross-entropy result: \frac{\partial L_i}{\partial \mathbf z} = \mathbf p - \mathbf y \quad\Rightarrow\quad \frac{\partial L_i}{\partial z_k} = (p_k - y_k).

Therefore: \frac{\partial L_i}{\partial \mathbf w_k} = (p_k - y_k)\mathbf x.

Interpretation: every class’s weight vector gets updated in the direction of \mathbf x, scaled by its own residual (p_k - y_k).

Step C: Stack the K column gradients into a matrix

If the gradient with respect to each column is: \frac{\partial L_i}{\partial \mathbf w_k} = (p_k - y_k)\mathbf x, then the full gradient matrix is: \frac{\partial L_i}{\partial \mathbf W} = \big[\,(p_1-y_1)\mathbf x \;\; (p_2-y_2)\mathbf x \;\; \cdots \;\; (p_K-y_K)\mathbf x\,\big].

Factor out \mathbf x:

\frac{\partial L_i}{\partial \mathbf W} = \mathbf x \big[(p_1-y_1)\;\; (p_2-y_2)\;\; \cdots\;\; (p_K-y_K)\big].

But that bracketed row vector is exactly (\mathbf p - \mathbf y)^\top.

So we get: \frac{\partial L_i}{\partial \mathbf W} = \mathbf x(\mathbf p - \mathbf y)^\top.

Why this is called an outer product

\mathbf x is a column vector in \mathbb R^{d \times 1}
(\mathbf p - \mathbf y)^\top is a row vector in \mathbb R^{1 \times K}

Multiplying them produces a matrix: (d \times 1)\cdot(1 \times K) = (d \times K), whose (j,k) entry is: \left[\mathbf x(\mathbf p - \mathbf y)^\top\right]_{jk} = x_j (p_k - y_k).

So each weight W_{jk} (feature j, class k) gets the update signal:

feature value \times class residual.

This is the key pattern in neural networks

This same outer-product structure reappears in the backward pass of a linear layer:

input activations \;\times\; upstream error signal

Multinomial logistic regression is the first place you can see it in a clean, closed-form setting.

Shape check

\mathbf x: (d \times 1)
(\mathbf p - \mathbf y)^\top: (1 \times K)
Result: (d \times K), matching \mathbf W.

This is the key tensor-gradient pattern that appears throughout deep learning.

Batch gradient (matrix form)

For a dataset of n observations:

\mathbf X \in \mathbb R^{n \times d}
\mathbf P, \mathbf Y \in \mathbb R^{n \times K}

Summing per-observation gradients yields:

\nabla_{\mathbf W} L = \sum_{i=1}^n \mathbf x_i (\mathbf p_i - \mathbf y_i)^\top = \mathbf X^\top (\mathbf P - \mathbf Y).

Shape check

[d \times n] \cdot [n \times K] = [d \times K].

This is the exact gradient used in practice.

Tensor gradient interpretation

The multinomial gradient can be viewed as a sum of outer products:

\nabla_{\mathbf W} L = \sum_{i=1}^n \mathbf x_i \otimes (\mathbf p_i - \mathbf y_i).

Each observation contributes a rank-1 update.

This perspective is essential for understanding:

vectorized implementations
GPU efficiency
linear layers in neural networks

Hessian structure (high-level view)

The parameter vector now has dK entries. The Hessian therefore has dimension (dK) \times (dK).

It has a block structure: each block corresponds to interactions between class-specific parameter vectors.

The (j,k) block is

\mathbf H_{jk} = \mathbf X^\top \mathbf S_{jk} \mathbf X,

where \mathbf S_{jk} is diagonal with entries

p_{ij}(\delta_{jk} - p_{ik}).

This coupling across classes is what makes the Hessian expensive to store and invert.

Scaling issue

For d = 100 and K = 10, the Hessian has size 1000 \times 1000.

Inverting it costs O((dK)^3) operations.

This is why full Newton methods are rarely used for multinomial models or neural networks.

Key takeaways

Multinomial logistic regression uses the same computational graph logic as the binary case
Gradients become matrices built from outer products
The batch gradient has the clean form
\nabla_{\mathbf W} L = \mathbf X^\top (\mathbf P - \mathbf Y)
This is exactly the backward pass of a linear layer followed by softmax
Hessians scale poorly, motivating first-order methods in deep learning

Self-check questions

What is the shape of \mathbf W if d=10 and K=3?
Why does \sum_k p_k = 1 imply \sum_k \frac{\partial L}{\partial z_k} = 0?
How does this gradient generalize the binary logistic regression gradient?

--- title: "05: Multinomial Logistic Derivatives" format: html: code-fold: true jupyter: python3 --- # Why this matters Multinomial logistic regression (softmax regression) is the first place where **parameter tensors** appear naturally. The parameters are no longer vectors, but matrices. The gradients are no longer vectors, but matrices built from **outer products**. This is not a new idea — it is exactly the same backpropagation logic we used in the binary case, now applied to a higher-dimensional computational graph. > This section is the bridge from classical statistics to neural network layers. # Model setup (single observation) For a single observation $\mathbf x \in \mathbb R^d$ with $K$ possible classes (note that the intercept will be handled by adding a 1 to the front of the feature vector): - **Parameters** $$ \mathbf W \in \mathbb R^{d \times K} \quad\text{(one column per class)} $$ - **Logits** $$ \mathbf z = \mathbf W^\top \mathbf x \in \mathbb R^K $$ - **Probabilities (softmax)** $$ p_k = \frac{e^{z_k}}{\sum_{j=1}^K e^{z_j}}, \qquad \mathbf p \in \mathbb R^K $$ - **Targets** $$ \mathbf y \in \mathbb R^K \quad\text{(one-hot encoded)} $$ - **Loss** $$ L_i = -\sum_{k=1}^K y_k \log p_k $$ The computation pipeline is: $$ \mathbf x \;\longrightarrow\; \mathbf z = \mathbf W^\top \mathbf x \;\longrightarrow\; \mathbf p = \text{softmax}(\mathbf z) \;\longrightarrow\; L_i(\mathbf p, \mathbf y). $$ # Computational graph perspective This graph mirrors the binary case: - a **linear node** producing logits - a **nonlinear normalization** (softmax) - a **loss node** The only difference is dimensionality. Backpropagation proceeds exactly the same way: we compute local derivatives at each node and multiply them using the chain rule. # Gradient derivation via the chain rule We now walk backward through the graph for a single observation. --- ## Step 1: Loss derivative with respect to logits ### Loss derivative (output node) For a single observation, the cross-entropy loss is $$ L_i = -\sum_{k=1}^K y_k \log p_k, $$ where $\mathbf y$ is one-hot encoded and $\mathbf p = \text{softmax}(\mathbf z)$. Differentiating with respect to the probability vector $\mathbf p$ gives $$ \frac{\partial L_i}{\partial p_k} = -\frac{y_k}{p_k}, \qquad k = 1,\dots,K. $$ In vector form: $$ \frac{\partial L_i}{\partial \mathbf p} = -\frac{\mathbf y}{\mathbf p}, $$ where the division is understood elementwise. **Shape check** $$ \frac{\partial L_i}{\partial \mathbf p} \in \mathbb R^{K}. $$ ### Softmax derivative (activation node) The softmax function is $$ p_k = \frac{e^{z_k}}{\sum_{j=1}^K e^{z_j}}. $$ Its Jacobian is a $K \times K$ matrix with entries $$ \frac{\partial p_k}{\partial z_j} = p_k(\delta_{kj} - p_j), $$ where $\delta_{kj}$ is the Kronecker delta. In matrix form, the Jacobian is $$ \frac{\partial \mathbf p}{\partial \mathbf z} = \mathbf J_{\text{softmax}} = \operatorname{diag}(\mathbf p) - \mathbf p \mathbf p^\top. $$ ### Chain rule at the logits Applying the chain rule: $$ \frac{\partial L_i}{\partial \mathbf z} = \left( \frac{\partial \mathbf p}{\partial \mathbf z} \right)^\top \frac{\partial L_i}{\partial \mathbf p}. $$ Substituting the expressions above: $$ \frac{\partial L_i}{\partial \mathbf z} = \left(\operatorname{diag}(\mathbf p) - \mathbf p \mathbf p^\top\right) \left(-\frac{\mathbf y}{\mathbf p}\right). $$ Carrying out the multiplication (using the fact that $\mathbf y$ is one-hot encoded) yields the remarkable simplification: $$ \frac{\partial L_i}{\partial \mathbf z} = \mathbf p - \mathbf y. $$ **Shape check** $$ \mathbf p - \mathbf y \in \mathbb R^{K}. $$ ### Interpretation As in binary logistic regression, the final gradient depends only on a **residual**: $$ \text{prediction} - \text{truth}. $$ The complexity of the softmax Jacobian collapses once it is paired with the cross-entropy loss. This is why softmax and cross-entropy are almost always used together in practice. ## Step 2: Logits with respect to the weight matrix The logits are defined by $$ \mathbf z = \mathbf W^\top \mathbf x. $$ Each logit is $$ z_k = \mathbf w_k^\top \mathbf x, $$ where $\mathbf w_k$ is the $k$-th column of $\mathbf W$. The derivative of $\mathbf z$ with respect to $\mathbf W$ produces an **outer product**: $$ \frac{\partial L_i}{\partial \mathbf W} = \mathbf x (\mathbf p - \mathbf y)^\top. $$ ::: {.callout-tip} ## Why an outer product shows up here At first glance, the gradient $$ \frac{\partial L_i}{\partial \mathbf W} = \mathbf x(\mathbf p - \mathbf y)^\top $$ can feel like it came out of nowhere. It didn’t. It is the most natural consequence of: 1. each class $k$ having its **own weight vector** $\mathbf w_k$, and 2. the linear node producing logits by **dot products** with the same input $\mathbf x$. --- ### Step A: Think column-by-column Recall that $\mathbf W$ is a matrix with one column per class: $$ \mathbf W = [\mathbf w_1 \;\; \mathbf w_2 \;\; \cdots \;\; \mathbf w_K], \qquad \mathbf w_k \in \mathbb R^d. $$ The logits are: $$ \mathbf z = \mathbf W^\top \mathbf x \quad\Longrightarrow\quad z_k = \mathbf w_k^\top \mathbf x. $$ So the loss depends on $\mathbf W$ through the $K$ separate dot products $\mathbf w_k^\top \mathbf x$. --- ### Step B: Apply the chain rule to one column Backpropagation tells us: $$ \frac{\partial L_i}{\partial \mathbf w_k} = \frac{\partial L_i}{\partial z_k} \cdot \frac{\partial z_k}{\partial \mathbf w_k}. $$ But $$ z_k = \mathbf w_k^\top \mathbf x \quad\Rightarrow\quad \frac{\partial z_k}{\partial \mathbf w_k} = \mathbf x. $$ So: $$ \frac{\partial L_i}{\partial \mathbf w_k} = \left(\frac{\partial L_i}{\partial z_k}\right)\mathbf x. $$ Now use the standard softmax + cross-entropy result: $$ \frac{\partial L_i}{\partial \mathbf z} = \mathbf p - \mathbf y \quad\Rightarrow\quad \frac{\partial L_i}{\partial z_k} = (p_k - y_k). $$ Therefore: $$ \frac{\partial L_i}{\partial \mathbf w_k} = (p_k - y_k)\mathbf x. $$ **Interpretation:** every class’s weight vector gets updated in the direction of $\mathbf x$, scaled by its own residual $(p_k - y_k)$. --- ### Step C: Stack the $K$ column gradients into a matrix If the gradient with respect to each column is: $$ \frac{\partial L_i}{\partial \mathbf w_k} = (p_k - y_k)\mathbf x, $$ then the full gradient matrix is: $$ \frac{\partial L_i}{\partial \mathbf W} = \big[\,(p_1-y_1)\mathbf x \;\; (p_2-y_2)\mathbf x \;\; \cdots \;\; (p_K-y_K)\mathbf x\,\big]. $$ Factor out $\mathbf x$: $$ \frac{\partial L_i}{\partial \mathbf W} = \mathbf x \big[(p_1-y_1)\;\; (p_2-y_2)\;\; \cdots\;\; (p_K-y_K)\big]. $$ But that bracketed row vector is exactly $(\mathbf p - \mathbf y)^\top$. So we get: $$ \frac{\partial L_i}{\partial \mathbf W} = \mathbf x(\mathbf p - \mathbf y)^\top. $$ --- ### Why this is called an outer product - $\mathbf x$ is a **column vector** in $\mathbb R^{d \times 1}$ - $(\mathbf p - \mathbf y)^\top$ is a **row vector** in $\mathbb R^{1 \times K}$ Multiplying them produces a matrix: $$ (d \times 1)\cdot(1 \times K) = (d \times K), $$ whose $(j,k)$ entry is: $$ \left[\mathbf x(\mathbf p - \mathbf y)^\top\right]_{jk} = x_j (p_k - y_k). $$ So each weight $W_{jk}$ (feature $j$, class $k$) gets the update signal: > **feature value** $\times$ **class residual**. --- ### This is the key pattern in neural networks This same outer-product structure reappears in the backward pass of a linear layer: - input activations $\;\times\;$ upstream error signal Multinomial logistic regression is the first place you can see it in a clean, closed-form setting. ::: **Shape check** - $\mathbf x$: $(d \times 1)$ - $(\mathbf p - \mathbf y)^\top$: $(1 \times K)$ - Result: $(d \times K)$, matching $\mathbf W$. This is the key tensor-gradient pattern that appears throughout deep learning. # Batch gradient (matrix form) For a dataset of $n$ observations: - $\mathbf X \in \mathbb R^{n \times d}$ - $\mathbf P, \mathbf Y \in \mathbb R^{n \times K}$ Summing per-observation gradients yields: $$ \nabla_{\mathbf W} L = \sum_{i=1}^n \mathbf x_i (\mathbf p_i - \mathbf y_i)^\top = \mathbf X^\top (\mathbf P - \mathbf Y). $$ **Shape check** $$ [d \times n] \cdot [n \times K] = [d \times K]. $$ This is the exact gradient used in practice. # Tensor gradient interpretation The multinomial gradient can be viewed as a sum of outer products: $$ \nabla_{\mathbf W} L = \sum_{i=1}^n \mathbf x_i \otimes (\mathbf p_i - \mathbf y_i). $$ Each observation contributes a rank-1 update. This perspective is essential for understanding: - vectorized implementations - GPU efficiency - linear layers in neural networks # Hessian structure (high-level view) The parameter vector now has $dK$ entries. The Hessian therefore has dimension $(dK) \times (dK)$. It has a **block structure**: each block corresponds to interactions between class-specific parameter vectors. The $(j,k)$ block is $$ \mathbf H_{jk} = \mathbf X^\top \mathbf S_{jk} \mathbf X, $$ where $\mathbf S_{jk}$ is diagonal with entries $$ p_{ij}(\delta_{jk} - p_{ik}). $$ This coupling across classes is what makes the Hessian expensive to store and invert. ::: {.callout-warning} ## Scaling issue For $d = 100$ and $K = 10$, the Hessian has size $1000 \times 1000$. Inverting it costs $O((dK)^3)$ operations. This is why full Newton methods are rarely used for multinomial models or neural networks. ::: # Key takeaways - Multinomial logistic regression uses the **same computational graph logic** as the binary case - Gradients become matrices built from **outer products** - The batch gradient has the clean form $$ \nabla_{\mathbf W} L = \mathbf X^\top (\mathbf P - \mathbf Y) $$ - This is exactly the backward pass of a linear layer followed by softmax - Hessians scale poorly, motivating first-order methods in deep learning # Self-check questions 1. What is the shape of $\mathbf W$ if $d=10$ and $K=3$? 2. Why does $\sum_k p_k = 1$ imply $\sum_k \frac{\partial L}{\partial z_k} = 0$? 3. How does this gradient generalize the binary logistic regression gradient?