Derive Backpropagation for Matrix-Product Layers

Consider a neural-network block whose output is produced by multiplying a sequence of trainable weight matrices together, then applying the resulting composite matrix to an input. There is no nonlinearity inside this block — it is a pure chain of matrix multiplications.

Let the trainable matrices be $W_1, W_2, \ldots, W_n$. Their cumulative product defines a single effective linear map

$$ C = W_1 W_2 \cdots W_n . $$

Given an input vector or mini-batch $X$, the forward pass is

$$ Z = C,X = W_1 W_2 \cdots W_n, X . $$

There is a scalar loss $\mathcal{L}$, and the upstream gradient

$$ G = \frac{\partial \mathcal{L}}{\partial Z} $$

is supplied by the loss or by later layers (so $G$ has the same shape as $Z$).

Derive the backward pass for this block. Specifically:

1. Express the gradient $\dfrac{\partial \mathcal{L}}{\partial W_j}$ with respect to each individual matrix $W_j$, for every $1 \le j \le n$.
2. Show explicitly how the multivariate chain rule (or the matrix differential) applies to the matrix product — i.e. why the answer has the form it does, not just the final formula.
3. Verify that the resulting gradient $\dfrac{\partial \mathcal{L}}{\partial W_j}$ has the same shape as $W_j$.
4. Describe an efficient implementation that avoids recomputing the same prefix and suffix matrix products for every $j$.

Constraints & Assumptions

• Shapes: let $W_j \in \mathbb{R}^{d_{j-1} \times d_j}$, so the composite $C \in \mathbb{R}^{d_0 \times d_n}$. For a mini-batch, $X \in \mathbb{R}^{d_n \times B}$ and $Z, G \in \mathbb{R}^{d_0 \times B}$; for a single example take $B = 1$.
• The matrices need not be square; adjacent inner dimensions just have to match so the product is defined.
• There is no nonlinearity or bias term inside this block — it is a pure chain of matrix multiplications. Any activation lives in the layers that produce $G$.
• "Efficient" means linear in the number of matrices $n$: the whole backward pass should cost $O(n)$ matrix products, not $O(n^2)$, and should avoid materializing large $d_0 \times d_{j-1}$ prefix matrices when $B$ is small.

Clarifying Questions to Ask

• Is $X$ a single vector or a mini-batch, and is the loss summed or averaged over the batch (which scales every $dW_j$)?
• Should I treat $X$ as fixed input, or is the gradient $\partial\mathcal{L}/\partial X$ also required?
• Are the matrices guaranteed non-singular / square, or fully general rectangular?
• Is there any nonlinearity or bias term between the matrices, or is the block a pure linear chain as stated?
• What's the target metric — closed-form gradient formula, an $O(n)$ algorithm, or working code/pseudocode?
• Do we ever need to materialize the composite $C$ at inference, or only its action $CX$?