Learning stochastic eigenvalues

Roman Andreev

► To cite this version:

Roman Andreev. Learning stochastic eigenvalues. 2016. hal-01313404v2

HAL Id: hal-01313404

https://hal.science/hal-01313404v2

Preprint submitted on 12 Dec 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# LEARNING STOCHASTIC EIGENVALUES

ROMAN ANDREEV†

**ABSTRACT.** We train an artificial neural network with one hidden layer on realizations of the first few eigenvalues of a partial differential operator that is parameterized by a vector of independent random variables. The eigenvalues exhibit “crossings” in the high-dimensional parameter space. The training set is constructed by sampling the parameter either at random nodes or at the Smolyak collocation nodes. The performance of the neural network is evaluated empirically on a large random test set. We find that training on random or quasi-random nodes is preferable to the Smolyak nodes. The neural network outperforms the Smolyak interpolation in terms of error bias and variance on nonsimple eigenvalues but not on the simple ones.

## 1. THE STOCHASTIC EIGENVALUE PROBLEM

In this note we compare the out-of-sample prediction of the Smolyak sparse grid interpolant with that of a simple artificial neural network for a stochastic eigenvalue problem. Our model problem is the parametric partial differential eigenvalue problem

$$ (1) \qquad -\nabla \cdot (c(x, y)\nabla u(x, y)) + u(x, y) = \lambda(y)u(x, y), \quad x \in D, \ y \in Y, $$

where $D := (-1, 1)^2$ is the spatial domain and $Y := [0, 1]^M$ is the parameter space. We impose homogeneous Neumann spatial boundary conditions. The conductivity coefficient has the form

$$ (2) \qquad c(x, y) := 1 + \sum_{m=1}^{M} \psi_m(x) y_m. $$

We think of the $y_m \in [0, 1]$ as independent uniformly distributed random variables. Throughout, the stochastic dimension is $M = 16$ and the $\psi_m$ are the indicator functions of the subdomains obtained by partitioning the spatial domain into $4 \times 4$ equal parts. The problem (1) will be discretized with P1 finite elements with $42^2$ unknowns using the Matlab PDE toolbox.

Our quantity of interest is the set of the first 9 eigenvalues,

$$ (3) \qquad \mathrm{QI}(y) := \{1 = \lambda_1(y) \le \lambda_2(y) \le \dots \le \lambda_9(y)\}. $$

For $y = 0$, the nontrivial ones are

$$ \lambda_{2,3} = \frac{0+1}{4}\pi^2 + 1, \quad \lambda_4 = \frac{1+1}{4}\pi^2 + 1, \quad \lambda_{5,6} = \frac{0+4}{4}\pi^2 + 1, \quad \lambda_{7,8} = \frac{1+4}{4}\pi^2 + 1, \quad \lambda_9 = \frac{4+4}{4}\pi^2 + 1, $$
and the last one is simple.

However, as we vary the parameter, the eigenvalues split, leading to a nonsmooth (Lipschitz) dependence on the parameter. For example (figure on the right and Figure 1), interpolating linearly between $c(\cdot, y_A) = 1 + \mathbb{1}_{x_1>0}$ and $c(\cdot, y_B) = 1 + \mathbb{1}_{x_2>0}$, the eigenvalues #2–3 remain almost identical (slight discretization asymmetry), #5–6 switch midway, while #7–8 (and also #9) remain distinct. These eigenvalue “crossings” or “switchings” make it difficult to capture the parameter dependence over the whole high-dimensional parameter space.

<figure><img src="image_2.png" /></figure>

In the next section we discuss two regression methods for that purpose, both based on sampling in the parameter domain: Smolyak interpolation and artificial neural networks.

Date: December 12, 2016.

2010 Mathematics Subject Classification. 15B51, 35R60, 62M45, 62M40, 65N25, 65N30, 65Y05, 65Y20, 68T37.

Key words and phrases. uncertainty quantification, stochastic eigenvalues, neural networks, sparse grids.
## 2. TWO REGRESSION METHODS

2.1. **The Smolyak interpolant** [2]. For each integer level $l \ge 0$ let $i_l$ be a univariate interpolation operator on the interval $[0, 1]$. We will use the polynomial interpolation operator based on the Clenshaw–Curtis nodes

$$ (4) \qquad N_0 = \{\frac{1}{2}\} \quad \text{and} \quad N_l = \{\frac{1}{2}(1 - \cos(2^{-l}k\pi)) : 0 \le k \le 2^l\}, \quad l \ge 1. $$

For a multiindex set $\Lambda \subset \mathbb{N}_0^M$, the multivariate Smolyak interpolation operator $l_\Lambda$ is defined by

$$ (5) \qquad l_\Lambda := \sum_{\nu \in \Lambda} (j_{\nu_1} \otimes j_{\nu_2} \otimes \dots \otimes j_{\nu_M}), $$

where $j_l := i_l - i_{l-1}$ is the univariate increment and $i_{-1} := 0$. We confine the discussion to the multiindex sets of a given total level $L \ge 0$,

