Toward Robust Uncertainty Estimation with Random Activation Functions

TL;DR

This paper introduces Random Activation Functions (RAFs) Ensemble, a novel method that enhances uncertainty estimation in neural networks by diversifying models with different random activation functions, leading to improved robustness in out-of-distribution scenarios.

Contribution

The paper proposes RAFs Ensemble, a new ensemble approach that increases diversity and robustness in uncertainty quantification by assigning different random activation functions to each network.

Findings

Outperforms state-of-the-art ensemble uncertainty methods

Demonstrates robustness on synthetic and real-world regression datasets

Improves uncertainty estimation in out-of-distribution cases

Abstract

Deep neural networks are in the limelight of machine learning with their excellent performance in many data-driven applications. However, they can lead to inaccurate predictions when queried in out-of-distribution data points, which can have detrimental effects especially in sensitive domains, such as healthcare and transportation, where erroneous predictions can be very costly and/or dangerous. Subsequently, quantifying the uncertainty of the output of a neural network is often leveraged to evaluate the confidence of its predictions, and ensemble models have proved to be effective in measuring the uncertainty by utilizing the variance of predictions over a pool of models. In this paper, we propose a novel approach for uncertainty quantification via ensembles, called Random Activation Functions (RAFs) Ensemble, that aims at improving the ensemble diversity toward a more robust estimation, by accommodating each neural network with a different (random) activation function. Extensive empirical study demonstrates that RAFs Ensemble outperforms state-of-the-art ensemble uncertainty quantification methods on both synthetic and real-world datasets in a series of regression tasks.

Tables

Table 1. Table 1: Performance of methods on all datasets w.r.t. NLL, including 95% confidence intervals. The best scores are in bold.

			NLL
	DE	HDE	AE	NTKGP-p.	RP-p.	RAFs
He et al. 1D	$>$100 $\pm$ 0.18	71.31 $\pm$ 0.51	38.75 $\pm$ 0.12	4.48 $\pm$ 0.18	13.05 $\pm$ 0.43	2.21 $\pm$ 0.18
Forrester et al. 1D	$>$100 $\pm$ 0.53	$>$100 $\pm$ 0.51	50.82 $\pm$ 0.52	$>$100 $\pm$ 0.50	13.7 $\pm$ 0.58	0.64 $\pm$ 0.74
Schaffer N.4 2D	0.29 $\pm$ 0.01	-0.71 $\pm$ 0.01	2.15 $\pm$ 0.01	-0.55 $\pm$ 0.01	-0.36 $\pm$ 0.01	-0.79 $\pm$ 0.01
Double pendulum 2D	2.95 $\pm$ 0.05	2.18 $\pm$ 0.84	-0.36 $\pm$ 0.05	-0.58 $\pm$ 0.05	-0.47 $\pm$ 0.05	-0.49 $\pm$ 0.04
Rastrigin 3D	29.24 $\pm$ 1.30	3.09 $\pm$ 1.15	35.94 $\pm$ 0.74	28.38 $\pm$ 0.64	4.35 $\pm$ 1.24	3.44 $\pm$ 1.05
Ishigami 3D	6.01 $\pm$ 0.08	$>$100 $\pm$ 0.08	8.73 $\pm$ 0.08	1.53 $\pm$ 0.08	-0.01 $\pm$ 0.08	0.06 $\pm$ 0.07
Environmental 4D	64.72 $\pm$ 0.23	7.84 $\pm$ 0.13	1.65 $\pm$ 0.20	4.5 $\pm$ 0.27	3.94 $\pm$ 0.21	0.81 $\pm$ 0.17
Griewank 4D	28.29 $\pm$ 2.43	5.50 $\pm$ 1.62	4.64 $\pm$ 3.06	10.21 $\pm$ 2.37	4.29 $\pm$ 2.93	4.79 $\pm$ 2.40
Roos & Arnold 5D	-2.02 $\pm$ 0.01	-2.21 $\pm$ 0.00	-1.89 $\pm$ 0.01	-1.71 $\pm$ 0.01	-1.70 $\pm$ 0.01	-2.1 $\pm$ 0.01
Friedman 5D	96.94 $\pm$ 0.41	$>$100 $\pm$ 0.51	15.04 $\pm$ 0.50	41.69 $\pm$ 0.39	4.22 $\pm$ 0.44	1.78 $\pm$ 0.39
Planar arm torque 6D	9.58 $\pm$ 0.07	4.11 $\pm$ 0.08	3.07 $\pm$ 0.05	-0.32 $\pm$ 0.08	-0.05 $\pm$ 0.07	-0.16 $\pm$ 0.06
Sum of powers 6D	$>$100 $\pm$ 0.41	$>$100 $\pm$ 0.62	55.03 $\pm$ 0.43	$>$100 $\pm$ 0.41	41.59 $\pm$ 0.40	35.22 $\pm$ 0.35
Ackley 7D	7.11 $\pm$ 0.23	1.38 $\pm$ 0.16	2.50 $\pm$ 0.36	3.11 $\pm$ 0.27	2.09 $\pm$ 0.26	1.16 $\pm$ 0.08
Piston simulation 7D	-2.19 $\pm$ 0.00	14.06 $\pm$ 0.00	3.50 $\pm$ 2.40	2.87 $\pm$ 2.93	2.67 $\pm$ 0.42	3.63 $\pm$ 0.57
Robot arm 8D	10.71 $\pm$ 0.03	6.87 $\pm$ 0.01	7.11 $\pm$ 0.01	0.27 $\pm$ 0.03	0.80 $\pm$ 0.06	0.25 $\pm$ 0.02
Borehole 8D	$>$100 $\pm$ 1.01	$>$100 $\pm$ 1.01	4.89 $\pm$ 1.87	5.48 $\pm$ 3.54	4.06 $\pm$ 1.20	4.36 $\pm$ 1.26
Styblinski-Tang 9D	$>$100 $\pm$ 3.05	$>$100 $\pm$ 0.00	40.80 $\pm$ 5.33	$>$100 $\pm$ 3.03	15.82 $\pm$ 6.31	25.23 $\pm$ 4.12
PUMA560 9D	6.59 $\pm$ 0.15	1.62 $\pm$ 0.14	4.24 $\pm$ 0.14	5.93 $\pm$ 0.08	6.40 $\pm$ 0.14	2.14 $\pm$ 0.13
Adapted Welch 10D	$>$100 $\pm$ 0.81	$>$100 $\pm$ 0.75	$>$100 $\pm$ 0.55	$>$100 $\pm$ 0.75	$>$100 $\pm$ 0.57	78.53 $\pm$ 0.67
Wing weight 10D	$>$100 $\pm$ 0.00	27.31 $\pm$ 4.37	5.46 $\pm$ 4.36	67.30 $\pm$ 0.53	5.54 $\pm$ 4.15	5.39 $\pm$ 1.69
Boston housing	74.54 $\pm$ 1.06	$>$100 $\pm$ 1.04	71.53 $\pm$ 1.06	70.82 $\pm$ 1.06	$>$100 $\pm$ 1.10	40.67 $\pm$ 1.00
Abalone	$>$100 $\pm$ 0.10	$>$100 $\pm$ 0.10	47.67 $\pm$ 0.10	$>$100 $\pm$ 0.10	28.37 $\pm$ 0.10	28.90 $\pm$ 0.10
Naval propulsion	-2.27 $\pm$ 0.00	$>$100 $\pm$ 0.00	3.92 $\pm$ 0.10	2.28 $\pm$ 1.51	2.16 $\pm$ 0.16	1.91 $\pm$ 0.07
Forest fire	15.71 $\pm$ 0.05	3.14 $\pm$ 0.02	2.66 $\pm$ 0.69	3.10 $\pm$ 1.11	4.68 $\pm$ 0.14	2.15 $\pm$ 0.28
Parkinson’s	26.74 $\pm$ 0.02	$>$100 $\pm$ 0.10	$>$100 $\pm$ 0.16	$>$100 $\pm$ 0.03	$>$100 $\pm$ 0.16	45.69 $\pm$ 0.16

Table 2. Table 2: Performance of methods on all datasets w.r.t. RMSE, including 95% confidence intervals. The best scores are in bold.

			RMSE
	DE	HDE	AE	NTKGP-p.	RP-p.	RAFs
