We often train neural networks by optimizing the mean squared error (MSE), which is an equation of a parabola $y=x^2$, with gradient descent.
We also say that weight adjustment in a neural network by the gradient descent algorithm can hit a local minimum and get stuck in there.
How are multiple local minima on the equation of a parabola possible, if a parabola has only one minimum?
$g(x) = x^2$ is indeed a parabola and thus has just one optimum.
However, the $\text{MSE}(\boldsymbol{x}, \boldsymbol{y}) = \sum_i (y_i - f(x_i))^2$, where $\boldsymbol{x}$ are the inputs, $\boldsymbol{y}$ the corresponding labels and the function $f$ is the model (e.g. a neural network), is not necessarily a parabola. In general, it is only a parabola if $f$ is a constant function and the sum is over one element.
For example, suppose that $f(x_i) = c, \forall i$, where $c \in \mathbb{R}$. Then $\text{MSE}(\boldsymbol{x}, \boldsymbol{y}) = \sum_i (y_i - c)^2$ will only change as a function of one variable, $\boldsymbol{y}$, as in the case of $g(x) = x^2$, where $g$ is a function of one variable, $x$. In that case, $(y_i - c)^2$ will just be a shifted version (either to the right or left depending on the sign of $c$) of $y_i^2$, so, for simplicity, let's ignore $c$. So, in the case $f$ is a constant function, then $\text{MSE}(\boldsymbol{x}, \boldsymbol{y}) = \sum_i y_i^2$, which is a sum of parabolas $y_i^2$, which is called a paraboloid. In this case, the paraboloid corresponding to $\text{MSE}(\boldsymbol{x}, \boldsymbol{y}) = \sum_i y_i^2$ will only have one optimum, just like a parabola. Furthermore, if the sum is just over one $y_i$, that is, $\text{MSE}(\boldsymbol{x}, \boldsymbol{y}) = \sum_i y_i^2 = y^2$ (where $\boldsymbol{y} = y$), then the MSE becomes a parabola.
In other cases, the MSE might not be a parabola or have just one optimum. For example, suppose that $f(x) = x^2$, $y_i = 1$ ($\forall i$), then $h(x) = (1 - x^2)^2$ looks as follows
which has two minima at $x=-1$ and $x=1$ and one maximum at $x=0$. We can find the two minima of this function $h$ using calculus: $h'(x) = -4x(1 - x^2)$, which becomes zero when $x=-1$ and $x=1$.
In this case, we only considered one term of the sum. If we considered the sum of terms of the form of $h$, then we could even have more "complicated" functions.
To conclude, given that $f$ can be arbitrarily complex, then also $\text{MSE}(\boldsymbol{x}, \boldsymbol{y})$, which is a function of $f$, can also become arbitrarily complex and have multiple minima. Given that neural networks can implement arbitrarily complex functions, then $\text{MSE}(\boldsymbol{x}, \boldsymbol{y})$ can easily have multiple minima. Moreover, the function $f$ (e.g. the neural network) changes during the training phase, which might introduce more complexity, in terms of which functions the MSE can be and thus which (and how many) optima it can have.
How are multiple local minima on the equation of a parabola possible, if a parabola has only one minimum?
A parabola has one minimum, and no separate local minima. So it isn't possible.
However...
Gradient descent works on the equation of mean squared error, which is an equation of a parabola $y=x^2$
Just because the loss function is a parabola with respect to the direct input, does not mean that the loss function is a parabola with respect to the parameters that indirectly cause that error.
In fact it only remains true for linear functions. When considering linear regression $\hat{y} = \sum_i w_i x_i + b$, there is only one global minimum (with specific values of $w_i$ or specific vector $\mathbf{w}$), and your assertion is true.
Once you add nonlinear activations, as in neural networks, then the relationship between error function and parameters of the model becomes far more complex. For the last/output layer you can carefully choose a loss function so that this cancels out - you can keep your single global minimum for logistic regression and softmax regression. However, one or more hidden layers, and all bets are off.
In fact you can prove quite easily that a neural network with a hidden layer must have multiple stationary points (not necessarily local minima). The outline of the proof is to note that there must be multiple equivalent solutions, since in a fully-connected network you can re-arrange the nodes into any order, move the weights to match, and it will be a new solution with exactly the same behaviour, including the same loss on the dataset. So a neural network with one hidden layer with $n$ nodes must have $n!$ absolute minimums. There is no way for these to exist without other stationary points in-between them.
There is theory to suggest that most of the stationary points found in practice will not be local minima, but saddle points.
As an example, this is an analysis of saddle points in a simple XOR approximator.
Another view on this topic: think about the derivative of the MSE with respect to the inputs. You will need to apply the chain rule:
$$\frac{dMSE}{dx_i}=\sum_i2(y_i-f(x_i))\cdot(-1)\cdot\frac{df(x_i)}{dx_i}$$
which only resembles the derivative of the parabola when, as mentioned in nbro's answer, $f(x_i)=c, \forall\ i$ or, more generally, when they $f(x)$ is linear (Neil's answer).