$$ (6) \qquad \Lambda_L := \{\nu \in \mathbb{N}_0^M : \sum_m \nu_m \le L\}, $$

and write $l_L$ for the resulting multivariate Smolyak interpolation operator. This multiindex set is monotone: if $\mu \in \Lambda_L$ and if $\nu \le \mu$ coordinatewise then also $\nu \in \Lambda_L$. This, and the nestedness property $N_{l-1} \subset N_l$ of the Clenshaw–Curtis nodes imply the unisolvency of the multivariate interpolation operator on the sparse grid (observe the redundancy in the union)

$$ (7) \qquad \mathcal{N}_L := \bigcup_{\nu \in \Lambda_L} (N_{\nu_1} \times N_{\nu_2} \times \dots \times N_{\nu_M}) \subset Y, $$

because the number of collocation nodes in (7) matches the dimension of the range of the interpolation operator (5). With $M = 16$ stochastic dimensions, the number of collocation nodes is $\#\mathcal{N}_0 = 1$, $\#\mathcal{N}_1 = 33$, $\#\mathcal{N}_2 = 545$, $\#\mathcal{N}_3 = 6'049$, etc. We write

$$ (8) \qquad \mathrm{Sm}[L] := l_L \circ \mathrm{QI} : \mathbb{R}^M \to \mathbb{R}^9 $$

for the interpolated quantity of interest.

2.2. **Artificial neural networks.** We consider the simple class of neural networks of the form

$$ (9) \qquad \text{input } (M \text{ parameters}) \xrightarrow{\text{linear}} \text{ReLU } (H \text{ hidden units}) \xrightarrow{\text{linear}} \text{output } (9 \text{ real values)}, \\ \text{or as the composition} $$

$$ (10) \qquad \mathrm{NN} : \mathbb{R}^M \to \mathbb{R}^9, \quad \mathrm{NN} = (\mathbf{h} \mapsto \mathbf{Ch} + \mathbf{c}) \circ \mathrm{ReLU} \circ (\mathbf{y} \mapsto \mathbf{By} + \mathbf{b}), $$

with the $(9(H+1)+H(M+1))$ learnable parameters

$$ (11) \qquad \{\mathbf{C} \in \mathbb{R}^{9 \times H}, \mathbf{c} \in \mathbb{R}^9\} = \text{2nd all-to-all layer}, \quad \{\mathbf{B} \in \mathbb{R}^{H \times M}, \mathbf{b} \in \mathbb{R}^H\} = \text{1st all-to-all layer.} $$

The “rectified linear unit” is the function $\mathrm{ReLU}: x \mapsto \max\{x, 0\}$, acting componentwise when applied to a vector. The number of “hidden units” will be $H \in \{M, 4M, 16M\}$.

A training set or a test set is a collection of reference nodes $\mathcal{N} \subset Y$ together with the precomputed values of the quantity of interest (3) at those nodes. We refer to a pair $(y, \mathrm{QI}(y))$ as a sample. We consider two types of reference nodes for training the neural network:

(1) the Smolyak collocation nodes $\mathcal{N}_L$ from (7), and

(2) randomly generated nodes $\tilde{\mathcal{N}}_L$ of the same cardinality as $\mathcal{N}_L$.

The random training set is generated only once for each level $L$ and used for all configurations. By “average” we will mean the ensemble average over the training set or over the test set.

The learnable parameters are initialized uniformly at random from $[-s, s]$ where $s = 1/\sqrt{M}$ for the first linear layer and $s = 1/\sqrt{H}$ for the second one. Then, they are trained in 10k iterations of stochastic gradient descent. Each iteration is a pass over the training set in random order; for each sample $(y, q)$, the learnable parameters $\alpha$ are adjusted (all at once) by the gradient scheme

$$ (12) \qquad \alpha \leftarrow \alpha - r \frac{\partial}{\partial \alpha} \mathrm{D}(q, \mathrm{NN}(y)), \quad \mathrm{D}(q, p) := \frac{1}{9} |q - p|_1, $$

where $r > 0$ is the learning rate (initially, $r = 0.1$) and $| \cdot |_1$ is the vector 1-norm measuring the prediction discrepancy. The derivatives of $\mathrm{ReLU}$ and $| \cdot |_1$ are right-continuous. After each iteration, the learning rate is set to 1% of the average prediction discrepancy. The input/output data are offset by their average during the training phase. We write $\mathrm{NN}[H, \mathcal{N}]$ for the neural network with $H$ hidden units trained on the node set $\mathcal{N}$. We use the `torch7` environment for the training.
### 3. EVALUATION

Given a predictor $P = (8)$ or $P = (10)$, we evaluate its quality over a fixed test set of 100k samples on random test nodes $\mathcal{N}_{test} \subset Y$. The average and the standard deviation of the eigenvalues in the test set are as follows (rounded to two postcomma digits):