He et al. 1D	3.71 $\pm$ 0.18	5.70 $\pm$ 0.51	3.15 $\pm$ 0.12	3.64 $\pm$ 0.18	5.24 $\pm$ 0.43	3.80 $\pm$ 0.18
Forrester et al. 1D	5.00 $\pm$ 0.53	4.12 $\pm$ 0.51	4.09 $\pm$ 0.52	6.05 $\pm$ 0.50	5.70 $\pm$ 0.58	2.8 $\pm$ 0.74
Schaffer N.4 2D	0.23 $\pm$ 0.01	0.34 $\pm$ 0.01	0.30 $\pm$ 0.01	0.24 $\pm$ 0.01	0.31 $\pm$ 0.01	0.27 $\pm$ 0.01
Double pendulum 2D	0.46 $\pm$ 0.05	2.22 $\pm$ 0.84	0.71 $\pm$ 0.05	0.51 $\pm$ 0.05	0.74 $\pm$ 0.05	0.58 $\pm$ 0.04
Rastrigin 3D	18.41 $\pm$ 1.30	10.96 $\pm$ 1.15	25.58 $\pm$ 0.74	18.10 $\pm$ 0.64	12.87 $\pm$ 1.24	14.85 $\pm$ 1.05
Ishigami 3D	0.69 $\pm$ 0.08	1.05 $\pm$ 0.08	0.88 $\pm$ 0.08	0.69 $\pm$ 0.08	0.58 $\pm$ 0.08	0.57 $\pm$ 0.07
Environmental 4D	2.04 $\pm$ 0.23	2.51 $\pm$ 0.13	1.83 $\pm$ 0.20	2.34 $\pm$ 0.27	2.03 $\pm$ 0.21	1.68 $\pm$ 0.17
Griewank 4D	83.97 $\pm$ 2.43	45.68 $\pm$ 1.62	42.12 $\pm$ 3.06	78.47 $\pm$ 2.37	38.62 $\pm$ 2.93	78.79 $\pm$ 2.40
Roos & Arnold 5D	0.07 $\pm$ 0.01	0.01 $\pm$ 0.00	0.07 $\pm$ 0.01	0.09 $\pm$ 0.01	0.08 $\pm$ 0.01	0.08 $\pm$ 0.01
Friedman 5D	3.17 $\pm$ 0.41	3.63 $\pm$ 0.51	2.95 $\pm$ 0.50	3.39 $\pm$ 0.39	2.74 $\pm$ 0.44	3.1 $\pm$ 0.39
Planar arm torque 6D	0.65 $\pm$ 0.07	0.62 $\pm$ 0.08	0.71 $\pm$ 0.05	0.71 $\pm$ 0.08	1.08 $\pm$ 0.07	0.74 $\pm$ 0.06
Sum of powers 6D	22.81 $\pm$ 0.41	21.19 $\pm$ 0.62	21.87 $\pm$ 0.43	22.79 $\pm$ 0.41	22.22 $\pm$ 0.40	22.24 $\pm$ 0.35
Ackley 7D	8.92 $\pm$ 0.23	2.43 $\pm$ 0.16	7.28 $\pm$ 0.36	8.58 $\pm$ 0.27	4.03 $\pm$ 0.26	1.33 $\pm$ 0.08
Piston simulation 7D	0.02 $\pm$ 0.00	0.04 $\pm$ 0.00	29.1 $\pm$ 2.40	$>$100 $\pm$ 2.93	5.78 $\pm$ 0.42	7.40 $\pm$ 0.57
Robot arm 8D	0.92 $\pm$ 0.03	0.80 $\pm$ 0.01	0.88 $\pm$ 0.01	0.93 $\pm$ 0.03	1.09 $\pm$ 0.06	0.83 $\pm$ 0.02
Borehole 8D	32.11 $\pm$ 1.01	32.12 $\pm$ 1.01	48.75 $\pm$ 1.87	$>$100 $\pm$ 3.54	38.60 $\pm$ 1.20	41.35 $\pm$ 1.26
Styblinski-Tang 9D	$>$100 $\pm$ 3.05	$>$100 $\pm$ 0.00	$>$100 $\pm$ 5.33	$>$100 $\pm$ 3.03	$>$100 $\pm$ 6.31	$>$100 $\pm$ 4.12
PUMA560 9D	3.93 $\pm$ 0.15	3.23 $\pm$ 0.14	3.40 $\pm$ 0.14	3.93 $\pm$ 0.08	3.24 $\pm$ 0.14	3.4 $\pm$ 0.13
Adapted Welch 10D	99.51 $\pm$ 0.81	99.4 $\pm$ 0.75	$>$100 $\pm$ 0.55	99.79 $\pm$ 0.75	$>$100 $\pm$ 0.57	100.00 $\pm$ 0.67
Wing weight 10D	$>$100 $\pm$ 0.00	58.16 $\pm$ 4.37	63.1 $\pm$ 4.36	$>$100 $\pm$ 0.53	63.35 $\pm$ 4.15	$>$100 $\pm$ 1.69
Boston housing	11.28 $\pm$ 1.06	11.36 $\pm$ 1.04	11.42 $\pm$ 1.06	11.28 $\pm$ 1.06	11.56 $\pm$ 1.10	11.31 $\pm$ 1.00
Abalone	2.06 $\pm$ 0.10	2.09 $\pm$ 0.10	2.08 $\pm$ 0.10	2.05 $\pm$ 0.10	2.09 $\pm$ 0.10	2.08 $\pm$ 0.10
Naval propulsion	0.02 $\pm$ 0.00	0.02 $\pm$ 0.00	38.86 $\pm$ 0.60	62.61 $\pm$ 1.51	9.40 $\pm$ 0.16	3.45 $\pm$ 0.08
Forest fire	1.97 $\pm$ 0.05	1.87 $\pm$ 0.02	6.43 $\pm$ 0.69	10.43 $\pm$ 1.11	2.32 $\pm$ 0.14	3.32 $\pm$ 0.28
Parkinson’s	12.17 $\pm$ 0.02	12.40 $\pm$ 0.10	12.49 $\pm$ 0.16	11.97 $\pm$ 0.03	12.60 $\pm$ 0.16	12.78 $\pm$ 0.16

Table 3. Table 3: Rank of the methods corresponding to NLL (left) and RMSE (right). The best overall score is in bold (ties are possible in case of an overlap in confidence intervals).

	DE	HDE	AE	NTKGP-p.	RP-p.	RAFs
He et al. 1D	6,2	5,3	4,1	2,2	3,3	1,2
Forrester et al. 1D	4,2	5,2	3,2	6,2	2,2	1,1
Schaffer N.4 2D	5,1	2,4	6,3	3,1	4,3	1,2
Double pendulum 2D	3,1	3,3	2,2	1,1	1,2	1,1
Rastrigin 3D	2,2	1,1	3,3	2,2	1,1	1,1
Ishigami 3D	3,1	5,3	4,2	2,1	1,1	1,1
Environmental 4D	6,1	5,2	2,1	4,1	3,1	1,1
Griewank 4D	3,3	1,1	1,1	2,2	1,1	1,2
Roos & Arnold 5D	3,1	1,1	4,3	5,1	5,1	2,2
Friedman 5D	5,1	6,1	3,1	4,1	2,1	1,1
Planar arm torque 6D	5,1	4,1	3,1	1,1	2,2	2,1
Sum of powers 6D	5,1	6,1	3,1	4,1	2,1	1,1
Ackley 7D	4,5	1,2	2,4	3,5	2,3	1,1
Piston simulation 7D	1,1	3,2	2,5	2,6	2,3	2,4
Robot arm 8D	5,3	3,1	4,2	1,3	2,4	1,1
Borehole 8D	2,1	3,1	1,4	1,5	1,2	1,3
Styblinski-Tang 9D	3,2	5,5	2,1	4,3	1,1	1,4
PUMA560 9D	5,2	1,1	3,1	4,2	4,1	2,1
Adapted Welch 10D	6,1	2,1	4,3	5,1	3,2	1,1
Wing weight 10D	4,4	2,1	1,1	3,3	1,1	1,2
Boston housing	3,1	5,1	2,1	2,1	5,1	1,1
Abalone	5,1	6,1	3,1	4,1	1,1	2,1
Naval propulsion plant	1,1	4,1	3,4	2,5	2,3	2,2
Forest fire	3,2	1,1	1,5	1,6	2,3	1,4
Parkinson’s	1,2	6,3	4,3	5,1	3,3	2,3

Table 4. Table 4: Number of training data points and testing data points for each dataset.

	Training	Testing
He et al. 1D	20	50
Forrester et al. 1D	20	50
Schaffer N.4 2D	1000	2500
Double pendulum 2D	1000	2500
Rastrigin 3D	200	500
Ishigami 3D	2000	5000
Environmental 4D	200	500
Griewank 4D	200	500
Roos & Arnold 5D	200	500
Friedman 5D	200	500
Planar arm torque 6D	200	500
Sum of powers 6D	200	500
Ackley 7D	400	1000
Piston simulation 7D	200	500
Robot arm 8D	200	500
Borehole 8D	2000	5000
Styblinski-Tang 9D	400	1000
PUMA560 9D	3693	4499
Adapted Welch 10D	200	500
Wing weight 10D	2000	5000
Boston housing	354	152
Abalone	1880	2297
Naval propulsion plant	5370	6564
Forest fire	200	317
Parkinson’s	2643	3232

Equations

$${\mathbf{\sigma }}_p^2={\mathbf{\sigma }}_a^2+{\mathbf{\sigma }}_e^2.$$

$$\mathbf{y}=f(X)+ϵ,$$

$$\mathcal{L}({\mathbf{\theta }}_j)=\frac{1}{n}||\mathbf{y}-{\hat{\mathbf{y}}}_j|{|}_2^2+\frac{1}{n}||{\Gamma }^{1/2}({\mathbf{\theta }}_j-{\mathbf{\theta }}_{0,j})|{|}_2^2,$$

