The original question asked whether the error function needs to be convex. No, it does not. The analysis presented below is intended to provide some insight and intuition about this and the modified question, which asks whether the error function could have multiple local minima.
Intuitively, there doesn't have to be any mathematically necessary relationship between the data and the training set. We should be able to find training data for which the model initially is poor, gets better with some regularization, and then gets worse again. The error curve cannot be convex in that case--at least not if we make the regularization parameter vary from $0$ to $\infty$.
Note that convex is not equivalent to having a unique minimum! However, similar ideas suggest multiple local minima are possible: during regularization, first the fitted model might get better for some training data while not appreciably changing for other training data, and then later it will get better for other training data, etc. A suitable mix of such training data ought to produce multiple local minima. To keep the analysis simple I won't attempt to show that.
Edit (to respond to the changed question)
I was so confident in the analysis presented below and the intuition behind it that I set about finding an example in the crudest possible way: I generated small random datasets, ran a Lasso on them, computed the total squared error for a small training set, and plotted its error curve. A few attempts produced one with two minima, which I will describe. The vectors are in the form $(x_1,x_2,y)$ for features $x_1$ and $x_2$ and response $y$.
Training data
$$(1,1,-0.1),\ (2,1,0.8),\ (1,2,1.2),\ (2,2,0.9)$$
Test data
$$(1,1,0.2),\ (1,2,0.4)$$
The Lasso was run using glmnet::glmmet
in R
, with all arguments left at their defaults. The values of $\lambda$ on the x axis are the reciprocals of the values reported by that software (because it parameterizes its penalty with $1/\lambda$).
An error curve with multiple local minima
Analysis
Let's consider any regularization method of fitting parameters $\beta=(\beta_1, \ldots, \beta_p)$ to data $x_i$ and corresponding responses $y_i$ that has these properties common to Ridge Regression and Lasso:
(Parameterization) The method is parameterized by real numbers $\lambda \in [0, \infty)$, with the unregularized model corresponding to $\lambda=0$.
(Continuity) The parameter estimate $\hat\beta$ depends continuously on $\lambda$ and the predicted values for any features vary continuously with $\hat\beta$.
(Shrinkage) As $\lambda\to\infty$, $\hat\beta\to 0$.
(Finiteness) For any feature vector $x$, as $\hat\beta\to 0$, the prediction $\hat y(x) = f(x, \hat\beta) \to 0$.
(Monotonic error) The error function comparing any value $y$ to a predicted value $\hat y$, $\mathcal{L}(y, \hat y)$, increases with the discrepancy $|\hat y - y|$ so that, with some abuse of notation, we may express it as $\mathcal{L}(|\hat y - y|)$.
(Zero in $(4)$ could be replaced by any constant.)
Suppose the data are such that the initial (unregularized) parameter estimate $\hat\beta(0)$ is not zero. Let's construct a training data set consisting of one observation $(x_0, y_0)$ for which $f(x_0, \hat\beta(0))\ne 0$. (If it's not possible to find such an $x_0$, then the initial model won't be very interesting!) Set $y_0=f(x_0, \hat\beta(0))/2$.
The assumptions imply the error curve $e: \lambda \to \mathcal{L}(y_0, f(x_0, \hat\beta(\lambda))$ has these properties:
$e(0) = \mathcal{L}(y_0, f(x_0, \hat\beta(0)) = \mathcal{L}(y_0, 2y_0) = \mathcal{L}(|y_0|)$ (because of the choice of $y_0$).
$\lim_{\lambda\to\infty}e(\lambda) = \mathcal{L}(y_0, 0) = \mathcal{L}(|y_0|)$ (because as $\lambda\to\infty$, $\hat\beta(\lambda)\to 0$, whence $\hat{y}(x_0)\to 0$).
Thus, its graph continuously connects two equally high (and finite) endpoints.
Qualitatively, there are three possibilities:
The prediction for the training set never changes. This is unlikely--just about any example you choose will not have this property.
Some intermediate predictions for $0\lt \lambda \lt \infty$ are worse than at the start $\lambda=0$ or in the limit $\lambda\to\infty$. This function cannot be convex.
All intermediate predictions lie between $0$ and $2y_0$. The continuity implies there will be at least one minimum of $e$, near which $e$ must be convex. But since $e(\lambda)$ approaches a finite constant asymptotically, it cannot be convex for large enough $\lambda$.
The vertical dashed line in the figure shows where the plot changes from convex (at its left) to non-convex (to the right). (There is also a region of non-convexity near $\lambda\approx 0$ in this figure, but this won't necessarily be the case in general.)