Outlier-Robust Neural Network Training: Efficient Optimization of
  Transformed Trimmed Loss with Variation Regularization
In this study, we consider outlier-robust predictive modeling using
highly-expressive neural networks. To this end, we employ (1) a transformed
trimmed loss (TTL), which is a computationally feasible variant of the
classical trimmed loss, and (2) a higher-order variation regularization (HOVR)
of the prediction model. Note that using only TTL to train the neural network
may possess outlier vulnerability, as its high expressive power causes it to
overfit even the outliers perfectly. However, simultaneously introducing HOVR
constrains the effective degrees of freedom, thereby avoiding fitting outliers.
We newly provide an efficient stochastic algorithm for optimization and its
theoretical convergence guarantee. (*Two authors contributed equally to this
work.)
Akifumi Okuno, Shotaro Yagishita
# A stochastic optimization approach

###### Abstract

While highly expressive parametric models including deep neural networks have an advantage to model complicated concepts, training such highly non-linear models is known to yield a high risk of notorious overfitting. To address this issue, this study considers a \((k,q)\)th order variation regularization (\((k,q)\)-VR), which is defined as the \(q\)th-powered integral of the absolute \(k\)th order derivative of the parametric models to be trained; penalizing the \((k,q)\)-VR is expected to yield a smoother function, which is expected to avoid overfitting. Particularly, \((k,q)\)-VR encompasses the conventional (general-order) total variation with \(q=1\). While the \((k,q)\)-VR terms applied to general parametric models are computationally intractable due to the integration, this study provides a stochastic optimization algorithm, that can efficiently train general models with the \((k,q)\)-VR without conducting explicit numerical integration. The proposed approach can be applied to the training of even deep neural networks whose structure is arbitrary, as it can be implemented by only a simple stochastic gradient descent algorithm and automatic differentiation. Our numerical experiments demonstrate that the neural networks trained with the \((k,q)\)-VR terms are more "resilient" than those with the conventional parameter regularization. The proposed algorithm also can be extended to the physics-informed training of neural networks (PINNs).

Figure 1: Non-linear neural network models (defined as a single-hidden layer perception with \(L=200\) hidden units) are trained with \(N=100\) samples that follow a quadratic function. The neural network trained with the \((\mathrm{c})\)\((3,2)\)-VR term with \(\eta_{3}=10^{-5}\) is more “resilient” (i.e., the number of inflection points is smaller) than those with (a) no regularization and (b) \(L_{2}\) regularization. All the optimization settings including the parameter initialization, decay schedule of the learning rate, minibatch sizes (\(n=m=5\)), and the number of iterations (\(T=2\times 10^{4}\)) for minibatch stochastic gradient descent, are the same. See Section 3 for more details, including the same experiments for linear and cubic functions.
Introduction

It is needless to emphasize that highly non-linear parametric models including deep neural network have attracted considerable attention these days, by virtue of their high representation capability (Goodfellow et al., 2016; Sze et al., 2017; Miikkulainen et al., 2019; Samek et al., 2021). They are proven to approximate arbitrary continuous functions (see, e.g., Cybenko (1989) for shallow and Yarotsky (2017) for deep neural networks); the highly non-linear models are expected to fit the underlying target functions adaptively. While the high expressive power can become a strong merit to capture the complicated concepts (e.g., large language models (Brown et al., 2020; Liu et al., 2023) used to capture the human language structure), training such highly non-linear models with a relatively smaller number of samples is known to yield notorious overfitting issues. See Figure 1(a) for a neural network trained without any regularization. The neural network seems to have a prodigal number of inflection points that are not welcome for better prediction and interpretation.

For a better prediction, various approaches have been proposed to avoid overfitting in the context of neural network training. To name a few, dropout (Srivastava et al., 2014), early-stopping (see, e.g., references in Yao et al. (2007)), data-augmentation (see, e.g., Hernandez-Garcia and Konig (2018) and Zhang et al. (2021) for the effectiveness), batch-normalization (Luo et al., 2019), and so forth. They are considered to be variants of regularization, and all of the above regularization approaches are compatible with more direct parameter regularization. The \(L_{2}\) parameter regularization is also known as ridge regularization, Tikhonov regularization or weight decay (Kingma and Ba, 2014; Goodfellow et al., 2016), and it is further generalized to \(L_{p}\) regularization including lasso-type (Tibshirani, 1996) and elastic-type (Zou and Hastie, 2005) regularizations. See Figure 1(b) for a neural network trained with \(L_{2}\) parameter regularization. Unfortunately, the neural network trained with \(L_{2}\) penalization still may contain a prodigal number of inflection points as the parameter regularization is not directly related to the output variation.

For additive models \(f_{\theta}^{(\text{add.})}(x)=\sum_{j=1}^{N}\theta_{j}\phi_{j}(x)\) defined with user-specified basis functions \(\{\phi_{j}\}\), that includes kernel and spline regression models as special cases, the \(L_{2}\)-type regularization term \(\langle\theta,G\theta\rangle\) for some matrix \(G\) coincides with a function norm \(\|f_{\theta}^{(\text{add.})}\|_{L_{2}}^{2}=\int\{f_{\theta}^{(\text{add.})}( x)\}^{2}\mathrm{d}x\) (see Appendix A). It is further generalized to the \((k,q)\)th order variation regularization (\((k,q)\)-VR) considered in this study:

\[C_{k,q}(f_{\theta}):=\frac{1}{2}\int_{\Omega}\Big{|}\frac{\partial^{k}f_{ \theta}(x)}{\partial x^{k}}\Big{|}^{q}\mathrm{d}x.\]

\((k,q)\)-VR encompasses a total variation (TV) regularization (Rudin et al., 1992; Engl et al., 1996; Osher et al., 2005), which corresponds to \((k,q)=(1,1)\), and TV has been incorporated into the training of the regression models equipped with splines (Stone, 1994; Stone et al., 1997; Mammen and van de Geer, 1997), kernels (Zou and Hastie, 2005), triograms (Koenker and Mizera, 2004), and Delaunay triangles (Pourya et al., 2023). The 1st order TV (\(k=1,q=1\)) regularization has been further generalized to the TV regularization of 2nd order (\(k=2\)) (Koenker and Mizera, 2004; Hinterberger and Scherzer, 2006; Duan et al., 2015) and general order (\(k\in\mathbb{N}\)) (Bredies and Holler, 2020). While the TV regularization can be simply incorporated into the training of such additive models whose basis functions are user-specified, the optimization techniques cannot be straightforwardly generalized to neural networks which adaptively learn the basis functions in their training.

The main problem to introduce the \((k,q)\)-VR to neural network training is the computational intractability of the integral. To avoid the integration, neural splines (Williams et al., 2021) approximate the 2nd order TV terms by their finite approximation variants; a similar idea can be found in a variety of existing studies (see, e.g., Koenker and Mizera (2004) equipped with a theoretical justification shown in Natanson (1974) Theorem IX.4.8), and it can be also regarded as computing numerical integration instead of the exact integration itself. Generally speaking, numerical integration requires high computational complexity to obtain high precision. Another interesting direction is attempted by a purely theoretical work (Unser, 2019). Therein, 2nd order TV regarding the activation function is considered for deep neural network training and proves the optimality of a piece-wise linear activation function (see Theorem 4). However, as also noted in Unser (2019) section 3.4, finding the optimal piece-wise linear function still remains difficult practically. The representer theorem is also considered in Banach space (Parhi and Nowak, 2021) with TV terms measured in the Radon domain. Their theories have been further developed mathematically (Parhi and Nowak, 2022; Unser, 2023; Bartolucci et al., 2023). While these works have provided significant progresses in the theoretical understanding of neural network behaviors, practical algorithm for training neural networks with exact \((k,q)\)-VR including TV regularization as a special case, still remains lacking.
To address these issues, this study proposes a stochastic optimization algorithm in line with Robbins and Monro (1951) and Ghadimi and Lan (2013), by following the similar idea used to optimize the intractable likelihood (Geyer and Thompson, 1992), contrastive divergence (Carreira-Perpinan and Hinton, 2005) and robust divergence (Okuno, 2023). The proposed stochastic optimization algorithm efficiently minimizes the loss function equipped with the \((k,q)\)-VR terms without conducting numerical integration. See Figure 1(c). The neural network trained with \((k,q)\)-VR terms is more "resilient" than those with \(L_{2}\) regularization, and it is more suitable for prediction and interpretation purposes. The proposed algorithm can be applied to general parametric models including arbitrary-structured neural networks, as the algorithm is a slight modification of the conventional stochastic gradient descent. Furthermore, the proposed algorithm also can be applied to the physics-informed training of a neural network (PINN; Dissanayake and Phan-Thien, 1994; Berg and Nystrom, 2018; Cuomo et al., 2022). Our approach can minimize the PINN loss function including the integral-type constraint directly, without computing the explicit numerical integration. See Section 4 for a brief discussion, though PINN extension slightly goes beyond the main scope of this study.

## 2 Higher-order variation regularization

Section 2.1 describes the conventional regression framework and the remaining issues to be discussed in this study. Section 2.2 describes the higher-order variation regularization, and Section 2.3 describes the proposed stochastic algorithm. While this study considers univariate functions for simplicity, the below discussions can be generalized to the multivariate case straightforwardly.

### Conventional regression framework

Let \(D_{1},D_{2},N\in\mathbb{N}\), \(\Omega:=[D_{1},D_{2}]\) and let \((x_{i},y_{i})\in\Omega\times\mathbb{R}\) be a pair of observations \((i=1,2,\ldots,N)\). Let \(f_{\theta}:\Omega\rightarrow\mathbb{R}\) be a regression model equipped with the parameter \(\theta\in\Theta\subset\mathbb{R}^{r}\) (\(r\in\mathbb{N}\)). We consider training the regression model by estimating the parameter \(\theta\), so that \(f_{\theta}(x_{i})\) well approximates the observed outcome \(y_{i}\). While this study employs a single-hidden-layer perceptron:

\[f_{\theta}(x)=\sum_{\ell=1}^{L}a_{\ell}\sigma(b_{\ell}x+c_{\ell})+d,\quad \theta=(\{a_{\ell}\}_{\ell},\{b_{\ell}\}_{\ell},\{c_{\ell}\}_{\ell},d) \tag{1}\]

for simplicity, the below discussion can be extended to arbitrary regression models including deep neural networks and kernel regression models. \(\sigma:\mathbb{R}\rightarrow\mathbb{R}\) denotes an activation function. Typically, we employ a sigmoid function \(\sigma(z)=1/\{1+\exp(-z)\}\) or a hyperbolic tangent function \(\sigma(z)=\{\exp(z)-\exp(-z)\}/\{\exp(z)+\exp(-z)\}\).