$$\mathbf{y}=L_1\sin ({\mathbf{\theta }}_1)+L_2\sin ({\mathbf{\theta }}_2)+ϵ$$

$$\mathbf{y}=\sqrt{4\pi }\cdot C(X)+ϵ,$$

$$C(X)=\frac{\mathbf{M}}{\sqrt{4\pi \mathbf{D}t}}\exp \left(\frac{-s^2}{4\mathbf{D}t}\right)+C^{\prime }(X)$$

$$C^{\prime }(X)=\frac{\mathbf{M}}{\sqrt{4\pi \mathbf{D}(t-\mathbf{\tau })}}\exp \left(\frac{-(s-\mathbf{L}{)}^2}{4\mathbf{D}(t-\mathbf{\tau })}\right)I(\mathbf{\tau }<t)$$

$${\mathbf{y}}_1=(M(\mathbf{q})\ddot{\mathbf{q}}+C(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}}{)}^T+ϵ,$$

$$M(\mathbf{q})=\left[\begin{array}{cc}0.2083+0.1250\cos ({\mathbf{q}}_2)& m({\mathbf{q}}_2)\\ m({\mathbf{q}}_2)& 0.0417\end{array}\right]$$

$$m({\mathbf{q}}_2)=0.0417+0.0625\cos ({\mathbf{q}}_2)$$

$$C(\mathbf{q},\dot{\mathbf{q}})=\left[\begin{array}{cc}-a\sin ({\mathbf{q}}_2){\dot{\mathbf{q}}}_2& a\sin ({\mathbf{q}}_2)({\dot{\mathbf{q}}}_1+{\dot{\mathbf{q}}}_2)\\ a\sin ({\mathbf{q}}_2){\dot{\mathbf{q}}}_1& 0\end{array}\right],$$

$$\mathbf{y}=2\pi \sqrt{\frac{\mathbf{M}}{\mathbf{k}+{\mathbf{S}}^2\frac{{\mathbf{P}}_0{\mathbf{V}}_0}{{\mathbf{T}}_0}\frac{{\mathbf{T}}_a}{{\mathbf{V}}^2}}}+ϵ,\text{where}$$

$$\mathbf{V}=\frac{\mathbf{S}}{2\mathbf{k}}\left(\sqrt{{\mathbf{A}}^2+4\mathbf{k}\frac{{\mathbf{P}}_0{\mathbf{V}}_0}{{\mathbf{T}}_0}{\mathbf{T}}_a}-\mathbf{A}\right)$$

$$\mathbf{A}={\mathbf{P}}_0\mathbf{qS}+19.62\mathbf{M}-\frac{{\mathbf{kV}}_0}{\mathbf{S}}.$$

$$\mathbf{y}=\sqrt{{\mathbf{u}}^2+{\mathbf{v}}^2}+ϵ,\text{where}$$

$$\mathbf{u}=\sum _{i=1}^4{\mathbf{L}}_i\cos \left(\sum _{j=1}^i{\mathbf{\theta }}_j\right)$$

$$\mathbf{v}=\sum _{i=1}^4{\mathbf{L}}_i\sin \left(\sum _{j=1}^i{\mathbf{\theta }}_j\right).$$

$$\mathbf{y}=\frac{2\pi {\mathbf{T}}_u({\mathbf{H}}_u-{\mathbf{H}}_l)}{\text{ln}(\mathbf{r}/{\mathbf{r}}_w)\left(1+\frac{2\mathbf{L}{\mathbf{T}}_u}{\text{ln}(\mathbf{r}/{\mathbf{r}}_w){\mathbf{r}}_w^2{\mathbf{K}}_w}+\frac{{\mathbf{T}}_u}{{\mathbf{T}}_l}\right)}+ϵ,$$

$${\mathbf{y}}_1=A(\mathbf{q}{)}^{-1}(\mathbf{\tau }-n(\mathbf{q},\dot{\mathbf{q}})-g(\mathbf{q}))+ϵ,$$

$$\mathbf{y}=0.036{\mathbf{S}}_w^{0.758}{\mathbf{W}}_{fw}^{0.0035}{\left(\frac{\mathbf{A}}{{\cos }^2(\mathbf{\Lambda })}\right)}^{0.6}{\mathbf{q}}^{0.006}{\mathbf{y}}^{\prime }$$

$${\mathbf{y}}^{\prime }={\mathbf{\lambda }}^{0.04}{\left(\frac{100{\mathbf{t}}_c}{\cos (\mathbf{\Lambda })}\right)}^{-0.3}({\mathbf{N}}_z{\mathbf{W}}_{dg}{)}^{0.49}+{\mathbf{y}}^{\prime \prime }$$

$${\mathbf{y}}^{\prime \prime }={\mathbf{S}}_w{\mathbf{W}}_p+ϵ,$$

$$\mathbf{y}=0.5+\frac{{\cos }^2(\sin (|{\mathbf{x}}_1^2-{\mathbf{x}}_2^2))-0.5}{(1+0.001({\mathbf{x}}_1^2+{\mathbf{x}}_2^2){)}^2}+ϵ,$$

$$\mathbf{y}=10d+\sum _{i=1}^d({\mathbf{x}}_i^2-10\cos (2\pi {\mathbf{x}}_i))+ϵ,$$

$$\mathbf{y}=\sum _{i=1}^d\frac{{\mathbf{x}}_i^2}{4000}-\prod _{i=1}^d\cos \left(\frac{{\mathbf{x}}_i}{\sqrt{i}}\right)+1+ϵ,$$

$$\mathbf{y}=-a\exp \left(-b\sqrt{\frac{1}{d}\sum _{i=1}^d{\mathbf{x}}_i^2}\right)-{\mathbf{y}}^{\prime }$$

$${\mathbf{y}}^{\prime }=\exp \left(\frac{1}{d}\sum _{i=1}^d\cos (c{\mathbf{x}}_i)\right)+a+\exp (1)+ϵ,$$

$$\mathbf{y}=\mathbf{x}\sin (\mathbf{x})+ϵ,$$

$$\mathbf{y}=(6\mathbf{x}-2{)}^2\sin (12\mathbf{x}-4)+ϵ,$$

$$\mathbf{y}=\sin ({\mathbf{x}}_1)+a{\sin }^2({\mathbf{x}}_2)+b{\mathbf{x}}_3^4\sin ({\mathbf{x}}_1)+ϵ,$$

Full text

Introduction

Recent advances in deep neural networks have demonstrated remarkable performance in a wide variety of applications, ranging from recommendation systems and improving user experience to natural language processing and speech recognition (Abiodun et al. 2018). Nevertheless, blindly relying on the outcome of these models can have harmful effects, especially in high-stake domains such as healthcare and autonomous driving, as models can provide inaccurate predictions when queried in out-of-distribution data points (Amodei et al. 2016). Consequently, correctly quantifying the uncertainty of models’ predictions is an admissible mechanism to distinguish where a model can or cannot be trusted, and thus, increases the transparency of models about their capabilities and limitations (Abdar et al. 2021). Uncertainty Quantification (UQ) is important for a variety of reasons. For instance, in order to preserve the model’s credibility, it is essential to report and communicate the encountered uncertainties regularly (Volodina and Challenor 2021). Additionally, models’ predictions are inevitably uncertain in most cases, which has to be addressed to increase their transparency, trustworthiness, and reliability.

In the machine learning literature, uncertainty is usually decomposed into two different types, namely aleatoric uncertainty and epistemic uncertainty (Kiureghian and Ditlevsen 2009). Aleatoric uncertainty, aka data uncertainty, refers to the inherent uncertainty that stems from the data itself, e.g., noise. On the other hand, epistemic uncertainty, also called model uncertainty, is the type of uncertainty that occurs due to the lack of sufficient data. While data uncertainty cannot be alleviated, model uncertainty can be addressed by e.g., acquiring more data. Let $\boldsymbol{\sigma}_{a}^{2}$ and $\boldsymbol{\sigma}_{e}^{2}$ denote the aleatoric and epistemic uncertainties, respectively. Since the distinction between the two is imprecise to some degree (Sullivan 2015), we focus on the predictive (total) uncertainty, which is defined as the sum of the two

[TABLE]

Accordingly, the approaches developed for uncertainty estimation can be categorized into three groups: Bayesian UQ methods, ensemble UQ methods, and a combination of both, i.e., Bayesian ensemble UQ (Abdar et al. 2021). In this paper, we focus on ensemble UQ techniques, either Bayesian or non-Bayesian, as this group is less explored compared to the solely Bayesian techniques. An ensemble model aggregates the predictions of multiple individual base-learners (or ensemble members), which in our case are neural networks (NNs), and the empirical variance of their predictions gives an approximate measure of uncertainty. The idea behind this heuristic is highly intuitive: the more the base-learners disagree on the outcome, the more uncertain they are. Therefore, the goal of ensemble members is to have a great level of disagreement (variability) in the areas where little or no data is available, and to have a high level of agreement in regions with abundance of data (Pearce et al. 2018).

