additive form over a set of training data $\{ (x_i, y_i) \}_{i=1}^n$ as follows: $$\begin{aligned}
    \rho(w) &= \frac{1}{n} \sum_{i=1}^n f_i( w ), \label{eq:fidelity}
    \end{aligned}$$ where $f_i( w )$ denotes a data fidelity incurred by a set of model parameters $w$ for a sample pair $(x_i, y_i)$. The data fidelity $f_i( w )$ is designed to measure the discrepancy between the prediction $h_w( x_i )$ with an input $x_i$ and its desired output $y_i$ for a given sample pair $(x_i, y_i)$. The regularization $\gamma( w )$ is designed to impose a smoothness constraint to the solution space, thus avoid undesirable over-fitting to the model. The weight-decay coefficient $\lambda \in \R$ is determined based on the relation between the underlying distribution of the data and the prior distribution of the model. In the optimization of the objective function $\mathcal{L}(w)$ defined by: $$\begin{aligned}
    \mathcal{L}(w) &= \frac{1}{n} \sum_{i=1}^n f_i( w ) + \lambda \, \gamma( w ), \label{eq:objective}
    \end{aligned}$$ where $f_i(w)$ and $\gamma(w)$ are assumed to be differentiable, we consider a first-order optimization algorithm leading to the following gradient descent step at each iteration $t$: $$\begin{aligned}
    w^{t+1} \coloneqq w^{t} - \eta^{t} \left( \frac{1}{n} \sum_{i=1}^n \nabla f_{i}(w^{t}) + \lambda \nabla \gamma(w^t) \right), \label{eq:update:vanilla}
    \end{aligned}$$ where 