We may estimate the parameter \(\theta\) by minimizing a loss function

\[A(\theta)=\frac{1}{N}\sum_{i=1}^{N}\nu(y_{i}-f_{\theta}(x_{i}))\]

defined with some non-negative function \(\nu:\mathbb{R}\rightarrow\mathbb{R}_{\geq 0}\). Typically, we employ a quadratic function \(\nu(z)=z^{2}\), and we may replace this loss function to, for instance, Tukey's biweight function (Huber and Ronchetti, 1981) and so forth depending on the regression purposes.

The single-hidden-layer perceptron (1) is known to have a universal approximation capability. Namely, (1) can approximate arbitrary continuous functions by increasing the number of hidden units \(L\) to infinity (see, e.g., Cybenko (1989)) as well as more recent deep neural networks (Yarotsky, 2017). While such high expressiveness is greatly welcomed to approximate the underlying functions adaptively, the highly expressive models usually contain a prodigal number of inflection points as illustrated in Figure 1(a), and it results in overfitting to the training samples. Owing to their observational errors, overfitting to only the observed training samples is known to degrade the prediction accuracy for unseen data. To address such an overfitting issue, we may consider a parameter regularization. Typically, we may employ a \(L_{2}\) regularization (also known as ridge regularization, Tikhonov regularization and weight decay)

\[B(\theta)=\frac{1}{2}\|\theta\|_{2}^{2}=\frac{1}{2}[\theta_{1}^{2}+\theta_{2} ^{2}+\cdots+\theta_{r}^{2}].\]
Then, we may simply minimize the regularized loss function \(\tilde{L}_{A}(\theta)=A(\theta)+\lambda B(\theta)\) to train the regression model \(f_{\theta}\). The minimizer is explicitly obtained for linear regression models, and even for non-linear neural networks, the loss function \(\tilde{L}_{A}(\theta)\) is easily minimized by leveraging gradient-based algorithms. Typically, a stochastic gradient descent algorithm is employed for optimization (see, e.g., adam optimizer (Kingma and Ba, 2014) for implementation). While such a parameter regularization is known to relieve the overfitting issue, the trained model still contains a number of inflection points as illustrated in Figure 1(b). While the parameter regularization can be regarded as penalization of the output variation for additive regression models (see Appendix A), the regularization seems rather indirect for more highly expressive models including neural networks; more direct regularization for the higher-order output variation is appreciated.

### Higher-order variation regularization

To address the aforementioned issues, this study employs a \((k,q)\)th order variation regularization (\((k,q)\)-VR) for the regression model \(f_{\theta}\) itself:

\[C_{k,q}(f_{\theta})=\frac{1}{2}\int_{\Omega}\big{|}f_{\theta}^{[k]}(x)\big{|} ^{q}\mathrm{d}x. \tag{2}\]

\(k\in\{0,1,2,\ldots,\}\) and \(q>0\) are user-specified parameters, and \(C_{k,q}(f)\) regularizes the \(q\)th powered variation of the \(k\)th derivative \(f^{[k]}(x)=\partial^{k}f(x)/\partial x^{k}\). \((1,1)\)-VR is also known as a total variation (Rudin et al., 1992; Engl et al., 1996; Osher et al., 2005), and it is further generalized to the 2nd order (\(k=2\)) (Koenker and Mizera, 2004; Hinterberger and Scherzer, 2006) and general order (\(k\in\mathbb{N}\)) (Bredies and Holler, 2020). Then, we may train the regression model \(f_{\theta}\) by minimizing the loss function equipped with the variation regularization:

\[L_{A,q}(\theta):=\underbrace{A(\theta)}_{\text{(loss func.)}}+\underbrace{ \lambda B(\theta)}_{\text{(param. neg.)}}+\underbrace{\sum_{k=0}^{K}\eta_{k}C_{k,q}(f_{ \theta})}_{\text{(variation reg.)}} \tag{3}\]

with user-specified hyperparameters \(\lambda\geq 0\) and \(\eta=(\eta_{0},\eta_{1},\eta_{2},\ldots,\eta_{K})\in\mathbb{R}_{\geq 0}^{K+1}\).

One important problem here is the computational intractability of the \((k,q)\)-VR term (2). Generally speaking, the explicit form of the integral term with highly non-linear models (such as single-hidden-layer perceptron (1) considered in this study) cannot be obtained. Although we may compute numerical integration for evaluating (2), it requires large number of samples for conducting numerical integration with high precision; numerical integration is incompatible with iterative algorithms such as gradient descent.

While the numerical integration is needed to compute the full-batch gradient descent, theories on stochastic optimization originated from Robbins and Monro (1951) suggest that only the unbiased estimate of the full-batch gradient is needed to optimize the loss function. Following a similar idea in the optimization of intractable likelihood (Geyer and Thompson, 1992), contrastive divergence (Carreira-Perpinan and Hinton, 2005), and robust divergence (Okuno, 2023), this study proposes a stochastic gradient descent that can be computed efficiently without conducting numerical integration.

### A stochastic gradient descent

To minimize the loss function (3) defined with the \((k,q)\)-VR (2), this study proposes a stochastic gradient descent algorithm described herein.