In this paper, we propose a novel method, called Random Activation Functions Ensemble (RAFs Ensemble), for a more robust uncertainty estimation in (deep) neural networks. RAFs Ensemble is developed on top of Anchored Ensemble technique, proposed by (Pearce et al. 2018), however, instead of initializing each NN member in the ensemble with the same activation function, the NNs in RAFs Ensemble are accommodated with different (random) activation functions in the hidden layers. This simple, yet crucial, modification greatly improves the overall diversity of the ensemble, which is one of the most important components in forming a successful ensemble. We empirically show that RAFs Ensemble provides high quality uncertainty estimates compared to five state-of-the-art ensemble methods, that is Deep Ensemble (Lakshminarayanan, Pritzel, and Blundell 2017), Neural Tangent Kernel Gaussian Process Parameter Ensemble (He, Lakshminarayanan, and Teh 2020), Anchored Ensemble (Pearce et al. 2018), Bootstrapped Ensemble of NNs Coupled with Random Priors (Osband, Aslanides, and Cassirer 2018), and Hyperdeep Ensemble (Wenzel et al. 2020). The comparisons are performed in a wide range of regression tasks on both synthetic and real-world datasets in terms of negative log-likelihood and root mean squared error.

Related Work

Uncertainty Quantification (UQ) is an active field of research and various methods have been proposed to efficiently estimate the uncertainty of machine learning models (see Abdar et al. 2021 for an extensive overview). While most research focuses on Bayesian deep learning (Srivastava et al. 2014; Blundell et al. 2015; Sensoy, Kandemir, and Kaplan 2018; Fan et al. 2020; Järvenpää, Vehtari, and Marttinen 2020; Charpentier, Zügner, and Günnemann 2020), deep ensemble methods, which benefit from the advantages of both deep learning and ensemble learning, have been recently leveraged for empirical uncertainty quantification (Egele et al. 2021; Hoffmann, Fortmeier, and Elster 2021; Brown, Bhuiyan, and Talbert 2020; Althoff, Rodrigues, and Bazame 2021). Although Bayesian UQ methods have solid theoretical foundation, they often require significant changes to the training procedure and are computationally expensive compared to non-Bayesian techniques such as ensembles (Egele et al. 2021; Rahaman and Thiery 2021; Lakshminarayanan, Pritzel, and Blundell 2017).

Lakshminarayanan, Pritzel, and Blundell (2017) are among the first to challenge Bayesian UQ methods by proposing Deep Ensemble, a simple and scalable technique, that demonstrates superb empirical performance on a variety of datasets. However, one of the challenges of ensemble techniques when quantifying uncertainty is that they tend to give overconfident predictions (Amodei et al. 2016). To address this, Pearce et al. (2018) propose to also regularize the model’s parameters w.r.t. the initialization values, instead of zero, leading to Anchored Ensembles, which additionally allows for performing Bayesian inference in NNs. He, Lakshminarayanan, and Teh (2020) relate Deep Ensembles to Bayesian inference using neural tangent kernels. Their method, i.e., Neural Tangent Kernel Gaussian Process Parameter Ensemble (NTKGP-param), trains all layers of a finite width NN, obtaining an exact posterior interpretation in the infinite width limit with neural tangent kernel parameterization and squared error loss. They prove that NTKGP-param is always more conservative than Deep Ensemble, yet, its advantages are generally not clear in practice.

A prominent advance to the Bayesian ensemble UQ methods is the bootstrapped ensemble of NNs coupled with random priors, proposed by (Osband, Aslanides, and Cassirer 2018), in which, the random prior function and neural models share an input and a summed output, but the networks are the only trainable parts, while the random prior remains untrained throughout the whole process. Furthermore, Wenzel et al. (2020) exploit an additional source of randomness in ensembles by designing ensembles not only over weights, but also over hyperparameters. Their method, called Hyperdeep Ensemble, demonstrates high accuracy for a number of different classification tasks. Nevertheless, despite the recent contributions in ensemble UQ methods, the research in this direction still needs further advancement.

Toward Robust Uncertainty Estimation

Preliminaries

Following the notations of (Lakshminarayanan, Pritzel, and Blundell 2017), let $S_{train}$ be a training dataset consisting of $n$ independently and identically drawn (i.i.d.) data points, $S_{train}=\{\boldsymbol{x}_{i},y_{i}\}_{i=1}^{n}$, where $\boldsymbol{x}_{i}\in\mathbb{R}^{d}$ denotes a $d$-dimensional feature vector and $y_{i}\in\mathbb{R}$ is a scalar output. Similarly, $S_{test}$ indicates the test set. Subsequently, $X$ represents the design matrix and $\boldsymbol{y}$ indicates the output vector, where $(S_{train}.X,S_{train}.\boldsymbol{y})$ and $(S_{test}.X,S_{test}.\boldsymbol{y})$ represent the train and test sets, respectively. Without the loss of generality, we consider the regression tasks of the form

[TABLE]

where $\epsilon$ is a normally distributed constant noise, i.e., $\epsilon\sim\mathcal{N}(0,\boldsymbol{\sigma}_{a}^{2})$, and is assumed to be known. The goal is hence to quantify the predictive uncertainty $\boldsymbol{\sigma}_{p}^{2}$ associated with $S_{test}.\boldsymbol{y}$, while optimizing $f$ on the training data.

We adapt the regularized loss function from the Anchored Ensemble technique (Pearce et al. 2018), in which, the regularization of the models’ parameters are carried out w.r.t. their initialization values instead of zero. Consequentially, given $\boldsymbol{\theta}_{j}$ as the parameters of the $j$th base-learner, the objective function is as follows

[TABLE]

where $\boldsymbol{\theta}_{0,j}$ is derived from the prior distribution, $\boldsymbol{\theta}_{0,j}\sim\mathcal{N}(\boldsymbol{\mu}_{0},\Sigma_{0})$, and $\Gamma$ is the regularization matrix. Furthermore, minimizing this objective allows for performing Bayesian inference in NNs. However, this technique only models the epistemic uncertainty, while aleatoric uncertainty is assumed to be constant (Pearce et al. 2018), which is a limitation, as it is not always possible to distinguish the different origins or types of uncertainty in practice (see Equation 1).

Therefore, in this paper, we aim at enhancing the performance of the ensemble toward a more robust uncertainty estimation. The literature suggests that diversifying the ensembles is effective in improving their predictive performance both theoretically and empirically (Zhou 2012; Zhang and Ma 2012; Hansen and Salamon 1990; Krogh and Vedelsby 1994). Ideally, diversity is achieved when the predictions made by each model in the ensemble are independent and uncorrelated. However, generating diverse ensemble members is not a straightforward task. The main impediment is the fact that each neural network is trained on the same training data to solve the same problem, which usually results in a high correlation among the individual base-learners (Zhou 2012). In the subsequent section, we introduce a simple technique to efficiently improve the overall diversity of the ensemble for a more reliable uncertainty quantification.

RAFs Ensemble

In this section, we present Random Activation Functions (RAFs) Ensemble for uncertainty estimation, which can be extended to all ensemble methods in terms of methodological modification. When a (Bayesian) ensemble is leveraged to estimate the uncertainty of a deep neural network model, we propose to increase the diversity of predictions among the ensemble members using varied activation functions (AFs), in addition to the random initialization of the parameters. To do so, instead of initializing the neural networks with the same activation function, each NN is accommodated with a different (random) activation function. Subsequently, distinct activation functions account for different non-linear properties introduced to each ensemble member, therewith improving the overall diversity of the ensemble.

As mentioned previously, the ensemble diversity is one of the most important building blocks when it comes to creating a successful ensemble (Hansen and Salamon 1990). Hence, it might be preferable to combine the predictions of top-performing base-learners with the predictions of weaker ones (Zhou 2012). Otherwise, stacking only strong models will likely result in a poor ensemble as the predictions made by the models will be highly correlated, and thus, the ensemble diversity will be greatly limited. Therefore, the choice of activation functions should be motivated purely by their variability and not their appropriateness for the task at-hand.

Let $\boldsymbol{\mu}_{0}$ be the prior means, $\Sigma_{0}$ be the prior covariance, $\boldsymbol{\hat{\sigma}}_{a}^{2}$ be an estimate of data noise, $m$ denote the number of base-learners, and $NN_{j}$ indicate the $j_{\text{th}}$ member, the entire procedure for both training and prediction is summarized in Algorithm 1. In this algorithm, a regularization matrix is first created and a set of activation functions is defined (line 1-2). Then, the NNs in the ensemble are trained to minimize the loss function defined in Equation 2 with stochastic gradient descent, using arbitrary optimizer and no early stopping (line 3-13). Note that if the size of the ensemble $m$ is smaller or equal to the cardinality of the AFs set $k$, then each NN is trained with a different activation function, and with random functions from the set, otherwise (line 7-11). Consequently, predictions are made with each ensemble member (line 14-16), which are then averaged and an estimate of the predictive uncertainty is computed (line 17-19).

Empirical Study

Experimental Setups