<table><thead><tr><th>eigenvalue:</th><th>#1</th><th>#2</th><th>#3</th><th>#4</th><th>#5</th><th>#6</th><th>#7</th><th>#8</th><th>#9</th></tr></thead><tbody><tr><td>average</td><td>1</td><td>4.55</td><td>4.68</td><td>8.24</td><td>5.24</td><td>15.47</td><td>18.77</td><td>19.62</td><td>29.6</td></tr><tr><td>std dev</td><td>0</td><td>0.2</td><td>0.2</td><td>0.39</td><td>0.74</td><td>0.73</td><td>0.97</td><td>0.98</td><td>1.47</td></tr></tbody></table>

Predictions along a non-axiparallel segment in the parameter space are shown in Figure 1. We examine the errors

$$ (13) \qquad \delta_k(y) := (\mathrm{QI}(y) - \mathrm{P}(y))_k, \quad y \in \mathcal{N}_{\mathrm{test}}, $$

in the prediction of the k-th eigenvalue. Unsorted predictions occurred at the following rates:

<table><thead><tr><th></th><th>Sm[L]</th><th>NN[M, N<sub>L</sub>]</th><th>NN[4M, N<sub>L</sub>]</th><th>NN[M, Ñ<sub>L</sub>]</th><th>NN[4M, Ñ<sub>L</sub>]</th></tr></thead><tbody><tr><td>L = 2</td><td>71.5%</td><td>3%</td><td>0.6%</td><td>0.31%</td><td>1.3%</td></tr><tr><td>L = 3</td><td>63.5%</td><td>2.8%</td><td>0.15%</td><td>0.02%</td><td>0.01%</td></tr></tbody></table>

We call the average $\bar{\delta}_k$ of $\delta_k$ the prediction bias. Figures 2–3 show the empirical cumulative distribution function of the bias-free prediction error ($\delta_k - \bar{\delta}_k$) along with the prediction bias $\bar{\delta}_k$ for different predictors. We observe the following (recall, $M = 16$ is the stochastic dimension):

(1) For the nonsimple eigenvalues #2, 3, 5, 6, 7, 8, the neural network $NN[4M, Ñ_L]$ trained on random nodes performs consistently better in terms of the prediction bias and prediction error variance than the Smolyak interpolation operator $Sm[L]$.

(2) This behavior is reversed for the simple eigenvalues #4 and #9. This is explained by the analytic dependence of simple eigenvalues on the parameter [1]. Increasing the number of hidden units from $H = 4M$ to $H = 16M$ is still not sufficient to surpass interpolation.

(3) Training the neural network on the Smolyak nodes as opposed to random nodes clearly decreases the prediction accuracy.

(4) The neural network with $M$ hidden units instead of $4M$ is less competitive on training level $L = 3$ than on $L = 2$.

(5) The results are essentially the same for the quasi-random Sobol points (as generated by the Matlab command `sobolset(M)` instead of the random ones (not shown).

Finally, we note that neural networks and sparse grid interpolation have different distributions of offline/online costs. While the offline training of the neural network is lengthy, the online prediction is much quicker than sparse grid interpolation, especially for high-dimensional parameters. However, both effects can be offset by parallelization.

## ACKNOWLEDGMENT

Supported by the Swiss NSF Advanced Postdoc.Mobility grant #164616.

## REFERENCES

[1] Roman Andreev and Christoph Schwab. Sparse tensor approximation of parametric eigenvalue problems. In Ivan G. Graham, Thomas Y. Hou, Omar Lakkis, and Robert Scheichl, editors, *Numerical Analysis of Multiscale Problems*, volume 83 of *Lecture Notes in CSE*, pages 203–241. Springer Berlin Heidelberg, 2012. 3

[2] Sergey A. Smolyak. Quadrature and interpolation formulas for tensor products of certain classes of functions. *Sov. Math. Dokl.*, 4:240–243, 1963. 2

<footnote>†UNIVERSITÉ PARIS DIDEROT, SORBONNE PARIS CITÉ, LJLL (UMR 7598 CNRS), F-75205, PARIS, FRANCE
E-mail address: roman.andreev@upmc.fr</footnote>
<figure><img src="image_3.png" /><figcaption>Above: The eigenvalues #5-6 switch (left). The eigenvalues #7-8 come close (right). The Smolyak interpolant for #7-8 is outside the range shown.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Above: The eigenvalues #2-3 are slightly distinct (due to discretization asymmetry).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Above: The simple eigenvalues #4 and #9. Eigenvalue #9 comes close to #10 midway.</figcaption></figure>

FIGURE 1. Exact and predicted eigenvalues along the parameter space segment from $c(\cdot, y_A) = 1 + 1_{x_1>0}$ to $c(\cdot, y_B) = 1 + 1_{x_2>0}$. See #2 for the annotation.
<figure><img src="image_6.png" /><figcaption>FIGURE 2. The empirical cumulative distribution function of the bias-free prediction errors for the eigenvalues #2, 3, 4, 5 (left to right). The prediction bias is shown as the bar in the lower right corner. Training levels $L = 3$ (top) and $L = 2$ (bottom).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 3. As in Figure 2 for the eigenvalues #6, 7, 8, 9 (left to right).</figcaption></figure>