According to wikipedia of backpropagation:
In fitting a neural network, backpropagation computes the gradient of the loss function during supervised learning with respect to the weights of the network for a single input–output example, and does so efficiently, unlike a naive direct computation of the gradient with respect to each weight individually.
Backpropagation is a special case of reverse accumulation of automatic differentiation, and it was announced in Rumelhart, Hinton & Williams (1986). Automatic differentiation has 2 modes: forward accumulation and backward accumulation.
Automatic differentiation is distinct from symbolic differentiation and numerical differentiation (the method of finite differences). Symbolic differentiation can lead to inefficient code and faces the difficulty of converting a computer program into a single expression, while numerical differentiation can introduce round-off errors in the discretization process and cancellation. Both classical methods have problems with calculating higher derivatives, where complexity and errors increase. Finally, both classical methods are slow at computing partial derivatives of a function with respect to many inputs, as is needed for gradient-based optimization algorithms. Automatic differentiation solves all of these problems.
Usually, two distinct modes of AD are presented, forward accumulation (or forward mode) and reverse accumulation (or reverse mode). Forward accumulation specifies that one traverses the chain rule from inside to outside, while reverse accumulation has the traversal from outside to inside...
Forward accumulation is more efficient than reverse accumulation for functions $ f:ℝ^n → ℝ^m$ with $m ≫ n$ as only $n$ sweeps are necessary, compared to $m$ sweeps for reverse accumulation.
Reverse accumulation is more efficient than forward accumulation for functions $ f:ℝ^n → ℝ^m$ with $m ≪ n$ as only $m$ sweeps are necessary, compared to $n$ sweeps for forward accumulation.
Under this historical computational differentiation techniques context, since the loss function in a MLP maps $ℝ^n \rightarrow ℝ^m$ with $m ≪ n$ for most cases, that's one of the main reasons why backpropagation is most often used to train a MLP as a special case of reverse accumulation introduced above.