In the experiments, the base-learners of RAFs Ensemble are multilayer perceptrons that consist of one hidden layer of 100 neurons. The ensemble size $m$ is set to five. This is standard for the implementations of all methods in this paper, as $m=5$ proved to be empirically sufficient for obtaining predictive uncertainty estimates in the experiments. In addition, we choose a set of seven activation functions which is comprised of (i) Gaussian Error Linear Unit (GELU) (Hendrycks and Gimpel 2016), (ii) Softsign (Turian, Bergstra, and Bengio 2009), (iii) Swish (Ramachandran, Zoph, and Le 2018), (iv) Scaled Exponential Linear Unit (SELU) (Klambauer et al. 2017), (v) hyperbolic tangent (tanh), (vi) error activation function, and (vii) linear (identity) activation function. Furthermore, the number of testing samples is set to be always larger than the number of training points $n$ to detail the uncertainty. Moreover, to account for epistemic uncertainty, the synthetic testing feature vectors $\boldsymbol{x}\in S_{test}$ range over wider intervals compared to $\boldsymbol{x}\in S_{train}$ and both are sampled uniformly at random. The code is available at https://github.com/YanasGH/RAFs˙code for reproducibility.

Baselines.

We include five state-of-the-art methods as baselines for empirical comparison with RAFs Ensemble as follows. (i) DE (Lakshminarayanan, Pritzel, and Blundell 2017), (ii) AE (Pearce et al. 2018), (iii) HDE (Wenzel et al. 2020), (iv) RP-param (Osband, Aslanides, and Cassirer 2018), and (v) NTKGP-param (He, Lakshminarayanan, and Teh 2020), on both synthetic and real-world datasets with different dimensionalities (see the Technical appendix for a detailed overview). To ensure fair comparison between the UQ techniques, roughly the same amount of time has been put into hyperparameter tuning for each method.

Synthetic Data.

We generate multiple synthetic datasets that fall into four categories: physical models (PM), many local minima (MLM), trigonometric (T), and others (O). Each set in the PM category is generated from a physical mathematical model, such that all values in $S_{train}$ and $S_{test}$ are achievable in the real world. Generally, the PM datasets have complex modeling dynamics and can be characterized as having predominant epistemic uncertainty due to the considerably wider testing sampling regions by design. Similarly, the MLM data, generated from functions with many local minima, are also designed so that the model uncertainty is higher than the aleatoric one. These datasets are usually hard to approximate due to their inherent high-nonlinearity and multimodality. Another category with higher epistemic uncertainty is trigonometric, such as data generated by (He, Lakshminarayanan, and Teh 2020) and (Forrester, Sobester, and Keane 2008), where the training data is partitioned into two equal-sized clusters in order to detail uncertainty on out-of-distribution data (see Figure 1). In contrast, the predominant type of uncertainty in the O category is aleatoric. This category includes datasets generated from various functions such as rational and product integrand functions. It is distinguished from the rest of the categories by its high interaction effects. The dimensionality of all datasets can range from one to ten and we consider two datasets per dimension, thus, the total number of synthetic data is 20. More detail on how the data is created can be found in the Technical appendix.

Real-world Data.

Additionally, we use five real-world datasets for evaluation: Boston housing, Abalone shells (Nash et al. 1994), Naval propulsion plant (Coraddu et al. 2014), Forest fire (Cortez and de Jesus Raimundo Morais 2007), and Parkinson’s disease dataset (Little et al. 2007). To account for aleatoric uncertainty (some) context factors are disregarded, such that this type of uncertainty is characteristically high (see Technical appendix for more details).

Evaluation Criteria.

We employ two evaluation criteria to gauge the overall performance of the trained models, namely calibration and robustness to the distribution shift. Both measures are inspired by the practical applications of NNs, as generally there is no theoretical evidence for evaluating uncertainty estimates (Abdar et al. 2021). Calibration is defined as the analytical process of adjusting the inputs with the purpose of making the model to predict the actual observations as precisely as possible (Bijak and Hilton 2021). The quality of calibration can be measured by proper scoring rules such as negative log-likelihood (NLL). NLL is a common choice when it comes to evaluating UQ estimates, as it depends on predictive uncertainty (Lakshminarayanan, Pritzel, and Blundell 2017). Additionally, due to its practical applicability in a wide spectrum of regression tasks, root mean squared error (RMSE) is measured, although it does not depend on the estimated uncertainty (Lakshminarayanan, Pritzel, and Blundell 2017), but serves as a proxy and a secondary assessor of the performance. Moreover, to measure the robustness/generalization of methods to distributional shift, we test the models in out-of-distribution settings, such as the synthetic datasets by (Forrester, Sobester, and Keane 2008; He, Lakshminarayanan, and Teh 2020).

Performance Results

Qualitative Comparison.

Figure 1 shows the performance of different methods compared to a Gaussian process with a neural tangent kernel (NTKGP analytic) as a reference, on a 1D toy dataset generated from $\boldsymbol{y}=\boldsymbol{x}\text{sin}(\boldsymbol{x})+\boldsymbol{\epsilon}$ (dashed line). The plots demonstrate that DE, HDE, and AE provide narrow uncertainty bounds in areas where no data has been observed by the model, which translates to high confidence in OOD data. On the contrary, NTKGP-param, RP-param, and RAFs Ensemble bound their uncertainty estimates with wider intervals in areas with no data, accounting for adequate quantification of epistemic uncertainty, while also indicating robustness to OOD data. Among these methods, RAFs Ensemble provides the widest uncertainty which is reasonable considering the amount of data that is available to the methods over each area. Moreover, this observation is quantitatively validated as RAFs Ensemble achieves the lowest NLL compared to the other methods (see Table 1).

Overall Performance.

We evaluate the overall performance of all methods in terms of both NLL and RMSE. The outcomes of comparing RAFs Ensemble with five baseline methods on twenty synthetic and five real-world datasets are outlined in Table 1 and Table 2. The results illustrate that our approach outperforms the competitors in most scenarios. Furthermore, Table 3 summarizes the obtained results in terms of ranking, in which the methods are ranked based on their performance for a particular dataset. The left integer corresponds to NLL, while the right one points to RMSE, and the bold values indicate the best-performing method.

Discussion.

The obtained results in this section illustrate that DE has good uncertainty estimates with respect to NLL for the real-world datasets, and it takes the first place for Naval propulsion and Parkinson’s datasets. For the rest of the data categories, when compared to the other methods, DE fails to provide strong performance, usually scoring a very low NLL rank. Therefore, this indicates that Deep Ensemble might have difficulty quantifying epistemic uncertainty in general as displayed by the experiments in this paper, but seemingly manages to capture aleatoric uncertainty well.

Unlike DE, the HDE provides outwardly reliable estimates for datasets with many local minima, despite its unimpressive overall results when compared to the other methods. However, both DE and HDE can produce uncertainty bounds that are unreasonably narrow in areas with unobserved data, as shown in Figure 1 and noted by (Heiss et al. 2021).

Nonetheless, AE demonstrates good performance in the dataset categories that exhibit higher epistemic uncertainty such as the physical models. This is due to the fact that AE is designed for capturing model uncertainty, while aleatoric uncertainty is assumed to be constant. Accordingly, AE achieves inferior performance on the real-world datasets, as those generally have more data uncertainty appropriated.

On the other hand, NTKGP-param achieves its finest performance for datasets in the physical model category, which is normally associated with substantial model uncertainty. A credible rationale to explain this insight is the fact that NTKGP-param tends to be more conservative than Deep Ensemble. However, it is generally unclear in which situations this is beneficial since the ensemble members of NTKGP-param will always be misspecified in practice according to (He, Lakshminarayanan, and Teh 2020).

Furthermore, RP-param manages to rank comparatively high for real-world datasets as well as trigonometric data, that contain vast amounts of aleatoric and epistemic uncertainty, respectively, indicating that it does not quantify either type of uncertainty better than the other. This observation serves as a demonstration that RP-param generalizes well for different types of datasets that exhibit broad characteristics. However, this technique fails to deliver low NLL scores on some occasions, which might be attributed to the fact that RP-param is based on bootstrapping. While bootstrapping can be a successful strategy for inducing diversity, it can sometimes harm the performance when the base-learners have multiple local optima, as is a common case with NNs (Lakshminarayanan, Pritzel, and Blundell 2017).

Nevertheless, RAFs Ensemble outperforms RP-param, and every other method in the comparisons, for 13 out of 25 datasets. In terms of NLL, our approach does not rank below the second place for any data, which is consistent with the strong results from Table 1. Meanwhile, the RMSE scores of this method are altogether satisfactory, although not as prominent compared to the NLL scores. In agreement with the overall outstanding results, RAFs Ensemble holds the highest NLL rank for all data from MLM and T categories, which can be contemplated as a concluding statement regarding its capabilities to estimate epistemic uncertainty and challenging multimodality. Among all categories, the real-world datasets are least favored by RAFs Ensemble, primarily due to their high level of aleatoric uncertainty. This indicates that RAFs Ensemble captures model uncertainty slightly better than aleatoric uncertainty. Nonetheless, the empirical superiority of this technique is due to the exhaustively exploited added source of randomness via random activation functions, combined with method simplicity and Bayesian behavior, resulted from the anchored loss (Equation 2). This successful combination leads to greatly improved diversity among ensemble members, which can be further confirmed by a direct comparison between RAFs Ensemble and AE. Note that even though RAFs Ensemble does not provide as prominent results with respect to RMSE in the higher dimensional datasets as it does in datasets of lower dimensions, it still achieves better or on par results compared to the state-of-the-art methods. In addition, RAFs Ensemble can be deployed in both complex and straightforward settings. On a related note, while DE struggles when dealing with high multimodality and RP-param underperforms when the dataset has interaction effects (from “others” category), RAFs excels in both such settings.

