- 18 Sep, 2025*
TL;DR In this long and technical blog, I explained:
Forward propagation: How inputs flow through the network layer by layer (using matrix operations) to generate predictions.
Computation Graph: How simple scalar examples help visualize backpropagation and build intuition before scaling up.
Backward propagation: How errors flow backward (again with matrix operations) to compute gradients and update weights.
Backpropagation appears quite straightforward when working with scalars or even simple vectors. However, once we step into the world of matrices, things quickly become more complex and difficult to follow. There are extra details and notations that make it less intuitive. Personally, although I managed to understand this concept while…
- 18 Sep, 2025*
TL;DR In this long and technical blog, I explained:
Forward propagation: How inputs flow through the network layer by layer (using matrix operations) to generate predictions.
Computation Graph: How simple scalar examples help visualize backpropagation and build intuition before scaling up.
Backward propagation: How errors flow backward (again with matrix operations) to compute gradients and update weights.
Backpropagation appears quite straightforward when working with scalars or even simple vectors. However, once we step into the world of matrices, things quickly become more complex and difficult to follow. There are extra details and notations that make it less intuitive. Personally, although I managed to understand this concept while preparing for deep learning exams, I usually lose my intuition a few months later. Reviewing it again requires extra effort to bring everything back together. That’s why I decided to write this technical blog post, both as a reminder for myself and as a guide for anyone trying to understand backpropagation in matrix form.
1 - Forward Propagation
When training a deep learning model, the first step is to compute a linear combination of input features. Given an input vector x∈ℝM with M features and a corresponding weight vector w∈ℝM, we calculate a weighted sum and add a bias term b∈ℝ1 to produce a single output z:
z=w1·x1+w2·x2+…+wM·xM+b This is simply the dot product between w and x with an added bias term:
z=wTx+b To visualize this computation, we can represent it in matrix form:
[z]=[w1w2⋯wM][x1x2⋮xM]+[b] or in image representation:
[z]=[bw1w2⋯wM][1x1x2⋮xM] For notational clarity, we’ll denote the bias term as w0=b and the constant input as x0=1. Under this convention, both vectors have dimension M+1: w∈ℝM+1 and x∈ℝM+1. Our equation becomes:
z=wTx or in matrix form:
[z]=[w0w1w2⋯wM][x0x1x2⋮xM] or in image representation:
1.1 - Adding Non-Linearity
The linear combination alone is insufficient for learning complex patterns. To add non-linearity, we apply an activation function σ(·) to the output z.
a=σ(z)=σ(wTx) Although σ often denotes the Sigmoid function, here it represents a general activation function. For this explanation, we’ll use the ReLU activation function in hidden layers due to its simplicity, computational efficiency, resistance to vanishing gradients, and widespread popularity:
a=σ(z)=max(0,z) we can show this with visual representation:
1.2 - Multi-class Classification
In practice, deep learning models often solve multi-class problems where we need to predict one of C possible classes. This requires the output a to be C-dimensional vector, with each dimension representing the score for a particular class. The predicted class corresponds to the entry with the highest activation value.
To generate C outputs, we need C distinct weight vectors, each of dimension (M+1). To achieve this, we need a separate weight vector for each class. Stacking them gives a weight matrix W∈ℝC×(M+1):
a⏟[C×1]=σ(W⏟[C×(M+1)]x⏟[(M+1)×1]) or in expanded form:
[a1a2⋮aC]=σ([W10W11W12…W1MW20W21W22…W2M⋮⋮⋮⋱⋮WC0WC1WC2…WCM][x0x1x2⋮xM]) or in visual form:
1.3 - Batch Processing
The computations described so far process only a single input example (batch size = 1). In practice, we process multiple inputs in parallel to improve computational efficiency.
If we process batch size of N inputs simultaneously, then our input becomes X∈ℝN×(M+1) (note the uppercase X since we now have a matrix rather than a vector). To handle N outputs while maintaining proper matrix dimensions, we transpose our weight matrix to W∈ℝ(M+1)×C:
A⏟[N×C]=σ(X⏟[N×(M+1)]W⏟[(M+1)×C]) or in matrix form:
[A11A12…A1CA21A22…A2C⋮⋮⋱⋮AN1AN2…ANC]=σ([X10X11X12…X1MX20X21X22…X2M⋮⋮⋮⋱⋮XN0XN1XN2…XNM][W01W02…W0CW11W12…W1CW21W22…W2C⋮⋮⋱⋮WM1WM2…WMC]) When dealing with batches, it can sometimes be hard to picture how forward propagation works. We can visualize the process like this:
1.4 - Multi-layer Networks
The calculations presented so far describe a single layer. Deep learning models stack multiple layers to learn increasingly complex representations. We can denote each layer with a superscript [l] to represent the layer number:
Z[l]⏟[N×Hl]=A[l−1]⏟[N×Hl−1]W[l]⏟[Hl−1×Hl]A[l]⏟[N×Hl]=σ(Z[l]⏟[N×Hl]) Here Hl represents the number of hidden units in layer l. The input to layer l activated output from the previous layer l−1 . For simplicity, we usually define A[0]=X (the input layer).
1.5 - Softmax for Classification
For classification problems, we apply the softmax activation function to the final layer’s output to convert the raw scores into a probability distribution. The softmax function ensures that all outputs sum to 1, allowing us to interpret them as class probabilities.
The softmax function is applied row-wise (for each input example) across the C classes:
Y^ij=softmax(Zij)=eZij∑l=1CeZili=1,2,…,N andj=1,2,…,C Intuitively, for input example i, this computes the probability of class j by normalizing the exponential of its score by the sum of exponentials across all C possible classes.
The matrix representation shows this row-wise operation:
[Y^11Y^12…Y^1CY^21Y^22…Y^2C⋮⋮⋱⋮Y^N1Y^N2…Y^NC]=[softmax(Z11Z12…Z1C)softmax(Z21Z22…Z2C)⋮⋮⋱⋮softmax(ZN1ZN2…ZNC)]=[eZ11∑l=1CeZ1leZ12∑l=1CeZ1l…eZ1C∑l=1CeZ1leZ21∑l=1CeZ2leZ22∑l=1CeZ2l…eZ2C∑l=1CeZ2l⋮⋮⋱⋮eZN1∑l=1CeZNleZN2∑l=1CeZNl…eZNC∑l=1CeZNl] If we denote the last layer as L, then the model’s output for N inputs and C classes can be represented as follows:
1.6 - Cross-Entropy Loss Function
After obtaining the predicted probabilities Y^, we measure the model’s performance by comparing these predictions with the ground truth labels Y∈ℝN×C using a loss function. The ground truth is represented as one-hot encoded vectors, where each input example has exactly one correct class. For multi-class classification, we typically use Cross-Entropy (CE) loss:
ℒ=CE(Y,Y^)=−1N∑i=1N∑j=1CYijlogY^ij This formula computes the element-wise product between the true labels Y and the logarithm of predicted probabilities logY^ then averages across all examples to produce a single scalar loss value.
We can express this using matrix operations with the element-wise (Hadamard) product ⊙:
ℒ⏟[1×1]=CE(Y,Y^)=−1N1T⏟[1×N](Y⊙logY^)⏟[N×C]1⏟[C×1] Expanding this matrix operation:
ℒ=−1N[1112…1N]([Y11Y12…Y1CY21Y22…Y2C⋮⋮⋱⋮YN1YN2…YNC]⊙[logY^11logY^12…logY^1ClogY^21logY^22…logY^2C⋮⋮⋱⋮logY^N1logY^N2…logY^NC])[1112⋮1C] We can simplify this by recognizing that multiplying by vectors of ones simply sums all elements in the matrix. Therefore:
ℒ=−1Nsum([Y11logY^11Y12logY^12…Y1ClogY^1CY21logY^21Y22logY^22…Y2ClogY^2C⋮⋮⋱⋮YN1logY^N1YN2logY^N2…YNClogY^NC]) We can show this summation with following image (excluding scaling factor −1N):
ℒij={−logY^ijifYij=10otherwisefori=1,2,…,N This completes our mathematical framework for forward propagation in deep learning models, from single predictions to batch processing with multi-class classification and loss computation.
2 - A Simple Computation Graph
Before diving lots of matrix gradient calculation, I would like to show a computation graph for forward and backward propagation so we can see the cases we need to be careful when we compute gradients. In this case, I want to only show scalar calculations first instead of thinking about matrices so we can adapt this into matrix backpropagation.
a[0]=xz[1]=a[0]·w[1]a[1]=σ(z[1])z[2]=a[1]·w[2]a[2]=σ(z[2])⋮z[l]=a[l−1]·w[l]a[l]=σ(z[l])⋮z[L−1]=a[L−2]·w[L−1]a[L−1]=σ(z[L−1])z[L]=a[L−1]·w[L]y^=softmax(z[L])ℒ=CE(y^)=−∑i=1Cyilogy^i where [l] and [L] represents arbitrary layer l and last layer L, respectively. σ(·) is activation function (ReLU in this case). In this representation, everything is scalar except w[L] because the output must be multi-class for softmax so I just played in the last layer to show a meaningful example.
We can visualize the computation graph of this network as follows:
2.1 - Chain Rule
Backpropagation relies on the chain rule of calculus to compute gradients efficiently. For a composite function f(g(h(x))), the chain rule states:
∂f∂x=∂f∂g·∂g∂h·∂h∂x In neural networks, we apply this principle to decompose the loss gradient into a product of simpler derivatives based on any parameter.
2.2 - Scalar Backward Propagation
After computing loss, we can start computing gradients with respect to it. Since there are C classes, we need to calculate derivative of each prediction y^i:
∂ℒ∂y^i=∂(−yilogy^i)∂y^i=−yiy^i We can visualize this gradient flow like this:
yi=softmax(zi[L])=ezi[L]∑j=1Cezj[L], the derivative ∂ℒ∂zk[L] depends not only on the numerator but also on the denominator:
∂ℒ∂zk[L]=∑i=1C∂ℒ∂y^i·∂y^i∂zk[L] It might be easier to understand when you see visual representation:
∂ℒ∂zk[L]=∑i=1C−yiy^i·∂y^i∂zk[L] The most confusing part begins here, because we now need to carefully compute the local gradient ∂y^i∂zk[L]. To do this, we apply the quotient rule:
∂y^i∂zk[L]=∂(ezi[L]∑j=1Cezj[L])∂zk[L]=(∂ezi[L]∂zk[L])·∑j=1Cezj[L]−ezi[L]·(∑j=1Cezj[L]∂zk[L])(∑j=1Cezj[L])2 For this derivative, there are two distinct cases to consider:
i=k ori≠k—
Case 1: i=k
- The derivative of numerator: ∂ezi[L]∂zk[L]=ezi[L]
- The derivative of the denominator is: ∂(∑j=1Cezj[L])∂zk[L]=ezk[L]
- Substituting these into the quotient rule: ∂y^i∂zk[L]=zi[L]·∑j=1Czj[L]−zi[L]·zk[L](∑j=1Czj[L])2
- Factorizing into: ∂y^i∂zk[L]=ezi[L]∑j=1Cezj[L]·(1−ezk[L]∑j=1Cezj[L])
- Shortly, this is equivalent to: ∂y^i∂zk[L]=y^i·(1−y^k)
Case 2: i≠k
The derivative of the numerator is zero, because ezi does not depend on zk.
The derivative of the denominator is: ∂(∑j=1Cezj[L])∂zk[L]=ezk[L]
- Thus, we have: ∂y^i∂zk[L]=0·∑j=1Czj[L]−zi[L]·zk[L](∑j=1Czj[L])2
- Simplifying: ∂y^i∂zk[L]=−ezi[L]∑j=1Cezj[L]·ezk[L]∑j=1Cezj[L]
- Shortly, this is equivalent to: ∂y^i∂zk[L]=−y^i·y^k— Finally, we can summarize these two cases as
∂y^i∂zk[L]={y^i·(1−y^k) ifi=k−y^i·y^k ifi≠k In literature, it is pretty common to use the following form:
∂y^i∂zk[L]=y^i·(δik−y^k) where δik is 1 if i=k, 0 otherwise. When we combine the gradient of cross-entropy and local gradient, we have:
∂ℒ∂zk[L]=∑i=1C−yiy^i·y^i·(δik−y^k) where y^i cancels each other:
∂ℒ∂zk[L]=∑i=1C−yi·(δik−y^k)=−yk+y^k∑i=1Cyi Since ground truth value yi is one-hot vector (1 for true class and 0 for others), equation simplifies to:
∂ℒ∂zk[L]=y^k−yk This compact form is what makes softmax combined with cross-entropy loss so useful. Instead of dealing with complicated fractions, we can directly use this result to propagate gradients to the next layers.
From this point, gradient calculations across layers will start to follow a recurring pattern. Each layer essentially repeats the same process: we compute gradients with respect to its inputs and its weights.
For the layer output z[L], we have two key variables to differentiate with respect to:
The weights w[L]: their gradients are crucial because they are the trainable parameters we want to optimize during learning.
The activations a[L−1]: while not parameters themselves, their gradients are equally important since they ensure the flow of gradients backward through the network, enabling earlier layers to update as well. Thus, even though only the weight gradients directly influence optimization, the activation gradients play a important role in keeping backpropagation alive throughout the entire network:
- If we calculate derivatives with respect to weights w[l]: ∂ℒ∂w[l]=∂ℒ∂z[l]·∂(a[l−1]·w[l])∂w[l]=∂ℒ∂z[l]·a[l−1]
- If we calculate derivatives with respect to weights a[l−1]: ∂ℒ∂a[l−1]=∂ℒ∂z[l]·∂(a[l−1]·w[l])∂a[l−1]=∂ℒ∂z[l]·w[l] Lastly, we need to calculate the derivative of a[l] with respect to z[l] in σ(z[l]) (assuming σ(z[l])=ReLU(z[l])):
∂ℒ∂z[l]=∂ℒ∂a[l]·∂a[l]∂z[l]=∂ℒ∂a[l]·1(z[l]>0) where 1(z[l]>0) outputs 1 if (z[l]>0), 0 otherwise as a indicator function.
Therefore, we can summarize backpropagation for scalar values with following pattern:
∂ℒ∂z[L]=y^−y⋮∂ℒ∂z[l]=∂ℒ∂a[l]·∂a[l]∂z[l]=∂ℒ∂a[l]·1(z[l]>0)∂ℒ∂w[l]=∂ℒ∂z[l]·∂z[l]∂w[l]=∂ℒz[l]·a[l−1]∂ℒ∂a[l−1]=∂ℒ∂z[l]·∂z[l]∂a[l−1]=∂ℒz[l]·w[l]∂ℒ∂z[l−1]=∂ℒ∂a[l−1]·∂a[l−1]∂z[l−1]=∂ℒ∂a[l−1]·1(z[l−1]>0)⋮∂ℒ∂z[1]=∂ℒ∂a[1]·∂a[1]∂z[1]=∂ℒ∂a[1]·1(z[1]>0)∂ℒ∂w[1]=∂ℒ∂z[1]·∂z[1]∂w[1]=∂ℒz[1]·a[0]
3 - Backward Propagation
After computing the loss through forward propagation, we need to update the model’s weights to minimize this loss. Backpropagation is the algorithm that computes the gradients of the loss function with respect to each parameter in the network. We’ll derive these gradients step by step, working backwards from the loss to the input layer.
3.1 - Gradient of Cross-Entropy with Softmax
Starting with loss function, we calculate relative gradients and go back step by step to calculate further gradients. Since the model has predictions with Softmax activation and the loss is calculated with cross-entropy, the combination of these methods has a nice property which simplifies the gradient:
∂ℒ∂Z[L]=1N(Y^−Y)⏟[N×C] or in matrix representation:
∂ℒ∂Z[L]=1N[∂ℒ∂Z11[L]∂ℒ∂Z12[L]…∂ℒ∂Z1C[L]∂ℒ∂Z21[L]∂ℒ∂Z22[L]…∂ℒ∂Z2C[L]⋮⋮⋱⋮∂ℒ∂ZN1[L]∂ℒ∂ZN2[L]…∂ℒ∂ZNC[L]]=1N[Y^11−Y11Y^12−Y12…Y^1C−Y1CY^21−Y21Y^22−Y22…Y^2C−Y2C⋮⋮⋱⋮Y^N1−YN1Y^N2−YN2…Y^NC−YNC] or in visual representation:
3.2 - Gradient of Arbitrary Layer l
At this stage, we can begin calculating the gradients of the weights in the corresponding layer. Since the forward propagation is computed as
Z[l]=A[l−1]W[l] there are two possible local gradients we might consider:
∂Z[l]∂W[l] or∂Z[l]∂A[l−1] You may wonder why we would need the gradient with respect to A[l−1] since our goal is to update the weight matrix W[l]. The reason is that the derivative with respect to A[l−1] becomes necessary for updating the weights of the previous layer, because
A[l−1]=σ(A[l−2]W[l−1]) When we want to calculate the gradients of weight matrix W[l] with respect to the loss function ℒ, we use chain rule as follows:
∂ℒ∂W[l]=∂ℒ∂Z[l]⏟Upstream Gradient∂Z[l]∂W[l]⏟ Local Gradient The local gradient gives us:
∂Z[l]∂W[l]=∂(A[l−1]W[l])∂W[l]=A[l−1]⏟[N×Hl−1] Since there are Hl−1×Hl values in weight matrix W[l], we need to find each individual gradient so the shape of ∂ℒ∂W[l] should be Hl−1×Hl. However, the dimensions of matrices does not match when we try to calculate matrix multiplication:
∂ℒ∂W[l]⏟[Hl−1×Hl]=∂ℒ∂Z[l]⏟[N×Hl]A[l−1]⏟[N×Hl−1] If we take transpose of A[l−1] and place it to the left side of upstream gradient, we have the correct dimension matching:
∂ℒ∂W[l]⏟[Hl−1×Hl]=A[l−1]T⏟[Hl−1×N]∂ℒ∂Z[l]⏟[N×Hl] At first glance, this transformation may seem arbitrary. Why do we transpose? How can we be sure this result is correct? To clarify, let’s check the matrix representation:
[∂ℒ∂W11[l]∂ℒ∂W12[l]…∂ℒ∂W1Hl−1[l]∂ℒ∂W21[l]∂ℒ∂W22[l]…∂ℒ∂W2Hl−1[l]⋮⋮⋱⋮∂ℒ∂WHl1[l]∂ℒ∂WHl2[l]…∂ℒ∂WHlHl−1[l]]=[A11[l−1]TA12[l−1]T…A1N[l−1]TA21[l−1]TA22[l−1]T…A2N[l−1]T⋮⋮⋱⋮AHl−11[l−1]TAHl−12[l−1]T…AHl−1N[l−1]T][∂ℒ∂Z11[l]∂ℒ∂Z12[l]…∂ℒ∂Z1Hl[l]∂ℒ∂Z21[l]∂ℒ∂Z22[l]…∂ℒ∂Z2Hl[l]⋮⋮⋱⋮∂ℒ∂ZN1[l]∂ℒ∂ZN2[l]…∂ℒ∂ZNHl[l]] If we look closely only one gradient calculation, ∂ℒ∂W11[l]:
∂ℒ∂W11[l]=[A11[l−1]TA12[l−1]T…A1N[l−1]T][∂ℒ∂Z11[l]∂ℒ∂Z21[l]⋮∂ℒ∂ZN1[l]] Since single value of weight matrix is fixed for each input value, the gradient of this weight is just summation of all N input derivatives:
∂ℒ∂W11[l]=A11[l−1]T∂ℒ∂Z11[l]+A12[l−1]T∂ℒ∂Z21[l]+…+A1N[l−1]T∂ℒ∂ZN1[l] In the image below, you can see the visualization of this calculation for 4 weights in the network:
∂ℒ∂A[l−1]=∂ℒ∂Z[l]∂Z[l]∂A[l−1] The local gradient is:
∂Z[L]∂A[l−1]=∂(A[l−1]W[l])∂A[l−1]=W[l]⏟[Hl−1×Hl] Once again, the dimensions of matrices do not match:
∂ℒ∂A[l−1]⏟[N×Hl−1]=∂ℒ∂Z[l]⏟[N×Hl]W[l]⏟[Hl−1×Hl] We can transpose the weight matrix to align matrices:
∂ℒ∂A[l−1]⏟[N×Hl−1]=∂ℒ∂Z[l]⏟[N×Hl]W[l]T⏟[Hl×Hl−1] We can see the reasoning more easily if we check matrix representation again:
[∂ℒ∂A11[l−1]∂ℒ∂A12[l−1]…∂ℒ∂A1Hl−1[l−1]∂ℒ∂A21[l−1]∂ℒ∂A22[l−1]…∂ℒ∂A2Hl−1[l−1]⋮⋮⋱⋮∂ℒ∂AN1[l−1]∂ℒ∂AN2[l−1]…∂ℒ∂ANHl−1[l−1]]=[∂ℒ∂Z11[l]∂ℒ∂Z12[l]…∂ℒ∂Z1Hl[l]∂ℒ∂Z21[l]∂ℒ∂Z22[l]…∂ℒ∂Z2Hl[l]⋮⋮⋱⋮∂ℒ∂ZN1[l]∂ℒ∂ZN2[l]…∂ℒ∂ZNHl[l]][W11[l]TW12[l]T…W1Hl−1[l]TW21[l]TW22[l]T…W2Hl−1[l]T⋮⋮⋱⋮WHl1[l]TWHl2[l]T…WHlHl−1[l]T] By focusing on only first value, ∂ℒ∂A11[l−1], the gradient of A11l depends on the values that it affected during forward propagation (the same logic as ∂ℒ∂W11[l]) :
∂ℒ∂A11[l−1]=[∂ℒ∂Z11[l]∂ℒ∂Z12[l]…∂ℒ∂Z1Hl[l]][W11[l]TW21[l]T⋮WHl1[l]T] Now we want to compute the gradient ∂ℒ∂Z[l−1]. Recall that A[l−1] is obtained by applying an activation function to Z[l−1]:
A[l−1]⏟[N×Hl−1]=σ(Z[l−1]⏟[N×Hl−1]) In our examples, we use the ReLU activation function. Its derivative with respect to each element is given by:
∂ReLU(Zij[l−1])∂Zij[l−1]={1 ifZij[l]>0,0 otherwise For N×Hl−1 inputs, we must compute a gradient for each element, producing a matrix of 0s and 1s. The derivative of ReLU applied element-wise can be written as the indicator (mask) of positive entries:
∂ReLU(Z[l−1])∂Z[l−1]=1(Z[l−1]>0) where 1(·) is the element-wise indicator function (1 when the condition is true, 0 otherwise). In matrix form:
∂ReLU(Z[l−1])∂Z[l−1]=[1(Z11[l−1]>0)1(Z12[l−1]>0)…1(Z1Hl−1[l−1]>0)1(Z21[l−1]>0)1(Z22[l−1]>0)…1(Z2Hl−1[l−1]>0)⋮⋮⋱⋮1(ZN1[l−1]>0)1(ZN2[l−1]>0)…1(ZNHl−1[l−1]>0)] Because the derivative is applied element-wise, the upstream gradient ∂ℒ∂A[l−1] is combined with this mask element-wise, not by matrix multiplication:
∂ℒ∂Z[l−1]⏟[N×Hl−1]=∂ℒ∂A[l−1]⏟[N×Hl−1]⊙∂A[l−1]∂Z[l−1]⏟[N×Hl−1] or
∂ℒ∂Z[l−1]⏟[N×Hl−1]=∂ℒ∂A[l−1]⏟[N×Hl−1]⊙1(Z[l−1]>0)⏟[N×Hl−1]
3.3 - Repetition
At this point, we completed the backpropagation using the matrix multiplication step. From here on, the process is essentially a repetition: each layer applies the same logic and propagates the gradients backwards through the weights (also biases implicitly) and activations.
The only real exception was the output layer, where we combined the softmax function with cross-entropy loss. This required a more detailed derivation, but the model is consistent for all hidden layers and is repeated across the entire network. You can see the pattern as follows:
∂ℒ∂Z[L]⏟[N×C]=1N(Y^−Y)⏟N×C∂ℒ∂W[L]⏟[HL−1×C]=∂(A[L−1]W[L])∂W[L]⏟[HL−1×N]∂ℒ∂Z[L]⏟[N×C]=A[L−1]T⏟[HL−1×N]∂ℒ∂Z[L]⏟[N×C]∂ℒ∂A[L−1]⏟[N×HL−1]=∂ℒ∂Z[L]⏟[N×C]∂(A[L−1]W[L])∂A[L−1]⏟[C×HL−1]=∂ℒ∂Z[L]⏟[N×C]W[L]T⏟[C×HL−1]∂ℒ∂Z[L−1]⏟[N×HL−1]=∂ℒ∂A[L−1]⏟[N×HL−1]⊙∂(σ(Z[L−1]))∂Z[L−1]⏟[N×HL−1]=∂ℒ∂A[L−1]⏟[N×HL−1]⊙1(Z[L−1]>0)⏟[N×HL−1]⋮∂ℒ∂W[l]⏟[Hl−1×Hl]=∂(A[l−1]W[l])∂W[l]⏟[Hl−1×N]∂ℒ∂Z[l]⏟[N×Hl]=A[l−1]T⏟[Hl−1×N]∂ℒ∂Z[l]⏟[N×Hl]∂ℒ∂A[l−1]⏟[N×Hl−1]=∂ℒ∂Z[l]⏟[N×Hl]∂(A[l−1]W[l])∂A[l−1]⏟[Hl×Hl−1]=∂ℒ∂Z[l]⏟[N×Hl]W[l]T⏟[Hl×Hl−1]∂ℒ∂Z[l−1]⏟[N×Hl−1]=∂ℒ∂A[l−1]⏟[N×Hl−1]⊙∂(σ(Z[l−1]))∂Z[l−1]⏟[N×Hl−1]=∂ℒ∂A[l−1]⏟[N×Hl−1]⊙1(Z[l−1]>0)⏟[N×Hl−1]⋮∂ℒ∂Z[1]⏟[N×H1]=∂ℒ∂A[1]⏟[N×H1]⊙∂(σ(Z[1]))∂Z[1]⏟[N×H1]=∂ℒ∂A[1]⏟[N×H1]⊙1(Z[1]>0)⏟[N×H1]∂ℒ∂W[1]⏟[H0×H1]=∂(A[0]W[1])∂W[1]⏟[H0×N]∂ℒ∂Z[1]⏟[N×H1]=A[0]T⏟[H0×N]∂ℒ∂Z[1]⏟[N×H1] I know it looks ugly but if you pay some attention, you will see the pattern for gradient calculation. There are just some rules that you need to follow for each block.
For comments, please send me an email.