Is backpropagation a fancy way of saying "calculate gradient by taking partial derivative w.r.t. all x's" I understanding how gradient is calculated in the usual context--it is just taking partial derivative w.r.t. each element of the X vector. Say a function $f$ has $n$ independent variables, denoted by $x_1$ ... $x_n$, the gradient of $f$ is:
$$
\nabla f = \left<\frac{\partial f}{\partial x_1};\frac{\partial f}{\partial x_2};...;\frac{\partial f}{\partial x_n}\right>
$$
However, when it comes to the gradient of the loss function of a neural network, a new term, backpropagation, emerges. My confusion is I am not sure if backpropagation makes the above calculation different or is it just another name of it.
Let me use the following simple neural network as an example (courtesy of this very informative video)：

This NN has two nodes, four weights (denoted by $w_1$ ... $w_4$) and three biases (denoted by $b_1$ ... $b_3$). So the loss function is defined like this:
$L(w_1, w_2, w_3, w_4, b_1, b_2, b_3)$ = [a loss function, which can be extremely lengthy but differentiable]
Without knowing the exact meaning of backpropagation, I would say its gradient is just:
$$
\nabla L = \left<\frac{\partial L}{\partial w_1};\frac{\partial L}{\partial w_2};...;\frac{\partial L}{\partial b_3}\right>
$$
and the calculation of the partial derivative w.r.t. $w_1$, $w_2$ ... $b_2$, $b_3$ is no more than repetitive application of chain rules, etc--we don't need anything more than good old calculus to make it work. With gradient calculated, we can use gradient descent to minimize the loss function (well, the convexity/local minima is another big issue, let's ignore it for the purpose of this particular question).
So it begs the question: what is the role of "backpropagation" in the calculation? Does backpropagation make the calculation easier (i.e., backpropagation is faster but the resultant gradient will be exactly the same as the traditional partial derivative approach) or is it just a name of this calculation?

Here are some examples to illustrate how backpropagation (automatic differentiation) is quite nontrivial, and the tradeoffs involved.
Consider $x_n = f_n(x_{n-1}) = f_n(f_{n-1}(x_{n-2})) = f_n(f_{n-1}(\ldots f_1(x_0)))$. Let $J_i(x)$ denote the jacobian of $f_i$ at $x$. Then it's true that $\frac{d x_n}{d x_0} = J_n(x_{i-1}) J_{n-1}(x_{i-1}) \ldots J_1(x_0)$.
But what order should you do the multiplication in? For example, if $n = 3$, and the dimension of $x_0, x_1, x_2, x_3$ are 1000, 100, 10, 1 respectively, then it makes sense to multiply from left to right, so the most expensive matrix multiplication you'll do is a $1 \times 100$ by $100 \times 1000$. If you go right to left, you'll end up doing a $10 \times 100$ by $100 \times 1000$ matmul, which is 10 times slower. The former method is called "reverse mode", and the latter is called "forward mode". Backpropagation is typically implemented using reverse mode, because it is more suited to neural network training (and in general, when the number of inputs is greater than the number of outputs).
Now consider if the dimension of $x_0, x_1, x_2, x_3, x_4$ are 1000, 1, 1000, 1, 1000. Both forward and reverse mode accumulation will result in a $1000 \times 1000$ by $1000 \times 1$ product, whereas if you multiply the two middle jacobians first, you only end up doing $1000 \times 1$ by $1 \times 1$. So this alternate multiplication order is much more efficient. Ironically, there is no efficient algorithm for determining the most efficient order - this is an NP-complete problem.
Another consideration is that typical jacobians are extremely sparse: Applying pointwise sigmoid activation to an vector of size 1 million creates a $1000000 \times 1000000$ jacobian matrix, with only 1 million non-zero values. It's infeasible to multiply jacobians this large, but to carry out backpropagation, it's actually only necessary to compute the "vector jacobian product": $\text{vjp}_i(v,x) = v^TJ_i(x)$. Then if $x_n$ has dimension 1 (i.e any loss function), $\frac{d x_n}{d x_0} = \text{vjp}_1( \text{vjp}_2( \ldots\text{vjp}_n(1,x_{n-1})\ldots, x_1), x_0)$ carries out reverse mode backprop. For a pointwise function $g$, $\text{vjp}_g(v,x) = v\ \circ \ g'(x)$ is trivial to compute.
So now suppose we've implemented blazing fast reverse mode automatic differentiation with vjps, and we want to minimize a function using second order information from the hessian matrix. We could compute the hessian by first computing the gradient $\nabla(x)$, and then backpropagate the gradient calculation itself: $H_i = \text{vjp}_{\nabla_i}(1, x)$.
Note that because we're using vjps / reverse mode backprop, we can only compute one row of the hessian at a time - as noted above, reverse mode is poorly suited when the dimension of the output is large. Ideally, we'd want to use reverse mode for computing the gradient, but forward mode to differentiate the gradient. As it turns out, you can implement a "jacobian vector product": $\text{jvp}_i(x,v) = J_i(x)v$ as two invocations of vjp (link below). A vjp is to reverse mode as jvp is to forward mode, and this enables us to efficiently perform second order optimization.
https://j-towns.github.io/2017/06/12/A-new-trick.html