Scalability to higher dimensions and larger networks.

To test the scalability of RAFs Ensemble, we compare it with the strongest baseline, RP-param, on two additional real-world datasets, i.e., a 65-dimensional data with around 20k samples and a 40-dimensional data with almost 40k samples. The former is the superconductivity dataset, where the goal is to predict the critical temperature of superconductors (Hamidieh 2018). The latter summarizes features about articles, where the target is the number of shares in social networks (Fernandes, Vinagre, and Cortez 2015). Both methods utilize the same neural architecture for their base-learners, that is two hidden layers of 128 hidden neurons each, which is more complex than the previous experiments. The conclusion of these experiments is conclusive in favor of our approach. RAFs Ensemble scores NLL of 5.49 and 25.89 on the first and second dataset, respectively, while RP-param scores NLL of over 100 on both datasets.

Confidence vs. Error.

We further analyze the relation between the RMSE and the precision thresholds in order to examine the confidence of each method in the prediction task. Figure 2 displays the confidence versus error plots for one synthetic and one real-world dataset, i.e., Friedman and Abalone (see the Technical Appendix for more detail). In this figure, for each precision threshold $\tau$, the RMSE is plotted for examples where the predicted precision $\boldsymbol{\sigma}_{p}^{-2}$ is larger than the threshold $\tau$, demonstrating confidence. In general, reliable estimates are expected to have decreasing error when the confidence is increasing. For Friedman dataset, it is clear that RAFs Ensemble delivers well-calibrated estimates, which is especially in contrast with DE, NTKGP-param, and HDE (Figure 2(a)). However, for the Abalone data, RP-param demonstrates the most reliable behavior, although RAFs Ensemble meets its performance at the last precision threshold (Figure 2(b)). Overall, our approach sustains lower error over most precision thresholds compared to the majority of the other methods, and this contrast in performance is emphasized as the predictions get more confident.

Ablation.

We study the effect of number of base-learners in the ensemble on the quality of UQ, which also measures the sensitivity of the results to the cardinality of the set of AFs $k$. We conduct an experiment on two different datasets, one synthetic (PUMA590) and one real-world (Abalone), where the results in terms of NLL are represented in Figure 3. Note that Figure 3(b) is shown in log-scale for better visibility. According to the theory, in the limit of infinite number of ensemble members, the ensemble error converges to zero (Hansen and Salamon 1990). However, practically speaking, five NNs in the ensemble provide optimal results regarding the trade-off between empirical performance and computational time (Lakshminarayanan, Pritzel, and Blundell 2017), which is also the case in our experiments. This is further confirmed by the plot in Figure 3(a). In addition, for the PUMA590 dataset, it seems that RAFs Ensemble’s performance is not impacted negatively by the number of NNs in the ensemble. Moreover, an interesting observation is the steep through for seven NNs (equal to $k$) in Figure 3(b), which is an indication that there might be a correlation between $k$ and the performance in some cases. A plausible reason for this is the fact that the additional source of randomness is utmostly exploited via a different activation function.

To further confirm the effectiveness of the random activation functions, we evaluate the performance of RAFs Ensemble (of five NNs) in terms of NLL w.r.t. different cardinalities $k$ of the set of AFs. The dataset used for this experiment is the superconductivity data. As the results in Figure 4 clearly suggest, by increasing the cardinality $k$, NLL has a decreasing pattern, which shows that having more random AFs significantly improves the performance of the ensemble.

Moreover, we combine our approach with RP-param instead of AE to show that RAFs can be methodologically applied to any ensemble technique. We evaluate the performance of this combination on the Parkinson’s dataset, using the same network architecture for fair comparison. The obtained results demonstrate that applying RAFs to RP-param leads to reducing the original NLL score of $>100$ to 48.66, which is in line with the results we get when comparing AE with RAFs Ensemble and is a further proof that the methodology indeed increases the performance.

Conclusions

We introduced a novel method, Random Activation Functions Ensemble, for a more robust uncertainty estimation in approaches based on neural networks, in which, each network in the ensemble is accomodated with a different (random) activation function to increase the diversity of the ensemble. The empirical study illustrates that our approach achieves excellent results in quantifying both epistemic and aleatoric uncertainty compared to five state-of-the-art ensemble uncertainty quantification methods on a series of regression tasks across 25 datasets, which proved there does not have to be a trade-off between simplicity and strong empirical performance. Furthermore, the properties of datasets such as dimensionality or complexity of modeling dynamics do not appear to affect RAFs Ensemble negatively, which also demonstrates robustness in out-of-distribution settings.

Appendix A Appendix

Appendix B Synthetic datasets

The number of training data points and testing data points for each dataset is shown in Table 4.

Physical models

Both the training data and testing data of every dataset in this dataset category are sampled from realistic ranges, so that all values are also possible in a real world setting.

Double pendulum 2D

A double pendulum, is a dynamical system, which consists of a pendulum with another pendulum attached to its end (Levien and Tan 1993). For the purposes of this work, the length of both pendulum ropes $L_{1}$ and $L_{2}$ is kept constant $L_{1}=L_{2}=1$ and the response variable $y_{i}$ that is being modeled is the horizontal position of the lower pendulum mass given $\boldsymbol{\theta}_{1}$ and $\boldsymbol{\theta}_{2}$:

[TABLE]

For generating the training dataset, the inputs $\boldsymbol{\theta}_{1}$ and $\boldsymbol{\theta}_{2}$ range over $[\frac{-2\pi}{3},\frac{\pi}{6}]$, while for testing both inputs are sampled from $[-\pi,\pi]$. Additionally, $\epsilon\sim\mathcal{N}(0,0.1^{2})$.

Environmental Model 4D

The Environmental Model function is a pollutant diffusion problem that models a pollutant spill at two locations caused by a chemical accident (Bliznyuk et al. 2008):

[TABLE]

where $\epsilon\sim\mathcal{N}(0,0.1^{2})$, the response variable $\boldsymbol{y}$ is the scaled the concentration of the pollutant $C(X)$ at the space-time vector $\boldsymbol{(s,t)}$:

[TABLE]

[TABLE]

where $\boldsymbol{M}$ is the mass of the pollutant spill, $\boldsymbol{D}$ is the diffusion rate in the channel, $\boldsymbol{L}$ is the location of the second spill, $\boldsymbol{\tau}$ is the time of the second spill and $I$ is the indicator function. For generating this dataset, $s$ and $t$ are fixed: $s=1$ and $t=40.1$. The ranges of the input values for training are as follows: $\boldsymbol{M}\in[7,13],\boldsymbol{D}\in[0.02,0.12],\boldsymbol{L}\in[0.01,3],\boldsymbol{\tau}\in[30.01,30.295]$. For testing, those ranges are wider: $\boldsymbol{M}\in[5,15],\boldsymbol{D}\in[0,0.15],\boldsymbol{L}\in[0.01,3.2],\boldsymbol{\tau}\in[23.71,31]$.

Planar arm torque 6D

Planar arm torque dataset approximates the first motor’s torque in the inverse dynamics of a Planar 2D Arm (Cully et al. 2018):

[TABLE]

where $\epsilon\sim\mathcal{N}(0,0.1^{2})$, $\boldsymbol{q}$ is a 2-dimensional vector denoting the articular position, $\dot{\boldsymbol{q}}$ is a 2-dimensional vector denoting the articular velocity, $\ddot{\boldsymbol{q}}$ is a 2-dimensional vector denoting the articular acceleration. $M(\boldsymbol{q})$ is the mass matrix and $C(\boldsymbol{q},\dot{\boldsymbol{q}})$ is the matrix of Coriolis and centrifugal forces:

[TABLE]

[TABLE]

[TABLE]

where $a=0.0625$. The features span though the following intervals for training: $\boldsymbol{q}_{1},\boldsymbol{q}_{2}\in[-\frac{\pi}{2},\frac{\pi}{2}]$, $\dot{\boldsymbol{q}}_{1},\dot{\boldsymbol{q}}_{2}\in[-\pi,\pi]$, $\ddot{\boldsymbol{q}}_{1},\ddot{\boldsymbol{q}}_{2}\in[-\pi,\pi]$. For generating the testing feature values $\boldsymbol{q}_{1},\boldsymbol{q}_{2}\in[-\pi,\pi]$, $\dot{\boldsymbol{q}}_{1},\dot{\boldsymbol{q}}_{2}\in[-2\pi,2\pi]$, $\ddot{\boldsymbol{q}}_{1},\ddot{\boldsymbol{q}}_{2}\in[-2\pi,2\pi]$ are used.

Piston simulation 7D

The Piston Simulation function models the circular motion of a piston within a cylinder. The piston’s linear motion is transformed into circular motion by connecting a linear rod to a disk. Thus, the faster the piston moves inside the cylinder, the quicker the disk rotation and thus, the faster the engine runs. The response variable $\boldsymbol{y}$ is the cycle time in seconds (Ben-Ari and Steinberg 2007; Freitas 1999), which is affected by the features via a chain of nonlinear functions:

