How are the gradients of a Neural Network calculated just by matrix multiplication? I would have expected some kind of derivative solving equation to be at work in order to back propagate the loss to each neuron. I hope my question is not too confused to answer.
In the network below, why does h_relu.T @ grad_y_pred (y hat - y) produce the gradient for w2? I understand sometimes the derivative is just something you multiply something by. Like the derivative of y = h_relu * w2 being h_relu, which doesn't make much sense to me either. I understand why the derivative of x^2 is 2x and that's about as far as my understanding of derivatives goes.

import numpy as np

# N is batch size; D_in is input dimension;
# H is hidden dimension; D_out is output dimension.
N, D_in, H, D_out = 64, 1000, 100, 10

# Create random input and output data
x = np.random.randn(N, D_in)
y = np.random.randn(N, D_out)

# Randomly initialize weights
w1 = np.random.randn(D_in, H)
w2 = np.random.randn(H, D_out)

learning_rate = 1e-6
for t in range(500):
    # Forward pass: compute predicted y
    h = x.dot(w1)
    h_relu = np.maximum(h, 0)
    y_pred = h_relu.dot(w2)

    # Compute and print loss
    loss = np.square(y_pred - y).sum()
    print(t, loss)

    # Backprop to compute gradients of w1 and w2 with respect to loss
    grad_y_pred = 2.0 * (y_pred - y)
    grad_w2 = h_relu.T.dot(grad_y_pred)
    grad_h_relu = grad_y_pred.dot(w2.T)
    grad_h = grad_h_relu.copy()
    grad_h[h &lt; 0] = 0
    grad_w1 = x.T.dot(grad_h)

    # Update weights
    w1 -= learning_rate * grad_w1
    w2 -= learning_rate * grad_w2
Like the derivative of y = h_relu * w2 being h_relu, which doesn't make much sense to me 


These are partial derivative and depends on w.r.t variable. In the end, we have to complete the chain of partial derivative from output to the respective weights. 

Your Network is -
$$x \xrightarrow[\text{}]{\text{}} w_1 \xrightarrow[\text{}]{\text{  h  }} ReLU \xrightarrow[\text{}]{\text{  h_relu  }} w_2 \xrightarrow[\text{}]{\text{ }} y\_pred$$
By definition ReLU is - $f(x)=\max(0, x).$
Derivative of ReLU - $\frac{d}{dx}\text{ReLU}(x)$
=$1 \forall x &gt; 0$
=$0 \forall x &lt; 0$

In backpropagation -
$\frac{\partial y\_pred}{\partial w_2} $ = $h\_relu$  ;since, y_pred = $h\_relu*w_2$
$\frac{\partial w_2}{\partial h} $ = $1$  ;derivative of ReLU=1,ignoring negative part for simplicity
$\frac{\partial h}{\partial w_1} $ = $x$ ;since, h = $x*w_1$
$\frac{\partial y\_pred}{\partial h\_relu} $ = $w_2$ ;since, y_pred = $h\_relu*w_2$

So, Loss gradient at different weights will be - 
 $L$ at the $y\_pred$
 $L*h\_relu$ at $w_2$ 
 $L*w_2$ at $h\_relu$ 
 $L*h\_relu*1*x$ at $w_1$