Firstly, we define a stochastic gradient, which is an unbiased estimate of the full-batch gradient. At the iteration \(t\in\mathbb{N}\), pick \((x_{1}^{(t)},y_{1}^{(t)}),(x_{2}^{(t)},y_{2}^{(t)}),\ldots,(x_{n}^{(t)},y_{n}^ {(t)})\) from \(\{(x_{i},y_{i})\}_{i=1}^{N}\) uniformly randomly; \(\{(x_{i}^{(t)},y_{i}^{(t)})\}\) is used for computing the gradient of the term \(A(\theta)\). Also pick \(z_{1}^{(t)},z_{2}^{(t)},\ldots,z_{m}^{(t)}\) uniformly randomly from \(\Omega=[D_{1},D_{2}]\); the \(\{z_{j}^{(t)}\}\) is used for unbiasedly estimating the integral-based \((k,q)\)-VR term (2). Specifically, we
define three terms

\[\alpha^{(t,n)}(\theta) :=-\frac{1}{n}\sum_{i=1}^{n}\nu^{[1]}(y_{i}^{(t)}-f_{\theta}(x_{i}^{( t)}))\frac{\partial}{\partial\theta}f_{\theta}(x_{i}^{(t)}), \tag{4}\] \[\beta(\theta) :=\frac{\partial}{\partial\theta}B(\theta),\] (5) \[\gamma_{k,q}^{(t,m)}(f_{\theta}) :=\frac{D_{2}-D_{1}}{m}\sum_{j=1}^{m}\text{sign}(f_{\theta}^{[k]}( z_{j}^{(t)}))|f_{\theta}^{[k]}(z_{j}^{(t)})|^{q-1}\frac{\partial}{\partial \theta}f_{\theta}^{[k]}(z_{j}^{(t)}), \tag{6}\]

where \(\nu^{[1]}\) is the 1st order derivative of the function \(\nu\) and \(\text{sign}(z)=\mathbb{1}\left(z>0\right)-\mathbb{1}\left(z<0\right)\). Then, the terms \(\alpha^{(t,m)}(\theta),\beta(\theta),\gamma_{k,q}^{(t,m)}(f_{\theta})\) are unbiased estimators of the gradients \(\partial A(\theta)/\partial\theta,\partial B(\theta)/\partial\theta, \partial C_{k}(f_{\theta})/\partial\theta\), respectively. See Appendix B for the explicit forms of the gradients of the \(k\)th order derivative \(f_{\theta}^{[k]}(x)\) for the single-hidden-layer perceptron (1). Note that the gradient also can be computed by an automatic differentiation implemented for training deep neural network (see, e.g., Paszke et al. (2017)). Then, it holds for a stochastic gradient

\[g_{\lambda,\eta}^{(t,n,m)}(\theta):=\alpha^{(t,n)}(\theta)+\lambda\beta( \theta)+\sum_{k=0}^{K}\eta_{k}Y_{k,q}^{(t,m)}(f_{\theta}) \tag{7}\]

that

\[\mathbb{E}^{(t)}\left(g_{\lambda,\eta}^{(t,n,m)}(\theta)\right)=\frac{ \partial}{\partial\theta}A(\theta)+\lambda\frac{\partial}{\partial\theta}B( \theta)+\sum_{k=0}^{K}\eta_{k}\frac{\partial}{\partial\theta}C_{k,q}(f_{ \theta})=\frac{\partial}{\partial\theta}L_{\lambda,\eta}(\theta). \tag{8}\]

\(\mathbb{E}^{(t)}\) represents the expectation with respect to \(X^{(t)}=\{x_{i}^{(t)}\}\), \(Y^{(t)}=\{y_{i}^{(t)}\}\), \(Z^{(t)}=\{z_{j}^{(t)}\}\). Note that the unbiased-ness (8) holds for any \(n,m\in\mathbb{N}\).

Using this stochastic gradient (7), we employ a stochastic gradient descent

\[\theta^{(t)}=\theta^{(t-1)}-\omega_{t}g_{\lambda,\eta}^{(t-1,n,m)}(\theta^{( t-1)}),\quad(t=1,2,\ldots,T). \tag{9}\]

The stochastic gradient descent (9) equipped with the unbiased estimator of the gradient (7) and the decreasing learning rate \(\omega_{t}\searrow 0\) is proved to optimize the loss function \(L_{\lambda,\eta}(\theta)\). See Proposition 1 for more rigorous descriptions showing a sufficient condition of the convergence.