[TABLE]

[TABLE]

[TABLE]

In the above equations $\boldsymbol{M}$ is the piston weight (kg), $\boldsymbol{S}$ is the piston surface area ($\text{m}^{2}$), $\boldsymbol{V}_{0}$ is the initial gas volume ($\text{m}^{3}$), $\boldsymbol{k}$ is the spring coefficient (N/m), $\boldsymbol{P}_{0}$ is the atmospheric pressure (N/$\text{m}^{2}$), $\boldsymbol{T}_{a}$ is the ambient temperature (K), $\boldsymbol{T}_{0}$ is the filling gas temperature (K) and the error term $\epsilon\sim\mathcal{N}(0,0.1^{2})$. The training features are from the following intervals: $\boldsymbol{M}\in[30,60]$, $\boldsymbol{S}\in[0.005,0.020]$ , $\boldsymbol{V}_{0}\in[0.002,0.010]$, $\boldsymbol{k}\in[1000,5000]$, $\boldsymbol{P}_{0}\in[90000,110000]$, $\boldsymbol{T}_{a}\in[290,296]$ , $\boldsymbol{T}_{0}\in[340,360]$. Comparably, the testing input values are: $\boldsymbol{M}\in[0,90]$, $\boldsymbol{S}\in[0.005,0.03]$ , $\boldsymbol{V}_{0}\in[0,0.015]$, $\boldsymbol{k}\in[10,6000]$, $\boldsymbol{P}_{0}\in[80000,120000]$, $\boldsymbol{T}_{a}\in[285,300]$ , $\boldsymbol{T}_{0}\in[300,400]$.

Robot arm 8D

The Robot Arm function models the position of a four-segment robot arm and the response is the distance from the end of the robot arm to the origin (An and Owen 2001):

[TABLE]

[TABLE]

[TABLE]

The shoulder of the robot arm is fixed at the origin, however, each of the four segments is positioned at angle $\boldsymbol{\theta}_{j}$ and has length $\boldsymbol{L}_{i}$. Each input variable for the training set is generated from $\boldsymbol{\theta}_{j}\in[0,\pi]$ and $\boldsymbol{L}_{i}\in[0,0.5]$, while the test input variables $\boldsymbol{\theta}_{j}$ and $\boldsymbol{L}_{i}$ range over $[0,2\pi]$ and $[0,1]$ respectively. Finally, $\epsilon\sim\mathcal{N}(0,0.1^{2})$.

Borehole 8D

The Borehole function models water flow through a borehole and thus, the response variable is the water flow rate ($\text{m}^{3}$/yr):

[TABLE]

where $\epsilon\sim\mathcal{N}(0,0.1^{2})$, $\boldsymbol{r}_{w}$ is the radius of a borehole (m), $\boldsymbol{r}$ is the radius of influence (m), $\boldsymbol{T}_{u}$ and $\boldsymbol{T}_{l}$ are the transmissivities of respectively upper and lower aquifers ($\text{m}^{2}$/yr), $\boldsymbol{H}_{u}$ and $\boldsymbol{H}_{l}$ are the potentiometric heads of respectively upper and lower aquifers (m), $\boldsymbol{L}$ is the length of a borehole (m) and $\boldsymbol{K}_{w}$ is the hydraulic conductivity of borehole (m/yr) (An and Owen 2001). The features for training are sampled from $\boldsymbol{r}_{w}\in[0.05,0.15]$, $\boldsymbol{r}\in[100,50000]$, $\boldsymbol{T}_{u}\in[63070,115600]$, $\boldsymbol{T}_{l}\in[63.1,116]$, $\boldsymbol{H}_{u}\in[990,1110]$, $\boldsymbol{H}_{l}\in[700,820]$, $\boldsymbol{L}\in[1120,1680]$, $\boldsymbol{K}_{w}\in[9855,12045]$, while the testing input - $\boldsymbol{r}_{w}\in[0.01,0.2]$, $\boldsymbol{r}\in[90,50010]$, $\boldsymbol{T}_{u}\in[63020,115650]$, $\boldsymbol{T}_{l}\in[60,120]$, $\boldsymbol{H}_{u}\in[950,1150]$, $\boldsymbol{H}_{l}\in[650,900]$, $\boldsymbol{L}\in[1100,1700]$, $\boldsymbol{K}_{w}\in[9800,12100]$.

PUMA560 9D

PUMA560 dataset is generated from a realistic simulation of the dynamics of a Puma 560 robot arm (Ghahramani 1996). The task is to predict the articular acceleration of one of the links of the robot arm:

[TABLE]

where $\boldsymbol{q}$ and $\dot{\boldsymbol{q}}$ are 3-dimensional vectors denoting respectively the angular positions and angular velocities of each of the three links, $\boldsymbol{\tau}$ is a 3-dimensional vector representing the torques at the three joints, $n(\boldsymbol{q},\dot{\boldsymbol{q}})$ is the Coriolis and centrifugal effects, $A(\boldsymbol{q})$ is the inertia matrix, $g$ is the gravity and $\epsilon\sim\mathcal{N}(0,0.4^{2})$ is the Gaussian noise. The test and train features are sampled from $\boldsymbol{q}_{1},\boldsymbol{q}_{2},\boldsymbol{q}_{3}\in\beta[\frac{-\pi}{2},\frac{\pi}{2}]$, $\dot{\boldsymbol{q}}_{1},\dot{\boldsymbol{q}}_{2},\dot{\boldsymbol{q}}_{3}\in\beta[\frac{-\pi}{2},\frac{\pi}{2}]$ and $\boldsymbol{\tau}_{1},\boldsymbol{\tau}_{2},\boldsymbol{\tau}_{3}\in\beta[\frac{-1}{2},\frac{1}{2}]$ with fixed $\beta=1.2$ to control the nonlinearity of the dataset. Therefore, this dataset can be considered as highly nonlinear and noisy.

Wing weight 10D

The Wing Weight function models a light aircraft wing, where the response is the wing’s weight (Forrester, Sobester, and Keane 2008):

[TABLE]

[TABLE]

[TABLE]

where $\boldsymbol{S}_{w}$ is the wing area ($\text{ft}^{2}$), $\boldsymbol{W}_{fw}$ is the weight of fuel in the wing (lb), $\boldsymbol{A}$ is the aspect ratio, $\boldsymbol{\Lambda}$ is the quarter-chord sweep (degrees), $\boldsymbol{q}$ is the dynamic pressure at cruise (lb/$\text{ft}^{2}$), $\boldsymbol{\lambda}$ is the taper ratio, $\boldsymbol{t}_{c}$ is the aerofoil thickness to chord ratio, $\boldsymbol{N}_{z}$ is the ultimate load factor, $\boldsymbol{W}_{dg}$ is the flight design gross weight (lb), $\boldsymbol{W}_{p}$ is the paint weight (lb/$\text{ft}^{2}$) and $\epsilon\sim\mathcal{N}(0,0.1^{2})$. The ranges of the features for training for each value are $\boldsymbol{S}_{w}\in[150,200]$, $\boldsymbol{W}_{fw}\in[220,300]$, $\boldsymbol{A}\in[6,10]$, $\boldsymbol{\Lambda}\in[-10,10]$, $\boldsymbol{q}\in[16,45]$, $\boldsymbol{\lambda}\in[0.5,1]$, $\boldsymbol{t}_{c}\in[0.08,0.18]$, $\boldsymbol{N}_{z}\in[2.5,6]$, $\boldsymbol{W}_{dg}\in[1700,2500]$ and $\boldsymbol{W}_{p}\in[0.025,0.08]$, whereas for testing the intervals contain values outside the usual ranges - $\boldsymbol{S}_{w}\in[100,250]$, $\boldsymbol{W}_{fw}\in[200,320]$, $\boldsymbol{A}\in[0,15]$, $\boldsymbol{\Lambda}\in[-20,20]$, $\boldsymbol{q}\in[0,60]$, $\boldsymbol{\lambda}\in[0,1.5]$, $\boldsymbol{t}_{c}\in[0.05,0.25]$, $\boldsymbol{N}_{z}\in[0.5,8]$, $\boldsymbol{W}_{dg}\in[1000,3000]$ and $\boldsymbol{W}_{p}\in[0,0.1]$.

Many local minima functions

The many local minima functions are extremely hard to be approximated due to their high nonlinearity and complexity.

Schaffer N.4 2D

Schaffer N.4 function, proposed in (Mishra 2006):

[TABLE]

where $\epsilon$ is the added noise and $\epsilon\sim\mathcal{N}(0,0.1^{2})$. The training inputs $\boldsymbol{x}_{1}$ and $\boldsymbol{x}_{2}$ are generated from $[-2,2]$, while for testing $\boldsymbol{x}_{1}$ and $\boldsymbol{x}_{2}$ range over $[-2.5,2.5]$.

Rastrigin 3D