**Proposition 1** (A simpler version of Ghadimi and Lan (2013) Theorem 2.1 (a)).: Let \(n,m\in\mathbb{N}\) be arbitrarily fixed. Assume that (i) \(f(\theta)=L_{\lambda,\eta}(\theta)\) is smooth, (ii) gradient of \(f(\theta)\) satisfies the Lipschitz property, i.e., \(\|\partial f(\theta)/\partial\theta-\partial f(\theta^{\prime})/\partial \theta\|\leq L\|\theta-\theta^{\prime}\|\), (iii) \(\mathbb{E}^{(t)}(g_{\lambda,\eta}^{(t,n,m)}(\theta^{(t)}))=\partial f(\theta^{ (t)})/\partial\theta\), and (iv) \(\mathbb{E}^{(t)}(\|g_{\lambda,\eta}^{(t,n,m)}(\theta^{(t)})-\partial f(\theta^ {(t)})/\partial\theta\|^{2})\leq\sigma^{2}\) for some \(\sigma\geq 0\), for any \(t\in 1,2,\ldots,T\). Assume that the learning rate \(\{\omega_{t}\}_{t=1}^{T}\) satisfies \(\sum_{t=1}^{T}\omega_{t}-\infty\) and \(\{\sum_{t=1}^{T}\omega_{t}^{-1}\sum_{t=1}^{T}\omega_{t}^{2}\to 0\) as \(T\to\infty\). Then, the sequence \(\{\theta^{(t)}\}\) obtained by the stochastic gradient descent (9) satisfies

\[\mathbb{E}_{\tau}\left(\left\|\frac{\partial}{\partial\theta}L_{\lambda,\eta}( \theta)\right\|^{2}\right)\overset{\text{in prob.}}{\rightharpoonup}0\quad(T \to\infty).\]

\(\mathbb{E}_{\tau}\) represents the expectation with respect to the step \(\tau\in\{1,2,\ldots,T\}\) randomly chosen with the probability \(\mathbb{P}(\tau=k\mid T)=\{2\omega_{k}-L\omega_{k}^{2}\}/\sum_{k=1}^{T}\{2 \omega_{k}-L\omega_{k}^{2}\}\) (\(k=1,2,\ldots,T\)).

Proposition 1 proves that the parameter estimated by the proposed stochastic algorithm is expected to reach to a (local) minima of the exact loss function \(L_{\lambda,\eta}(\theta)\), but not the finite approximation of the loss function via numerical integration. It is also noted that the convergence holds regardless of the sample sizes \(n,m\in\mathbb{N}\) (namely, there is no need to take the limit \(n,m\to\infty\) to guarantee the convergence; \(n=m=5\) is used in our numerical experiments). As stochastic algorithm includes randomness to compute the stochastic gradient, it is believed that stochastic gradient algorithms are more likely to escape from the local minima or a saddle point (compared with the full-batch gradient descent), i.e., the stochastic solution would be better than the fullbatch-based estimators for some cases. See, e.g., Jin et al. (2021).
\begin{table}
\begin{tabular}{l c c c} \hline  & \(f_{1}(x)\): linear & \(f_{2}(x)\): quadratic & \(f_{3}(x)\): cubic \\ \hline No regularization & \(0.943\pm 0.0092\) & \(0.938\pm 0.0106\) & \(0.863\pm 0.0271\) \\ Smaller NN (\(L=20\)) & \(0.971\pm 0.0095\) & \(0.916\pm 0.0289\) & \(0.689\pm 0.1890\) \\ \(L_{2}\) (\(\lambda=10^{-1}\)) & \(0.975\pm 0.0054\) & \(0.971\pm 0.0069\) & \(0.928\pm 0.0234\) \\ \(L_{2}\) (\(\lambda=10^{-3}\)) & \(0.933\pm 0.0121\) & \(0.909\pm 0.0248\) & \(0.782\pm 0.0443\) \\ \(L_{2}\) (\(\lambda=10^{-5}\)) & \(0.960\pm 0.0114\) & \(0.930\pm 0.0211\) & \(0.849\pm 0.0528\) \\ \hline (1,2)-VR (\(\eta_{1}=10^{-1}\)) & \(0.996\pm 0.0007\) & \(0.967\pm 0.0271\) & \(0.823\pm 0.1720\) \\ (2,2)-VR (\(\eta_{2}=10^{-3}\)) & \(0.999\pm 0.0010\) & \(0.994\pm 0.0031\) & \(0.842\pm 0.1000\) \\ (3,2)-VR (\(\eta_{3}=10^{-5}\)) & \(0.999\pm 0.0007\) & \(0.996\pm 0.0036\) & \(0.966\pm 0.0301\) \\ \hline \end{tabular}
\end{table}
Table 2: Predictive correlation for \(N=10^{2}\).

\begin{table}
\begin{tabular}{l c c c} \hline  & \(f_{1}(x)\): linear & \(f_{2}(x)\): quadratic & \(f_{3}(x)\): cubic \\ \hline No regularization & \(0.943\pm 0.0092\) & \(0.938\pm 0.0106\) & \(0.863\pm 0.0271\) \\ Smaller NN (\(L=20\)) & \(0.971\pm 0.0095\) & \(0.916\pm 0.0289\) & \(0.689\pm 0.1890\) \\ \(L_{2}\) (\(\lambda=10^{-1}\)) & \(0.975\pm 0.0054\) & \(0.971\pm 0.0069\) & \(0.928\pm 0.0234\) \\ \(L_{2}\) (\(\lambda=10^{-3}\)) & \(0.933\pm 0.0121\) & \(0.909\pm 0.0248\) & \(0.782\pm 0.0443\) \\ \(L_{2}\) (\(\lambda=10^{-5}\)) & \(0.960\pm 0.0114\) & \(0.930\pm 0.0211\) & \(0.849\pm 0.0528\) \\ \hline (1,2)-VR (\(\eta_{1}=10^{-1}\)) & \(0.996\pm 0.0007\) & \(0.967\pm 0.0271\) & \(0.823\pm 0.1720\) \\ (2,2)-VR (\(\eta_{2}=10^{-3}\)) & \(0.999\pm 0.0010\) & \(0.994\pm 0.0031\) & \(0.842\pm 0.1000\) \\ (3,2)-VR (\(\eta_{3}=10^{-5}\)) & \(0.999\pm 0.0007\) & \(0.996\pm 0.0036\) & \(0.966\pm 0.0301\) \\ \hline \end{tabular}
\end{table}
Table 3: Predictive correlation for \(N=10^{3}\).
Figure 2: Neural networks trained with \(N=100\) samples.
Experiments

This section conducts numerical experiments to demonstrate the proposed approach. Source codes to reproduce numerical results are provided in [https://github.com/oknakfm/HOVR](https://github.com/oknakfm/HOVR).

### Experimental settings

Dataset generation:we generate \(x_{i}\sim U(-1/2,1/2)\) uniformly randomly from the interval \(\Omega=[-1/2,1/2]\) and normal random numbers \(y_{i}\sim N(f(x_{i}),\sigma^{2})\) with \(\sigma=0.1\). We consider three settings: (i) linear: \(f(x)=x\), (ii) quadratic: \(f(x)=4x^{2}\), (iii) cubic: \(f(x)=\frac{6q}{7}(x+3/8)x(x-3/8)\).

Neural network architecture:we employ a single-hidden-layer perceptron (1) with \(L=200\) hidden units. Hyperbolic tangent function \(\tanh(z)=\{\exp(z)-\exp(-z)\}/\{\exp(z)+\exp(-z)\}\) is employed for the activation function \(\sigma(z)\).

Neural network training:we compute a stochastic gradient descent described in Section 2.3. We employ \(\nu(z)=z^{2}\) for the loss function \(A(\theta)\) and \(q=2\) for the variation regularization. We randomly pick \(n=m=5\) samples for computing the stochastic gradient \(\alpha^{(t,m)}(\theta)\) and \(\gamma^{(t,m)}_{k}(f_{0})\). The parameters of the neural network (1) is initialized by normal random numbers: we randomly pick \(a^{(0)}_{\ell},c^{(0)}_{\ell}\) from \(N(0,1)\), \(b^{(0)}_{\ell}\) from \(N(0,1/\hat{\psi}(\{x_{i}\}))\), and specify \(d^{(0)}=0\). The learning rate is designed to be cyclic (by following the similar idea as Smith (2017)): the learning rate is initialized by \(\omega_{0}=10^{-3}\), and the rate is multiplied by \(0.9\) for each \(25\) iterations. For each \(10^{3}\) iteration, the learning rate is pull-back to \(10^{-3}\). Overall, we compute \(T=2\times 10^{4}\) iterations, meaning that our SGD repeats \(20\) cycles of decreasing learning rate.

Baselines:we train the neural network (i) without regularization, (ii) with a smaller number of hidden units (\(L=20\)), and (iii) with \(L_{2}\) regularization as baselines.

Evaluation metric:after training the neural network using the proposed approach and baselines, we compute the predicted \(\hat{y}_{i}\) for \(x_{i}\) (\(i=1,2,\ldots,10^{4}\)) that are regularly placed over the interval \([-1/2+0.1,1/2-0.1]\). Using \(y^{*}_{i}=f(x_{i})\), we compute the (predictive) correlation coefficient. Larger scores are better. We compute the predictive correlation \(5\) times with different random seeds. We summarize the results with the mean and the standard deviation for comparison purposes.

### Results

Experimental results are shown in Tables 1-3. Also see Figure 2 for the illustration of the neural networks trained with \(N=10^{2}\) samples.

Overall, variation regularization demonstrates the highest predictive correlation for all the settings (\(N=10^{1}\), \(10^{2}\), \(10^{3}\) and \(f_{1},f_{2},f_{3}\)) while all the remaining neural networks are trained with the same optimization setting (except for the type of regularization). There are several observations through these experiments.

First, regularizing higher-order variation yields "smoother" functions. See Figure 2: regularizing the \(3\)rd order derivative (\(\eta_{3}=10^{-5}\)) yields smoother function than those with \(2\)nd (\(\eta_{2}=10^{-3}\)) and \(1\)st (\(\eta_{1}=10^{-1}\)) order derivatives. Particularly, regularizing the \(1\)st order derivative (\(\eta_{1}=10^{-1}\)) forces the estimated functions to be closer to piece-wise constant functions, and regularizing the \(2\)nd order derivative (\(\eta_{2}=10^{-3}\) forces the functions to be closer to piece-wise linear ones. Second, the \(L_{2}\) regularization provides functions that still contain many inflection points while \((k,q)\)-VR provides smoother functions. Third, reducing the number of hidden units in NN is effective when the underlying function is simple. For instance, the prediction accuracy is improved by reducing the number of hidden units \(L\) if the underlying function is linear. However, if the underlying function is more complicated, for instance, quadratic and cubic, the prediction accuracy is degraded as the smaller NN has less expressive power. The smaller NN is not capable of fully approximating the non-linear functions.
Application to the physics-informed training of a neural network

Besides traditional numerical approaches such as finite element method (see, e.g., Oden and Reddy (2012)), neural networks have been leveraged to solve a user-specified differential equations (Dissanayake and Phan-Thien, 1994; Lagaris et al., 1998) so that the neural network model itself is trained to satisfy the differential equations. In recent years, this approach has been combined with automatic differentiation widely-used in modern deep neural network frameworks (Paszke et al., 2017). The neural network framework that aims to solve the differential equation with the aid of automatic differentiation is especially called physics-informed training of a neural network (or simply referred to as physics-informed neural network; PINN), and PINN has drawn considerable attention these days (Berg and Nystrom, 2018; Cuomo et al., 2022), especially in the field of natural science. See a comprehensive survey (Cuomo et al., 2022) for more details; herein, we summarize the concept of PINN and the relation to this study. While we describe the ordinary differential equation (ODE) for simplicity, the concept has been extended to partial differential equations.

Let \(D_{1},D_{2}\in\mathbb{N}\), and let \(\Omega:=[D_{1},D_{2}]\) be the support of the function to be considered. Let \(f_{\theta}\) be a neural network and let \(f_{\theta}^{[k]}\) denote its \(k\)th derivative (with respect to the input \(x\)). Let \(\mathcal{L}[f](x)=0\) (\(\forall x\in\Omega\)) be the ODE to be solved (for instance, \(\mathcal{L}[f]=f^{[2]}+\alpha f\) represents a simple harmonic motion). Then, PINN is formulated to minimize a loss function

\[\int_{\Omega}\mathcal{L}[f_{\theta}](x)^{2}\mathrm{d}x+\kappa \mathcal{B}[f_{\theta}], \tag{10}\]

where \(\mathcal{B}[f_{\theta}]\) represents the loss function for boundary conditions (e.g., \(\mathcal{B}[f_{\theta}]=\{f_{\theta}(D_{1})-f_{1}\}^{2}+\{f_{\theta}(D_{2})- f_{2}\}^{2}\) for conditions \(f(D_{1})=f_{1},f(D_{2})=f_{2}\)) and \(\kappa>0\) is a hyperparameter. Here, we can further consider the case that several labels for the underlying function \(f\) is available, i.e., we assume that the pairs \((x_{i},y_{i})\) satisfying \(y_{i}\approx f(x_{i})\) are available, we can introduce the label loss into the original form (10):

\[\frac{1}{N}\sum_{i=1}^{N}\nu(y_{i}-f_{\theta}(x_{i}))+\frac{\eta }{2}\int_{\Omega}\mathcal{L}[f_{\theta}](x)^{2}\mathrm{d}x+\kappa\mathcal{B} [f_{\theta}]. \tag{11}\]

While the original PINN is designed to solve the ODE, minimizing (11) also can be regarded as the regression problem with the ODE-based regularization (with not that large hyperparameter \(\eta>0\)). The \((k,q)\)-VR term \(C_{k,q}(f_{\theta})\) coincides with the second term in the ODE-constrained regression (11) by specifying \(\mathcal{L}[f_{\theta}]=|f_{\theta}^{[k]}|^{2/q}\), so the training of neural network with the \((k,q)\)-VR is closely related to the PINN formulation.

While there is a similarity, our approach is significantly different from the PINN framework, from the computational perspective. In the computation of PINN, explicit form of the integral term \(\int_{\Omega}\mathcal{L}[f_{\theta}](x)^{2}\mathrm{d}x\) is hard to be obtained (as well as the computational difficulty of the \((k,q)\)-VR term discussed so far), thus existing studies finitely approximate this integration as

\[\int_{\Omega}\mathcal{L}[f_{\theta}](x)^{2}\mathrm{d}x\approx \frac{D_{2}-D_{1}}{M}\sum_{j=1}^{M}\mathcal{L}[f_{\theta}](\bar{x}_{j})^{2} \tag{12}\]

with \(M\) training points \(\{\bar{x}_{j}\}_{j=1}^{M}\subset\Omega\) regularly picked from the interval \(\Omega=[D_{1},D_{2}]\). See, e.g., Section 2.3 in Cuomo et al. (2022) for the current surrounding situation in the computation of PINN, including the efficient approximation of the integral (12).

Due to the finite approximation, theoretical analyses indicate that the (practically-computed) finite approximation of PINN includes an approximation error depending on the parameter \(M\in\mathbb{N}\) adjusting the precision (Shin et al., 2020; De Ryck et al., 2021). See Section 2.4 in Cuomo et al. (2022) for more details. As long as the training points \(\{\bar{x}_{j}\}_{j=1}^{M}\) are first fixed, the approximation error does not vanish even if the finite approximation is minimized by the stochastic gradient descent.

Our proposed stochastic algorithm can be straightforwardly generalized to this PINN setting. See (13) in Appendix C; the PINN extension of our stochastic gradient is an unbiased estimator of the gradient for the original PINN loss function (11) which contains the integral term. Therefore, the stochastic gradient descent using the unbiased gradient (13) directly minimizes the loss function (11) without going through the finite approximation (12). Although the PINN extension goes beyond the main scope of this study, we certainly admit a huge potential for this PINN direction.
Conclusion and possible future direction

This study considered directly regularizing \(k\)th order derivative of the parametric models to be trained. While the \((k,q)\)-VR is represented by a computationally-intractable integral, the loss function equipped with the \((k,q)\)-VR can be efficiently minimized by a proposed stochastic optimization algorithm. Numerical experiments demonstrated the training of highly-nonlinear neural networks with the \((k,q)\)-VR terms. Although this study considered only the regression problems for simplicity, \((k,q)\)-VR terms can be straightforwardly incorporated into various statistical problems, such as classification. Furthermore, the proposed approach can be applied to general parametric models, including the arbitrary structure of deep neural networks.

The possible future direction of the proposed approach is diverse. One direction is the application to physics-informed neural networks described in Section 4. The proposed approach is free from the approximation error which is originated from the finite approximation of the integral. Another direction is the fusion of expressive neural networks and robust estimation (against outliers). While expressive models such as neural networks may fit both the underlying target function and the outlier distribution, regularization of the higher-order variation may prevent "overfitting" to the outlier distribution.

## Acknowledgement

A. Okuno was supported by JSPS KAKENHI (21K17718, 22H05106). We would like to thank Keisuke Yano for the helpful discussions.

## Appendix A Regularization for additive model

An additive model \(f_{\theta}^{(\text{add.})}(x)=\sum_{j=1}^{J}\theta_{j}\phi_{j}(x)\) with user-specified basis functions \(\{\phi_{j}\}\) encompasses kernel and spline regression models as special cases. With the gram matrix \(G=(g_{ij}),g_{ij}=\int_{\Omega}\phi_{i}(x)\phi_{j}(x)\mathrm{d}x\), we have

\[\|f_{\theta}^{(\text{add.})}\|_{L^{2}(\Omega)}^{2}=\int_{\Omega}f_{\theta}^{( \text{add.})}(x)^{2}\mathrm{d}x=\sum_{l=1}^{J}\sum_{j=1}^{J}\theta_{i}\theta_ {j}\underbrace{\int_{\Omega}\phi_{l}(x)\phi_{j}(x)\mathrm{d}x}_{-g_{ij}}=\sum_ {i=1}^{J}\theta_{i}\theta_{j}g_{ij}=\langle\theta,G\theta\rangle.\]

See, e.g., Smale and Zhou (2007) for more details of the regularization for kernel regression.

## Appendix B Gradient of \(k\)th order variation

The \(k\)th order variation of the single-hidden layer perceptron (1) is obtained as

\[\frac{\partial^{k}}{\partial x^{k}}f_{\theta}(x)=\sum_{\ell=1}^{L}\alpha_{ \ell}b_{\ell}^{k}\sigma^{[k]}(b_{\ell}x+c_{\ell})+d\mathbb{1}_{[k=0]},\]

where \(\sigma^{[k]}\) denotes the \(k\)th derivative of the activation function and \(\mathbb{1}_{\mathcal{A}}\) denotes the indicator function which outputs 1 if and only if the event \(\mathcal{A}\) occurs, and 0 otherwise. For the hyperbolic tangent activation function \(\sigma(z)=\tanh(z)=\{\exp(z)-\exp(-z)\}/\{\exp(z)+\exp(-z)\}\), we have \(\sigma^{[1]}(z)=1-\tanh^{2}z\), \(\sigma^{[2]}(z)=-2\tanh z(1-\tanh^{2}z)\), \(\sigma^{[3]}(z)=2(1-\tanh^{2}z)(3\tanh^{2}z-1)\), and \(\sigma^{[4]}(z)=8\tanh(z)(\tanh^{2}z-1)(3\tanh^{2}z-2)\). Then, the
gradient of the \(k\)th order variation is also obtained as

\[\frac{\partial}{\partial a_{\ell^{\prime}}}\frac{\partial^{k}}{ \partial x^{k}}f_{\theta}(x) =b_{\ell^{\prime}}^{k}\sigma^{|k|}(b_{\ell^{\prime}}x+c_{\ell^{ \prime}}),\] \[\frac{\partial}{\partial b_{\ell^{\prime}}}\frac{\partial^{k}}{ \partial x^{k}}f_{\theta}(x) =ka_{\ell^{\prime}}b_{\ell^{\prime}}^{k-1}\sigma^{|k|}(b_{\ell^{ \prime}}x+c_{\ell^{\prime}})+a_{\ell^{\prime}}b_{\ell^{\prime}}^{k}x\sigma^{|k |+1|}(b_{\ell^{\prime}}x+c_{\ell^{\prime}})\] \[\frac{\partial}{\partial c_{\ell^{\prime}}}\frac{\partial^{k}}{ \partial x^{k}}f_{\theta}(x) =a_{\ell^{\prime}}b_{\ell^{\prime}}^{k}\sigma^{|k|+1|}(b_{\ell^{ \prime}}x+c_{\ell^{\prime}}),\] \[\frac{\partial}{\partial d}\frac{\partial^{k}}{\partial x^{k}}f_{ \theta}(x) =\mathbb{1}_{\{k=0\}}.\]

## Appendix C An unbiased stochastic gradient of the PINN loss function

Pick \(\{(x_{i}^{(t)},y_{i}^{(t)})\}_{i=1}^{n}\) from \(\{(x_{i},y_{i})\}_{i=1}^{N}\), and also pick \(\{z_{j}^{(t)}\}_{j=1}^{m}\) from the interval \(\Omega=\{D_{1},D_{2}\}\) uniformly randomly. Then, we define

\[-\frac{1}{n}\sum_{i=1}^{n}\nu^{|1|}(y_{i}^{(t)}-f_{\theta}(x_{i}^{(t)}))\frac{ \partial}{\partial\theta}f_{\theta}(x_{i}^{(t)})+\eta\frac{D_{2}-D_{1}}{m} \sum_{j=1}^{m}\mathcal{L}[f_{\theta}](z_{j}^{(t)})\frac{\partial}{\partial \theta}\mathcal{L}[f_{\theta}](z_{j}^{(t)})+\kappa\frac{\partial}{\partial \theta}\mathcal{B}[f_{\theta}]. \tag{13}\]

(13) is an unbiased estimator of the gradient of the PINN loss (11) regardless of the number of samples \(n,m\in\mathbb{N}\). The same discussion shown in Section 2.3 can be applied to the stochastic gradient descent using (13).