The Rastrigin function is highly multimodal, but locations of the minima are regularly distributed. Rastrigin proposed the function in two dimensions (Rastrigin 1974), which was later generalized to $d$ dimensions by Rudolph (Rudolph 1990). However, for the purpose of this work, the function is used in 3 dimensions ($d=3$):

[TABLE]

where $\epsilon\sim\mathcal{N}(0,0.1^{2})$. Each training features $\boldsymbol{x}$ is sampled from $[-5.12,5.12]$, whereas the testing features are generated from $[-5.5,5.5]$.

Griewank 4D

The local minima of the Griewank function are widespread and regularly distributed. It is presented in $d$ dimensions (Griewank 1981):

[TABLE]

where $\epsilon\sim\mathcal{N}(0,0.1^{2})$ is homoscedastic noise. The dimenionality of this function chosen for the aims of this work is $d=4$. The training features $\boldsymbol{x}$ are generated from $[-500,500]$, while the testing feature vectors from $[-600,600]$.

Ackley 7D

Ackley function is introduced as a $d$-dimensional function (Ackley 1987):

[TABLE]

[TABLE]

where $a=20,b=0.2,c=2\pi$ and $\epsilon\sim\mathcal{N}(0,0.1^{2})$. For the purposes of this work, $d=7$, $\boldsymbol{x}$ for training are sampled from $[-30,30]$ and $\boldsymbol{x}$ for testing are generated from $[-32.768,32.768]$, as the latter is the usual test area (Molga and Smutnicki 2005).

Trigonometric

He et al. 1D

The generating function:

[TABLE]

where $\epsilon\sim\mathcal{N}(0,0.1^{2})$, is proposed by He et al. (He, Lakshminarayanan, and Teh 2020). In order to detail uncertainty on out-of-distribution test data, the training points are partitioned into two equal-sized clusters (Figure 5(a)). The clusters, that configure the training data $\boldsymbol{x}$, are generated from $[-2,-0.67]$ and $[0.67,2]$, while the testing points $\boldsymbol{x}$ range in $[-6,6]$.

Forrester et al. 1D

Forrester et al. function is a simple one-dimensional multimodal function (Forrester, Sobester, and Keane 2008):

[TABLE]

where $\epsilon\sim\mathcal{N}(0,0.1^{2})$. Similarly to He et al. dataset, training vector $\boldsymbol{x}$ is split into two clusters to detail uncertainty on OoD test data (Figure 5(b)). Thus, $\boldsymbol{x}$ ranges over $[0.2,0.4]$ and $[0.65,0.85]$, while testing $\boldsymbol{x}$ is sampled from $[0,1]$.

Ishigami 3D

Ishigami function, introduced in (Ishigami and Homma 1990), shows strong non-linearity and non-monotonicity as well as characteristic dependence on $x_{3}$ (Sobol and Levitan 1999):

[TABLE]

where $a=7$ and $b=0.1$, following (Crestaux, Le Maître, and Martinez 2009). Also, $\epsilon\sim\mathcal{N}(0,0.1^{2})$. The training values of $\boldsymbol{x}_{1},\boldsymbol{x}_{2},\boldsymbol{x}_{3}$ are sampled from $[-\frac{\pi}{2},\frac{\pi}{2}]$. Similarly, $\boldsymbol{x}_{1},\boldsymbol{x}_{2},\boldsymbol{x}_{3}$ for testing are sampled $[-\frac{2\pi}{3},\frac{2\pi}{3}]$ respectively.

Friedman 5D

Friedman et al. have proposed the following five-dimensional function (Friedman 1991; Friedman, Grosse, and Stuetzle 1983):

[TABLE]

where $\epsilon\sim\mathcal{N}(0,0.1^{2})$. The training data $\boldsymbol{x}_{1},\boldsymbol{x}_{2},\boldsymbol{x}_{3},\boldsymbol{x}_{4}$ and $\boldsymbol{x}_{5}$ are sampled from $[0,0.5]$, while the testing data $\boldsymbol{x}_{1},\boldsymbol{x}_{2},\boldsymbol{x}_{3},\boldsymbol{x}_{4}$ and $\boldsymbol{x}_{5}$ is generated from $[0,1]$.

Others

Roos & Arnold 5D

The Roos & Arnold function, proposed in (Roos and Arnold 1963), is formed from the products of one-dimensional functions:

[TABLE]

where $\epsilon\sim\mathcal{N}(0,0.1^{2})$ and $d=5$. It is described by Kucherenko et al. as a function with dominant high-order interaction terms and a high effective dimension (Kucherenko et al. 2011). $\boldsymbol{x}$ for training and $\boldsymbol{x}$ for testing are sampled respectively from $[0,0.8]$ and $[0,1]$.

Sum of powers 6D

This bowl-shaped $D$-dimensional function, introduced in (Molga and Smutnicki 2005), represents a sum of different powers:

[TABLE]

where $\epsilon\sim\mathcal{N}(0,0.1^{2})$ and $d=6$. The training data $\boldsymbol{x}_{i}$ ranges over $[-0.75,0.75]$, while the testing data $\boldsymbol{x}_{i}$ from $[-1,1]$.

Styblinski-Tang 9D

Styblinski-Tang function is proposed as a function in $d$ dimensions (Yi, Ahn, and Ji 2020):

[TABLE]

where $\epsilon\sim\mathcal{N}(0,0.1^{2})$ and $d=9$. The testing $\boldsymbol{x}_{i}$ and training $\boldsymbol{x}_{i}$ feature vectors are sampled from the intervals $[-5,5]$ and $[-6,6]$ respectively.

Adapted Welch et al. 10D

The original function, proposed by Welch et al. (Welch et al. 1992), contains 20 dimensions such that some input variables have a very high effect on the output, compared to others. This function is considered challenging, because of its interactions and nonlinear effects. To fit the needs of this work, the Welch et al. function is adapted and its new version has 10 dimension, while still preserving its characteristics:

[TABLE]

[TABLE]

[TABLE]

where $\epsilon\sim\mathcal{N}(0,0.1^{2})$. The ranges of training features $\boldsymbol{x}_{1},\boldsymbol{x}_{2},\boldsymbol{x}_{3},\boldsymbol{x}_{4},\boldsymbol{x}_{5},\boldsymbol{x}_{6},\boldsymbol{x}_{7},\boldsymbol{x}_{8},\boldsymbol{x}_{9}$ and $\boldsymbol{x}_{10}$ are $[-0.5,0.5]$. The testing features $\boldsymbol{x}_{1},\boldsymbol{x}_{2},\boldsymbol{x}_{3},\boldsymbol{x}_{4},\boldsymbol{x}_{5},\boldsymbol{x}_{6},\boldsymbol{x}_{7},\boldsymbol{x}_{8},\boldsymbol{x}_{9}$ and $\boldsymbol{x}_{10}$ and $[-1,1]$.

Appendix C Real-world datasets

The number of training data points and testing data points for each dataset is shown in Table 4.

Boston housing

The goal of the Boston housing dataset is to predict the price of a house given its number of rooms and other context factors. However, in this paper, the number of rooms is the only considered independent variable and context factors, such as house location and crime rate by town, are disregarded.

Abalone shells

The Abalone dataset contains data regarding abalone shells and the regression task is to predict the age of a shell, determined by the number of rings, given several physical measurements (Nash et al. 1994). The features used in this study are: length (denoting the longest shell measurement in mm), diameter in mm, height (with meat in shell in mm), whole weight (whole abalone in grams) and sucked weight (meat weight in grams).

Naval propulsion plant

The naval propulsion plant dataset has been generated from a sophisticated simulator of a Gas Turbine (GT) (Coraddu et al. 2014). The task is to predict the GT propulsion plant’s decay state coefficient. Originally, the features are given in a 16-dimensional feature vector containing the GT measures at steady state of the physical asset, but in this work only GT shaft torque (kN/m), GT rate of revolutions (rpm), high pressure turbine exit temperature (C) and GT exhaust gas pressure (bar), are being used as features.

Forest fire

This dataset concerns forest fire data from the Montesinho natural park of Portugal (Cortez and de Jesus Raimundo Morais 2007). The aim of this dataset is to predict the burnt area given Fine Fuel Moisture Cod (FFMC) index, Duff Moisture Code (DMC) index, Drought Code (DC) index, Initial Spread (ISI) index, temperature in Celsius degrees and relative humidity (in percentage). The full set of attributes of this data set includes also spatial coordinates within the park, day and month, wind and speed, but those are discarded as their addition provides too detailed context information contradicting this study’s goals. Additionally, the dependent variable was first transformed with a $\ln(x+1)$ function, just like in (Cortez and de Jesus Raimundo Morais 2007).

Parkinson’s disease

Parkinson’s disease dataset is composed of a range of biomedical voice measurements from people with Parkinson’s disease (PD) (Little et al. 2007). In total there are 23 features, but only five are being used: NHR and HNR, which are both measures of ratio of noise to tonal components in the voice status, DFA, which is a signal fractal scaling exponent, PPE, denoting three nonlinear measures of fundamental frequency variation and RPDE, which is a nonlinear dynamical complexity measure. Therefore, the goal of this regression task it to predict the total Unified Parkinson’s Disease Rating Scale (UPDRS) given the aforementioned five features.