# Deep neural networks for smooth approximation of physics with higher order
and continuity B-spline base functions

Kamil Doległo(1), Anna Paszyńska(2), Maciej Paszyński(1) and Leszek Demkowicz
(3) (1)Institute of Computer Science,
AGH University of Science and Technology, Kraków, Poland
e-mail<EMAIL_ADDRESS>
(2) Faculty of Physics, Astronomy and Applied Computer Science,
Jagiellonian University, Kraków, Poland
e-mail<EMAIL_ADDRESS>
(3) Oden Institute for Computational and Engineering Sciences,
The University of Texas at Austin, USA
e-mail<EMAIL_ADDRESS>
###### Abstract

This paper deals with the following important research question.
Traditionally, the neural network employs non-linear activation functions
concatenated with linear operators to approximate a given physical phenomenon.
They ”fill the space” with the concatenations of the activation functions and
linear operators and adjust their coefficients to approximate the physical
phenomena. We claim that it is better to ”fill the space” with linear
combinations of smooth higher-order B-splines base functions as employed by
isogeometric analysis and utilize the neural networks to adjust the
coefficients of linear combinations. In other words, the possibilities of
using neural networks for approximating the B-spline base functions’
coefficients and by approximating the solution directly are evaluated. Solving
differential equations with neural networks has been proposed by Maziar Raissi
et al. in 2017 [16] by introducing Physics-informed Neural Networks (PINN),
which naturally encode underlying physical laws as prior information.
Approximation of coefficients using a function as an input leverages the well-
known capability of neural networks being universal function approximators
[6]. In essence, in the PINN approach the network approximates the value of
the given field at a given point. We present an alternative approach, where
the physcial quantity is approximated as a linear combination of smooth
B-spline basis functions, and the neural network approximates the coefficients
of B-splines. This research compares results from the DNN approximating the
coefficients of the linear combination of B-spline basis functions, with the
DNN approximating the solution directly. We show that our approach is cheaper
and more accurate when approximating smooth physical fields.

###### keywords:

deep neural networks , physics informed neural networks , isogeometric
analysis , finite element method

##  1 Introduction

Isogeometric analysis (IGA) has been proposed in 2005 by T.J.R.Hughes et al.
[8] as a generalization of the Finite Element Method into higher order and
continuity B-spline basis functions. The method employs smooth B-spline basis
to approximate scalar or vector fields representing solutions of different
physical phenomena. The implementation aspects of IGA are summarized in [12].
There are a couple of open source libraries offering the isogeometric analysis
computational framework. One of them is PetIGA. PetIGA, a framework for high-
performance isogeometric analysis has been developed and described in an
article by L. Dalcin, et. al.[4]. The framework is based on PETSc
111https://www.mcs.anl.gov/petsc/, a high-performance library for the scalable
solution of partial differential equations. This library uses traditional
solvers for ordinary and partial differential equations. Another open source
package is GeoPDEs222http://rafavzqz.github.io/geopdes/, a package for
isogeometric analysis in MATLAB and Octave. There are also some MATLAB
implementations of stabilized finite element method computations available in
an e-book by M. Paszyński[14]. Approximation of physics-based fields with
splines requires selecting appropriate coefficients for the base functions. In
IGA the coefficients are obtained by solving appropriate systems of linear
equations.

The usage of neural networks for solving differential equations has been
proposed many times. An article by C. Michoski et al.[10] presents work on
solving partial differential equations with deep neural networks, reviews and
extends some of them and focuses on irregular solutions.

Deep neural networks have recently been shown [2] to solve ordinary
differential equations with an adaptive precision to speed ratio. In the paper
by Yulia Rubinova et. al., neural networks are used to parameterize the
derivative of the hidden state instead of the traditional approach of
specifying a discrete sequence of hidden layers. In other words, an ODE is
embedded into a neural network, by using differential equations solvers as the
layers and hidden layers are not predetermined beforehand, but rather the
number of layers depends on the desired accuracy. This new family of deep
neural network models could possibly replace residual networks. ODE-nets are
suitable for time-series data. While ordinary neural networks are discrete and
have problems with irregular data, ODE-nets are continuous and allow for
evaluation at any point, which results in better accuracy for time-series
data. They also have faster testing times at the expense of longer training
times. Deep Neural Network are also used for optimization of the computational
procedures of finite element method [1, 7].

An article by M. Raissi et al.[16] introduced Physics-informed Neural
Networks, which naturally encode underlying physical laws as prior information
and can be used for both continuous and discrete time models. They are
constructed with the help of automatic differentiation to differentiate neural
networks with respect to their input coordinates and model parameters. PINNs
is used for some frequently modelled physical phenomena, like heat transport
or flow simulation [13, 15].

Application of Physics-informed Neural Networks in electric power systems has
been recently researched in a preprint by George S. Misyris et. al.[11]. The
authors state that the usage of PINNs allows to accurately determine results
of differential equations up to 87 times faster than conventional methods. The
method used requires less initial training data and can result in smaller
neural networks while demonstrating high performance.

Physics-informed neural networks sometimes fail to be trained, however some
research on that topic is being conducted. In the preprint by S. Wang, et al.
[17] a novel gradient descent algorithm is being proposed to improve PINNs.
One of the PINN models’ disadvantages over regular deep neural networks, is
that they can only predict one specific solution of a Partial Differential
Equation. This may be possible to overcome by adapting the PINN learning
method to accommodate additional PDE solutions. Nevertheless, they could be
trained with minimal to no data from the actual solution if the appropriate
boundary conditions were given.

Approximation of coefficients using a function as an input leverages the well-
known capability of neural networks being universal function approximators[6].
This work establishes that the standard multilayer feedforward network
architectures can approximate virtually any function of interest to any
desired degree of accuracy, provided sufficiently many hidden units are
available. It does not, however, address the issue of how many units are
needed to attain a given accuracy of approximation. For this paper, it means
that a possible failure of the approach with approximating field values could
be attributed not to a misuse of neural networks, but rather to inadequate
learning, inadequate numbers of hidden units, or the presence of a stochastic
rather than a deterministic relation between the input and the target.

This paper proposes a new method to reduce the computational cost of the
approximation by applying Deep Neural Networks (DNN) to approximate the base
functions’ coefficients. We compare this method to the DNN approximating the
solution obtained directly from IGA solver. We also compare the method to
Physics Informed Neural Networks (PINN) [16] approach where we approximate the
field values directly.

There are the following open research questions:

  1. 1.

What is the approximation error of similar DNN in both approaches?

The approximation error may differ for the same problem solved with different
approaches. One method might tend to be more accurate than the other methods.

  2. 2.

What is the approximation errors difference between the IGA solution and the
solutions from DNN based approaches?

If applying DNN in IGA problems proves to be faster than the traditional IGA
method, it might still be unusable in practice, because of high approximation
errors.

  3. 3.

Which approach is better?

The Universal Approximation Theorem[6] states that neural networks can
represent a variety of functions when given appropriate weights, it does not,
however, provide a way to construct those networks and weights, only stating
that such a construction is possible.

There are the following research objectives:

  1. 1.

Investigating the viability of using DNN for approximating coefficients of
base B-spline functions used in IGA.

Neural networks may prove to be usable in IGA. DNN based approaches should be
inspected in terms of accuracy, performance and model size.

  2. 2.

Comparing the approaches in applying DNN in terms of accuracy and training
time.

Direct solution approximation with DNN, B-spline coefficient approximation
using DNN, and PINN approximations may differ in accuracy and training time,
due to inherent differences in the nature of all approaches.

Paper [5] presents an application of PINN into solid mechanics examples, and
it compares the PINN model learning based on the exact solution, on the finite
element method solution and on the isogeometric analysis solution. It
concludes that IGA smooth model results in superior convergence of the PINN
training. In our paper, we show that DNN can actually learn the coefficients
of the B-spline basis functions, and it results in a faster convergence and
more accurate approximation of the solution than either standard PINN or DNN
learning IGA solution.

The structure of our paper is the following. We start in Section 2 by
introducing the coefficients of the linear combination of B-splines and their
applicability to approximate smooth physical fields. Next, Section 3 compares
the three methods employed to incorporate the DNN for the solution of PDEs.
Section 4 presents the discussion on the conclusions and future work. We also
introduce three Appendixes where we derive the three compared methods on a
simple one-dimensional example.

##  2 Solving differential equations with neural networks and linear
combinations of higher-order and continuity B-spline base functions

Let assume we want to approximate a function $f(x,y)$ with a linear
combination of B-spline functions $u(x,y)$:

$f(x,y)\approx u(x,y)$
$u(x,y)=\sum_{i=1,j=1}^{N_{x},N_{y}}u_{i,j}B_{i}^{x}B_{j}^{y}$ (1)

where $B_{i}^{x}$, $B_{j}^{y}$ denote basis functions over the $x$ and $y$
axis respectively, while $N_{x}$ and $N_{y}$ denote numbers of basis
functions. An example 2D quadratic base is shown in Figure 1. The coefficients
of the best approximation can be obtained by solving

$\displaystyle\begin{bmatrix}\int
B_{1}^{x}B_{1}^{y}B_{1}^{x}B_{1}^{y}dxdy&\ldots&\int
B_{1}^{x}B_{1}^{y}B_{N_{x}}^{x}B_{N_{y}}^{y}dxdy\\\ \vdots&\vdots&\vdots\\\
\int B_{N_{x}}^{x}B_{N_{y}}^{y}B_{1}^{x}B_{1}^{y}dxdy&\ldots&\int
B_{N_{x}}^{x}B_{N_{y}}^{y}B_{N_{x}}^{x}B_{N_{y}}^{y}dxdy\\\ \end{bmatrix}$
$\displaystyle\begin{bmatrix}u_{1,1}\\\ \vdots\\\
u_{N_{x},N_{y}}\end{bmatrix}=\begin{bmatrix}\int
f(x,y)B_{1}^{x}B_{1}^{y}dxdy\\\ \vdots\\\ \int
f(x,y)B_{N_{x}}^{x}B_{N_{y}}^{y}dxdy\\\ \end{bmatrix}$ (2)

This system is solvable using traditional solvers but computation-heavy [3],
despite symmetries in the matrix, especially for a case of non-regular
geometry of the computational domain where the B-spline basis functions are
span.

Figure 1: An example of 2-dimensional quadratic basis generated from
$[0,0,0,1,2,2,2]$ knot vectors over $x$ and $y$ axes, and the smooth linear
combination of them.

The main goal of this project is to verify three different approaches for
training deep neural networks in order to provide solutions of IGA simulations
with partial differential equations.

  * 1.

Approximation of coefficients The first method is based on ”indirect”
approximation of the solution, where a DNN predicts the coefficients of the
B-spline basis functions that span over the computational mesh. Linear
combinations of those basis functions approximate the solution field. In this
approach, neural networks approximate the values of the base functions’
coefficients, $ANN_{i,j}(n)\approx u_{i,j}$, of the linear combination
$u_{h}(x,y)=\sum_{i=1,...,N_{x};j=1,...,N_{y}}u_{i,j}B_{i}^{x}B_{j}^{y}$
(where $B_{i}^{x}$ and $B_{j}^{y}$ are the B-spline base functions). This
linear combination approximates the exact solution
$u_{h}(x,y)\approx~{}u(x,y)$ of a given problem. The final solution $u_{h}$ is
obtained by multiplying the relevant base functions by the coefficients
predicted by the neural network. There may exist some correlation between the
coefficients, that the neural networks can learn, especially in the case of a
family of functions. The details of this method for simple one-dimensional
example are presented in Appendix A.

  * 2.

Direct approximation of the result The second method is based on direct
approximation of the solution scalar field by a DNN. Here we have a neural
network approximating the solution directly, $ANN(n,x,y)\approx u(x,y)$, we do
not use the intermediate linear combination $u_{h}(x,y)$ in contrast to the
previous case. The details of this method for simple one-dimensional example
are presented in Appendix B.

  * 3.

Physics Informed Neural Network The third method is based on classical PINN
approach [16]. The details of this method for simple one-dimensional example
are presented in Appendix C.

The input data for the DNNs in the first and second cases are the physical
model parameters (eg. boundary conditions, material data and forcing), but the
outputs from the DNNs black boxes are different in both cases. In the first
case, the output is the coefficients of the linear combination of B-spline
basis functions, which depend on the dimensions of the knot vectors and the
order and continuity of the basis functions. In the second case, it is the set
of point-wise values of the solution scalar field. The output is thus
parameterized with the number of point values along the X and Y-axes.

DNN for the first and second cases learned a family of functions. Physics-
informedneural networks can learn only a selected function from the family,
due to the training method used. The DNN for PINN takes as an input the $x$
and $y$ coordinates of the point, and the $n$ parameter is fixed for the PINN
neural network.

##  3 Comparison of three methods

Figure 2: Visualization of the heat transfer problem. The red (external) edge
has a Neumann boundary condition $\frac{\partial u}{\partial{\bf v}}(x)=g(x)$
and the black (internal) edge has a Dirichlet boundary condition $u=0$, where
$u$ is the heat function, ${\bf v}$ is the versor normal to the boundary, and
$g(x)$ is a given function.

Heat transfer experiments have been conducted for two-dimensional heat
equation on an L-shaped area. External sides are bound with a Neumann boundary
condition and internal sides are bound with a Dirichlet boundary condition.
This models the heating of the external boundary and fixed zero temperature at
the internal boundary.

The equation that we solve is $\Delta u=0$ where $\Delta$ is the Laplacian,
namely $\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial
y^{2}}=0$ the forcing term is equal to 0, the Dirichlet boundary condition is
equal to $0$, and the Neumann boundary condition is $\frac{\partial
u}{\partial{\bf v}}=g$, where we compute the directional derivative in the
direction ${\bf v}$ perpendicular to the boundary, and $g(x,y)$ is a given
”heating” function. We enforce the parameterized heating of the external
boundary by using a family of ”heating” functions, namely

$\displaystyle g(x,y)=v_{i}2\pi ncos(2\pi n\cdot x_{i})\cdot sin(2\pi n\cdot
x_{j})\textrm{ for }|x_{i}|=1$

where $i=1,2$, $[j=(i+1)mod1]+1$, and ${\bf v}=(v_{1},v_{2})$ is the versor
normal to the boundary. By changing the $n$ parameter we adjust the heating
function $g$ and we obtain different solutions to the heat transfer problem.

Following [14] (chapter 3), we transform the problem into the weak
formulation. We introduce the domain $\Omega=[-1,1]^{2}$, and we test and
integrate by parts to obtain

$\displaystyle b(u,v)=l(v)\quad\forall v\in V$ $\displaystyle
b(u,v)=\int_{\Omega}\nabla u\cdot\nabla vdxdy$ $\displaystyle
l(v)=\int_{\partial\Omega}gvdS$ (3)

We discretize with B-spline base functions. We introduce the knot vectors
along $x$ and $y$ axis of coordinates, where we repeat the knot at 0, to
obtain the tensor product basis function $B^{x}_{i}B^{y}_{j}$,
$i=1,...,N_{x}$, $j=1,...,N_{y}$. An example of knot vectors [0 0 0 0.5 1 1 1
1.5 2 2 2] $\times$ [0 0 0 0.5 1 1 1 1.5 2 2 2] and resulting B-spline base
functions with B-splines over lower-left corner of the domain set to 0, are
illustrated in Figure 3.

Figure 3: Knot vectors [0 0 0 0.5 1 1 1 1.5 2 2 2] $\times$ [0 0 0 0.5 1 1 1
1.5 2 2 2] and resulting B-spline base functions, with the splines over left-
lower quarter of the domain set to 0.

The coefficients of B-splines employed for training of the DNN are obtained by
solving the system of equations

$\displaystyle\begin{bmatrix}b\left(B_{1}^{x}B_{1}^{y},B_{1}^{x}B_{1}^{y}\right)&\ldots&b\left(B_{1}^{x}B_{1}^{y},B_{N_{x}}^{x}B_{N_{y}}^{y}\right)\\\
\vdots&\vdots&\vdots\\\
b\left(B_{N_{x}}^{x}B_{N_{y}}^{y},B_{1}^{x}B_{1}^{y}\right)&\ldots&b\left(B_{N_{x}}^{x}B_{N_{y}}^{y},B_{N_{x}}^{x}B_{N_{y}}^{y}\right)\\\
\end{bmatrix}$ $\displaystyle\begin{bmatrix}u_{1,1}\\\ \vdots\\\
u_{N_{x},N_{y}}\end{bmatrix}=\begin{bmatrix}l\left(B_{1}^{x}B_{1}^{y}\right)\\\
\vdots\\\ l\left(B_{N_{x}}^{x}B_{N_{y}}^{y}\right)\\\ \end{bmatrix}$ (4)

Additionally, for the rows related to B-splines located in the left-lower
quarter of the domain, we set the rows to zero, we put 1.0 on the diagonal and
zero on the right-hand side. An alternative way of setting boundary conditions
would be to add an additional non-trainable layer to the neural network
enforcing the Dirichlet boundary condition following the ideas described in
[9].

###  3.1 Coefficient approximation by DNN

We’d like the network to learn a family of solutions of the heat transfer
problem. The input to the network is the $n$ parameter. The output to the
network are the coefficients $u_{ij}$ of the linear combination of B-splines.

Figure 4: Visualization of the neural network architecture. The value of $i$
is between 100 and 1000, and the value of $j$ is equal to the number of
coefficients, as explained in Table 1.

The network architecture varies with input data - the output layer has to have
the same number of neurons as coefficients. This is the main drawback of this
approach, as the network has to be reconfigured with each mesh change. The
model was trained with a learning rate reduction on plateau.

Layer | Number of neurons | Activation
---|---|---
|  | function
input | 100 - 1000 | ReLU
hidden layer 1 | 100 - 1000 | ReLU
output 2 | equal to | none
| the number of coefficients |
Table 1: Network architecture of the DNN used to approximate the coefficients
of the heat transfer solution. The number of neurons in the hidden layer
varied with inputs. For 169 output coefficients the hidden layer had 100
neurons and the number rose to as much as 1000 neurons for 2025 output
coefficients. ReLU stands for Rectified Linear Unit.

Figure 5: Approximating the $B_{5,2}(x)B_{4,2}(y)$ and $B_{7,2}(x)B_{7,2}(y)$
coefficients with a neural network trained to approximate the coefficient
value directly

We present results obtained for 7 quadratic B-spline functions in each
direction spanned over the $[0,0,0,1,2,2,3,4,4,4]$ knot vector, which gives
$7\times 7=49$ coefficients in total. Coefficient approximation yielded
satisfactory results, see left panel in Figure 6.

Figure 6: Left panel: The IGA solution of (4) (green-blue ) and the
coefficients approximating DNN (orange-red) solutions. The solutions are
practically identical. The MSE of the coefficients is $2.27\text{e-}8$, the
pointwise MSE is $6.60\text{e-}9$. Right panel: The IGA solution of (4) (blue-
green) and the DNN approximating the solution directly (red-yellow). The
solutions are practically identical. The mean squared error is approximately
$8.42\text{e-}7$. Figure 7: Difference between the solutions constructed using
coefficients from IGA solver (4) and the coefficients obtained from trained
DNN. The mean squared error is approximately $6.60\text{e-}9$. Figure 8: Mean
squared error loss function averaged over an epoch for the first 800 out of 2
000 epochs. Each batch contained 19 entries.

The neural network could learn a family of heat functions and accurately
predict the coefficients for the base functions, as shown in Figure 5. The
difference between the solution obtained with IGA solver (4) and the solution
obtained from the DNN trained to approximate the coefficients is shown in the
Figure 7. Model convergence is visible in the Figure 8. Table 2 summarizes the
results for several different meshes and different base functions.

Mesh | # coeffs | Spline degree | MSE value | Training
---|---|---|---|---
|  |  |  | time [s]
10 x 10 | 121 | linear | 1.98e-07 | 9
10 x 10 | 169 | quadratic | 5.29e-07 | 10
10 x 10 | 225 | cubic | 2.18e-07 | 11
20 x 20 | 441 | linear | 6.99e-07 | 14
20 x 20 | 529 | quadratic | 6.96e-07 | 18
20 x 20 | 625 | cubic | 9.14e-08 | 19
40 x 40 | 1681 | linear | 1.46e-06 | 15
40 x 40 | 1849 | quadratic | 1.37e-06 | 16
40 x 40 | 2025 | cubic | 7.58e-06 | 18
Table 2: Summary of obtained training results for a coefficient-approximating
DNN. The training time rises with the number of coefficients, as expected. The
training times are only an indication of the time required to train the model,
because they vary significantly with batch size and the number of iterations.
The error tends to increase with the number of coefficients, however the
actual output is still perfectly usable.

###  3.2 DNN approximating the solution directly

A deep neural network with 2 hidden layers has been trained to directly
approximate the result. The network was designed to take 3 arguments - the
value of $n$, and $x$ and $y$ coordinates of the desired point of the
solution. The model was a fully-connected feed-forward neural network with 2
hidden layers, 100 neurons each. Direct approximation also yielded
satisfactory results. The neural network could learn a family of heat
functions and accurately predict the solutions, as shown in right panel in
Figure 6. The solution is slightly less accurate than the one obtained with
coefficient approximation, however it should still be usable. A difference
between the solution obtained from IGA solver (4) and the DNN approximating
the solution directly is shown in Figure 10. Model convergence is visible in
the Figure 11. Table 4 summarizes the learning times for several different
meshes and different base functions and compares them to the learning times of
the $n$ input network. It is worth mentioning here that the network
approximating the solution does not have to be retrained for different meshes,
in contrast to the coefficient approximation, which has to be retrained for
different meshes, so the training times here should be treated only as an
indication of the order of magnitude.

Layer | Number of neurons | Activation function
---|---|---
input | 100 | ReLU
hidden layer 1 | 100 | ReLU
hidden layer 2 | 100 | ReLU
output | 1 | none
Table 3: Summary of the layers and activation functions of the DNN used to directly approximate the heat transfer solution. ReLU stands for Rectified Linear Unit. Figure 9: Visualization of the architecture of the DNN used to directly approximate the heat transfer solution. Figure 10: Difference between the solution obtained with IGA solver (4) and the the DNN approximating the solution directly. The mean squared error is approximately $8.42\text{e-}7$. Figure 11: Mean squared error loss function over 40 training epochs. Mesh size | Spline degree | Coefficient | Direct
---|---|---|---
|  | network | network
|  | training | training
|  | time [s] | time [s]
10 x 10 | linear | 9 | 115
10 x 10 | quadratic | 10 | 114
10 x 10 | cubic | 11 | 48
20 x 20 | linear | 14 | 43
20 x 20 | quadratic | 18 | 47
20 x 20 | cubic | 19 | 42
40 x 40 | linear | 15 | 44
40 x 40 | quadratic | 16 | 42
40 x 40 | cubic | 18 | 48
Table 4: Comparison of training times of the coefficient-approximating neural
network and the solution approximating neural network.The MSE of achieved
outputs had the same order of magnitude.

###  3.3 Direct approximation with PINN

Deep neural networks in previous sections learned a family of functions.
Physics-informed neural networks can learn only a selected function from the
family, due to the training method used. The network was designed to take 2
arguments - the value of $x$ and $y$ coordinates of the desired point of the
solution. The model was a fully-connected feed-forward neural network with 2
hidden layers, 50 neurons each:

Layer | Number of neurons | Activation function
---|---|---
input | 50 | ReLU
hidden layer 1 | 50 | ReLU
hidden layer 2 | 50 | ReLU
output | 1 | none
Table 5: Summary of the layers and activation functions of the PINN used to
directly approximate the heat transfer solution. ReLU stands for Rectified
Linear Unit. Figure 12: Visualization of the architecture of the PINN used to
directly approximate the heat transfer solution.

A physics-informed neural network has been trained to directly approximate the
result, respecting the heat transport differential equation, the heat source
function, and boundary conditions presented in subsection 3.1. A difference
between the solution obtained with the IGA solver 4 and the solution obtained
from the PINN is shown in Figure 13. The results are satisfactory,
unfortunately PINNs can’t learn families of functions, at least not in their
original form.

###  3.4 A short note on overfitting and network capacity

Despite the loss function plots suggesting the overfitting of models, all
networks performed well on never–seen-before data with MSE values of the same
order of magnitude as in the training data. Coefficient approximating models
were also evaluated with a 5-fold cross validation. While the network
architectures might suggest reduced network capacity compared to the
architectures with the same number of input nodes as the number of input
parameters, adding dropout layers hasn’t improved the accuracy.

Figure 13: Visualization of the difference between the solution obtained with
classical solution and the solution from the PINN. The mean squared error is
approximately $2.43\text{e-}6$. Figure 14: Mean squared error loss function
over the 15 training epochs. Each batch contained 40 entries.

##  4 Conclusions and future work

In terms of research objectives, the results of this paper are as follows:

  1. 1.

Investigating the viability of using NN for approximating coefficients of base
B-spline functions used in IGA.

Application of neural networks in approximating coefficients of the base
functions proved to be possible. Deep Neural Networks were able to learn a
family of functions with more than acceptable accuracy for CAD applications
such as prototyping. More research is needed to check whether neural networks
can accurately approximate coefficients with arbitrary knot vectors as an
input, as opposed to a predefined knot vectors used in this paper.

  2. 2.

Comparing the two approaches in applying NNs to IGA in terms of accuracy and
training time.

Approximating base functions’ coefficients was a more accurate (figures 6, 7,
10, and 13) and quicker (figure 8, 11 and 14) to train approach, but the
network has to be retrained for different meshes. Direct approximation with
deep neural networks was slower to train (table 4) but was independent of the
underlying mesh, as the networks approximate the final solution. Physics-
informed Neural Networks could only predict one specific function, but it
required less data from solved problems, provided all the underlying equations
and boundary conditions are known at training time.

Further research is needed in some areas. The question whether the number of
B-splines and coefficients along with the size of the network impact the
approximation error is left unanswered. The data summarized in the table 2
suggest that there might exist a correlation between the problem size, the
desired accuracy and the approximation error, due to the more difficult
learning process which involves more variables.

## Acknowledgments

The paper was partially financed by AGH University of Science and Technology
Statutory Fund.

##  5 Appendix A: One-dimensional example of neural network learning
coefficients of B-splines

###  5.1 One dimensional heat-transfer problem

Let us introduce the knot vector [0 0 0 1 1 1] defining the quadratic B-spline
basis functions with $C^{0}$ separators

$B_{1,2}(x)=(1-x)^{2};\quad B_{2,2}(x)=2x(1-x);\quad B_{3,2}(x)=x^{2}$ (5)

Figure 15: Three B-splines over a single interval (element)

Let us introduce the problem

$-u^{\prime\prime}(x)=f(x)\quad x\in(0,0.5)$ (6)

defined over $x\in(0,0.5)$, with boundary conditions $u(0)=0$ and
$u^{\prime}(0.5)=g(x)$. We setup $g(x)=n\pi cos(n\pi x)$ and
$f(x)=n^{2}\pi^{2}sin(n\pi x)$. The family of solution of this problem are

$f_{n}(x)=sin(n\pi x)$ (7)

Figure 16: Plot of different solutions $f_{n}(x)$ for $n=0.1,0.2,...,1$.

We transform this problem into the weak form

$\displaystyle\int_{0}^{0.5}u^{\prime}(x)v^{\prime}(x)dx=\int_{0}^{0.5}f(x)v(x)dx+v(0.5)g(0.5)\quad\forall
v$ (8)

and we discretize with B-spline basis functions

$\displaystyle u_{h}=\sum_{i=1,2,3}u_{i}B_{i,2}(x)$ (9)

to obtain

$\displaystyle\begin{bmatrix}\int_{0,1/2}B_{1,2}^{\prime}(x)B_{1,2}^{\prime}(x)dx&\int_{0,1/2}B_{1,2}^{\prime}(x)B_{2,2}^{\prime}(x)dx&\int_{0,1/2}B_{1,2}^{\prime}(x)B_{3,2}^{\prime}(x)dx\\\
\int_{0,1/2}B_{2,2}^{\prime}(x)B_{1,2}^{\prime}(x)dx&\int_{0,1/2}B_{2,2}^{\prime}(x)B_{2,2}^{\prime}(x)dx&\int_{0,1/2}B_{2,2}^{\prime}(x)B_{3,2}^{\prime}(x)dx\\\
\int_{0,1/2}B_{3,2}^{\prime}(x)B_{1,2}^{\prime}(x)dx&\int_{0,1/2}B_{3,2}^{\prime}(x)B_{2,2}^{\prime}(x)dx&\int_{0,1/2}B_{3,2}^{\prime}(x)B_{3,2}^{\prime}(x)dx\\\
\end{bmatrix}\begin{bmatrix}u_{1}\\\ u_{2}\\\ u_{3}\end{bmatrix}=$
$\displaystyle\begin{bmatrix}\int_{0,1/2}B_{1,2}(x)f_{n}(x)dx\\\
\int_{0,1/2}B_{2,2}(x)f_{n}(x)dx\\\ \int_{0,1/2}B_{3,2}(x)f_{n}(x)dx+n\pi
cos(n\pi 0.5)\end{bmatrix}$ (10)

Figure 17: Coefficients of approximation
$u_{1}B_{1,2}(x)+u_{2}B_{2,2}+u_{3}B_{3,2}$ for $n\in(0,0.5)\subset{\cal R}$.

###  5.2 Artificial neural network for $u_{i}$

Let us introduce the artificial neural network

$ANN_{i}(n)=u_{i}$ (11)

where $n$ is the index of the $f_{n}$ function, $i=1,2,3$ (for three
coefficients of B-splines).

Given sinus family function index $n$, it returns the coefficient $u_{i}$ of
B-splines for approximation of this function over $(0,0.5)$.

$ANN_{i}(n)=c_{i}\sigma\left(a_{i}n+b_{i}\right)+d_{i}$ (12)

where the activation function

$\sigma(x)=\frac{1}{1+e^{-x}}$ (13)

###  5.3 Training

The goal of the training is to find values of the weights
$a_{i},b_{i},c_{i},d_{i}$

We prepare a set of samples

  * 1.

We randomly select $n\in(0,1)$

  * 2.

We solve the IGA problem (29)

  * 3.

Input data $(n)$, output data $(u_{1},u_{2},u_{3})$

In other words, we train the neural network for some selected functions
$f_{n}$, with a hope, that it will work for a given function of interest from
the family.

How to train the artificial neural network? We define the error function

$\displaystyle
e_{i}(n)=0.5\left(ANN_{i}(n)-u_{i}(n)\right)^{2}=0.5\left(c_{i}\sigma\left(a_{i}n+b_{i}\right)+d_{i}-u_{i}(n)\right)^{2}=$
(14) $\displaystyle
0.5\left(\left(\frac{c_{i}}{1+exp(-\left(a_{i}n+b_{i}\right)}+d_{i}\right)-u_{i}(n)\right)^{2}$
(15)

Now, we compute the derivatives

$\displaystyle\frac{\partial e_{i}(n)}{\partial
a_{i}}=\frac{c_{i}nexp(-a_{i}n-b_{i})(ANN_{i}(n)-u_{i}(n))}{(exp(-ax-b)+1)^{2}}$
(16) $\displaystyle\frac{\partial e_{i}(n)}{\partial
b_{i}}=\frac{c_{i}exp(-a_{i}n-b_{i})(ANN_{i}(n)-u_{i}(n))}{(exp(-ax-b)+1)^{2}}$
(17) $\displaystyle\frac{\partial e_{i}(n)}{\partial
c_{i}}=\frac{(ANN_{i}(n)-u_{i}(n))}{(exp(-ax-b)+1)}$ (18)
$\displaystyle\frac{\partial e_{i}(n)}{\partial d_{i}}=(ANN_{i}(n)-u_{i}(n))$
(19)

they say “how fast the error is changing if I modify a given coefficient”.

We loop through the data set $\\{n,(u_{1}(n),u_{2}(n),u_{3}(n))\\}_{n\in A}$
where $A$ is the set of selected points from $(0,0.5)$, and we train each of
the three $ANN_{1}$, $ANN_{2}$, and $ANN_{3}$

  1. 1.

Select $(n,(u_{1},u_{2},u_{3}))$

  2. 2.

Compute $u_{i}=ANN_{i}(n)=c_{i}\sigma\left(a_{i}n+b_{i}\right)+d_{i}$

  3. 3.

Compute $e_{i}(n)$

  4. 4.

Compute $\frac{\partial e_{i}(n)}{\partial a_{i}},\frac{\partial
e_{i}(n)}{\partial b_{i}},\frac{\partial e_{i}(n)}{\partial
c_{i}},\frac{\partial e_{i}(n)}{\partial d_{i}}$

  5. 5.

Correct

$\displaystyle a_{i}=a_{i}-\eta*\frac{\partial e_{i}(n)}{\partial a_{i}}$ (21)
$\displaystyle b_{i}=b_{i}-\eta*\frac{\partial e_{i}(n)}{\partial b_{i}}$ (22)
$\displaystyle c_{i}=c_{i}-\eta*\frac{\partial e_{i}(n)}{\partial c_{i}}$ (23)
$\displaystyle d_{i}=d_{i}-\eta*\frac{\partial e_{i}(n)}{\partial d_{i}}$ (24)

where $\eta\in(0,1)$. This is like a local gradient method.

###  5.4 MATLAB implementation

% Creation of dataset
A = [1/5 1/10 1/30; 1/10 2/15 1/10; 1/30 1/10 1/5];
i=1;
for n=0.01:0.01:0.5
rhs= [ (pi*pi*n*n+2*cos(pi*n)-2)/(pi*pi*pi*n*n*n);
(-2*pi*n*sin(pi*n)-4*cos(pi*n)+4)/(pi*pi*pi*n*n*n);
((2-pi*pi*n*n)*cos(pi*n)+2*pi*n*sin(pi*n)-2)/(pi*pi*pi*n*n*n) ];
u=A $\backslash$ rhs;
dataset_in(i)=n;
dataset_u1(i)=u(1);
dataset_u2(i)=u(2);
dataset_u3(i)=u(3);
i=i+1;
endfor
ndataset=i-1;

% Training
a1=1.0; b1=1.0; c1=1.0; d1=1.0;
eta1=0.1;
r = 0 + (1-0).*rand(ndataset,1);
r=r.*ndataset;
for j=1:ndataset
i=floor(r(j));
eval1 = c1*1.0/(1.0+exp(-(a1*dataset_in(i)+b1)))+d1;
error1 = 0.5*(eval1-dataset_u1(i))2;
derrorda = c1*dataset_in(i)*exp(-a1*dataset_in(i)-b1)*
(eval1-dataset_u1(i))/(exp(-a1*dataset_in(i)-b1)+1)2;
a1=a1-eta1* derrorda;
derrordb = c1*exp(-a1*dataset_in(i)-b1)*
(eval1-dataset_u1(i))/(exp(-a1*dataset_in(i)-b1); b1=b1-eta1* derrordb;
derrordc = (eval1-dataset_u1(i))/(exp(-a1*dataset_in(i)-b1)+1);
c1=c1-eta1* derrordc;
derrordd = (eval1-dataset_u1(i));
d1=d1-eta1* derrordd;

% evaluation of ANN approximation of sin(n*pi*x) for n=0.333
n=0.333;
u1 = c1*1.0/(1.0+exp(-(a1*n+b1)))+d1;
u2 = c2*1.0/(1.0+exp(-(a2*n+b2)))+d2;
u3 = c3*1.0/(1.0+exp(-(a3*n+b3)))+d3;
x=0:0.01:0.5;
y=sin(n*pi.*x);
z=u1*(1-x).2+u2*2*x.*(1-x)+u3*x.2;
plot(x,y,x,z);

We tried starting points $1.0,10.0,-1.0,-10.0$ for all the combinations of
$a_{i},b_{i},c_{i},d_{i}$ (256 runs) and the best result (smaller errors) we
obtain for

$a_{1}=b_{1}=c_{1}=d_{1}=1.0$; $a_{2}=b_{2}=1,c_{2}=10.0,d_{2}=-1.0$
$a_{3}=b_{2}=3,c_{3}=10.0,d_{3}=-1.0$

We used $\eta=0.1$. We coded the ANN and the training in hand-made MATLAB
code.

Figure 18: Training for the simple artificial nerual network (11) starting
from $a_{1}=b_{1}=c_{1}=d_{1}=1.0$; $a_{2}=b_{2}=1,c_{2}=10.0,d_{2}=-1.0$,
$a_{3}=b_{2}=3,c_{3}=10.0,d_{3}=-1.0$, for $\eta=0.1$. Figure 19: Convergence
of errors for training of ANN1, ANN2, ANN3

###  5.5 Verification

Figure 20: Verification of the neural network approximation of solution for
$n=0.333$ with u1 = c1*1.0/(1.0+exp(-(a1*n+b1)))+d1; u2 =
c2*1.0/(1.0+exp(-(a2*n+b2)))+d2; u3 = c3*1.0/(1.0+exp(-(a3*n+b3)))+d3;
z=u1*(1-x).2+u2*2*x.*(1-x)+u3*x.2;

In Figure 18 we present the training over 50 samples, and in Figure 19 we
present the convergence of the training.

We select $n=0.333$ and we compute

$\displaystyle
z(x)=ANN_{1}(n)*B_{1,2}(x)+ANN_{2}(n)*B_{2,2}(x)+ANN_{3}(n)*B_{3,2}(x)=$ (25)
$\displaystyle\left(\frac{c_{1}}{1+exp(-a_{1}*n-b_{1})}+d_{1}\right)*(1-x)^{2}+\left(\frac{c_{2}}{1+exp(-a_{2}*x-b_{2})}+d_{2}\right)*2x(1-x)^{2}+$
(26)
$\displaystyle+\left(\frac{c_{3}}{1+exp(-a_{3}*x-b_{3})}+d_{3}\right)*x^{2}$
(27)

we compare with $sin(0.333\pi x)$ and $sin(0.777\pi x)$ in Figure 20.

##  6 Appendix B: One-dimensional example of neural network learning IGA
solution

###  6.1 One dimensional heat-transfer problem

We focus again on the heat-transfer problem

$-u^{\prime\prime}(x)=f(x)\quad x\in(0,0.5)$ (28)

defined over $x\in(0,0.5)$, with boundary conditions $u(0)=0$ and
$u^{\prime}(0.5)=g(x)$, with $g(x)=n\pi cos(n\pi x)$ and
$f(x)=n^{2}\pi^{2}sin(n\pi x)$. We transform this problem into the weak form
and we discretize with B-spline basis functions

$\displaystyle\begin{bmatrix}\int_{0,1/2}B_{1,2}^{\prime}(x)B_{1,2}^{\prime}(x)dx&\int_{0,1/2}B_{1,2}^{\prime}(x)B_{2,2}^{\prime}(x)dx&\int_{0,1/2}B_{1,2}^{\prime}(x)B_{3,2}^{\prime}(x)dx\\\
\int_{0,1/2}B_{2,2}^{\prime}(x)B_{1,2}^{\prime}(x)dx&\int_{0,1/2}B_{2,2}^{\prime}(x)B_{2,2}^{\prime}(x)dx&\int_{0,1/2}B_{2,2}^{\prime}(x)B_{3,2}^{\prime}(x)dx\\\
\int_{0,1/2}B_{3,2}^{\prime}(x)B_{1,2}^{\prime}(x)dx&\int_{0,1/2}B_{3,2}^{\prime}(x)B_{2,2}^{\prime}(x)dx&\int_{0,1/2}B_{3,2}^{\prime}(x)B_{3,2}^{\prime}(x)dx\\\
\end{bmatrix}\begin{bmatrix}u_{1}\\\ u_{2}\\\ u_{3}\end{bmatrix}=$
$\displaystyle\begin{bmatrix}\int_{0,1/2}B_{1,2}(x)f_{n}(x)dx\\\
\int_{0,1/2}B_{2,2}(x)f_{n}(x)dx\\\ \int_{0,1/2}B_{3,2}(x)f_{n}(x)dx+n\pi
cos(n\pi 0.5)\end{bmatrix}$ (29)

###  6.2 Artificial neural network approximating solution

Let us introduce the artificial neural network

$ANN(n,x)=y$ (30)

where $n$ is the index of the $f_{n}$ function, and $x$ is the argument.

$ANN(n,x)=c\sigma\left(\begin{bmatrix}a_{1}a_{2}\end{bmatrix}\begin{bmatrix}n\\\
x\end{bmatrix}+b\right)+d$ (31)

where the activation function

$\sigma(x)=\frac{1}{1+e^{-x}}$ (32)

###  6.3 Training

The goal of the training is to find values of the weights $a_{1},a_{2},b,c,d$

We prepare a set of samples

  * 1.

We randomly select $n\in(0,1)$ and $x\in(0,1)$

  * 2.

We solve the IGA problem (29) to obtain $\left(u_{1},u_{2},u_{3}\right)$

  * 3.

Input data $(n,x)$, output data
$y(n,x)=u_{1}B_{1,2}+u_{2}B_{2,2}+u_{3}B_{3,2}$

We define the error function

$\displaystyle
e(n,x)=0.5\left(ANN(n,x)-y(n,x)\right)^{2}=0.5\left(c\sigma\left(a_{1}n+a_{2}x+b\right)+d-y(n,x)\right)^{2}=$
$\displaystyle
0.5\left(\left(\frac{c}{1+exp(-a_{1}n-a_{2}x-b)}+d\right)-y(n,x)\right)^{2}$
(33)

Now, we compute the derivatives

$\displaystyle\frac{\partial e(n,x)}{\partial
a_{1}}=\frac{cnexp\left(-a_{1}n-a_{2}x-b\right)\left(\frac{c}{exp\left(-a_{1}n-a_{2}x-b\right)+1}+d-y\right)}{\left(exp\left(-a_{1}n-a_{2}x-b\right)+1\right)^{2}}$
(34) $\displaystyle\frac{\partial e(n,x)}{\partial
a_{2}}=\frac{cxexp\left(-a_{1}n-a_{2}x-b\right)\left(\frac{c}{exp\left(-a_{1}n-a_{2}x-b\right)+1}+d-y\right)}{\left(exp\left(-a_{1}n-a_{2}x-b\right)+1\right)^{2}}$
(35) $\displaystyle\frac{\partial e(n,x)}{\partial
b}=\frac{cexp\left(-a_{1}n-a_{2}x-b\right)\left(\frac{c}{exp\left(-a_{1}n-a_{2}x-b\right)+1}+d-y\right)}{\left(exp\left(-a_{1}n-a_{2}x-b\right)+1\right)^{2}}$
(36) $\displaystyle\frac{\partial e(n,x)}{\partial
c}=\frac{\frac{c}{exp\left(-a_{1}n-a_{2}x-b\right)+1}+d-y}{exp\left(-a_{1}n-a_{2}x-b\right)+1}$
(37) $\displaystyle\frac{\partial e(n,x)}{\partial
d}=\frac{c}{exp\left(-a_{1}n-a_{2}x-b\right)+1}+d-y$ (38)

We loop through the data set $\\{(n,x),y\\}_{(n,x)\in A}$ where $A$ is the set
of selected points from $(0,0.5)\times(0,1)$, and we train $ANN$

  1. 1.

Select $((n,x),y)$

  2. 2.

Compute
$y=ANN(n,x)=c\sigma\left(\begin{bmatrix}a_{1}a_{2}\end{bmatrix}\begin{bmatrix}n\\\
x\end{bmatrix}+b\right)+d$

  3. 3.

Compute $e(n,x)$

  4. 4.

Compute $\frac{\partial e(n,x)}{\partial a_{1}},\frac{\partial
e(n,x)}{\partial a_{2}},\frac{\partial e(n,x)}{\partial b},\frac{\partial
e(n,x)}{\partial c},\frac{\partial e(n,x)}{\partial d}$

  5. 5.

Correct

$\displaystyle a_{1}=a_{1}-\eta*\frac{\partial e(n,x)}{\partial a_{1}}$ (39)
$\displaystyle a_{2}=a_{2}-\eta*\frac{\partial e(n,x)}{\partial a_{2}}$ (40)
$\displaystyle b=b-\eta*\frac{\partial e(n,x)}{\partial b}$ (41)
$\displaystyle c=c-\eta*\frac{\partial e(n,x)}{\partial c}$ (42)
$\displaystyle d=d-\eta*\frac{\partial e(n,x)}{\partial d}$ (43)

where $\eta\in(0,1)$.

Figure 21: Training for the simple artificial neuRal network (30) starting
from $a_{1}=a_{2}=b=c=d=1.0$, for $\eta=0.1$. Figure 22: Convergence of error
for training of ANN

###  6.4 MATLAB implementation

% Creation of dataset
A = [1/5 1/10 1/30; 1/10 2/15 1/10; 1/30 1/10 1/5];
i=1;
for n=0.01:0.001:0.5
for x=0.01:0.001:0.5
rhs= [ (pi*pi*n*n+2*cos(pi*n)-2)/(pi*pi*pi*n*n*n);
(-2*pi*n*sin(pi*n)-4*cos(pi*n)+4)/(pi*pi*pi*n*n*n);
((2-pi*pi*n*n)*cos(pi*n)+2*pi*n*sin(pi*n)-2)/(pi*pi*pi*n*n*n) ];
u=A $\backslash$ rhs;
y=u(1)*(1-x).2+u(2)*2*x.*(1-x)+u(3)*x.2;
dataset_in_n(i)=n;
dataset_in_x(i)=x;
dataset_y(i)=y;
i=i+1;
endfor
endfor
ndataset=i-1;

% Training
a1=1.0; a2=1.0; b=1.0; c=1.0; d=1.0;
eta=0.1;
r = 0 + (1-0).*rand(ndataset,1);
r=r.*ndataset;
for j=1:ndataset
i=floor(r(j));
eval = c*1.0/(1.0+exp(-(a1*dataset_in_n(i)+a2*dataset_in_x(i)+b)))+d;
error = 0.5*(eval-dataset_y(i))2;;
derrorda1 = ( c*dataset_in_n(i)*exp(-a1*dataset_in_n(i)-a2*dataset_in_x(i)-b)*
(c / (exp(-a1*dataset_in_n(i)-a2*dataset_in_x(i)-b)+1)+d-dataset_y(i))
) / power((exp(-a1*dataset_in_n(i)-a2*dataset_in_x(i)-b)+1),2);
a1=a1-eta* derrorda1;
derrorda2 = ( c*dataset_in_x(i)*exp(-a1*dataset_in_n(i)-a2*dataset_in_x(i)-b)*
(c / (exp(-a1*dataset_in_n(i)-a2*dataset_in_x(i)-b)+1)+d-dataset_y(i))
) / power((exp(-a1*dataset_in_n(i)-a2*dataset_in_x(i)-b)+1),2);
a2=a2-eta* derrorda2;
derrordb = ( c*exp(-a1*dataset_in_n(i)-a2*dataset_in_x(i)-b)*
(c / (exp(-a1*dataset_in_n(i)-a2*dataset_in_x(i)-b)+1)+d-dataset_y(i))
) / power((exp(-a1*dataset_in_n(i)-a2*dataset_in_x(i)-b)+1),2);
b=b-eta* derrordb;
derrordc = ( c /
(exp(-a1*dataset_in_n(i)-a2*dataset_in_x(i)-b)+1)+d-dataset_y(i)
) / (exp(-a1*dataset_in_n(i)-a2*dataset_in_x(i)-b)+1);
c=c-eta* derrordc;
derrordd = c /
(exp(-a1*dataset_in_n(i)-a2*dataset_in_x(i)-b)+1)+d-dataset_y(i);
d=d-eta* derrordd;

% evaluation of ANN approximation of sin(n*pi*x) for n=0.333
n=0.333;
x=0:0.01:0.5;
y=sin(n*pi.*x);
eval = c*1.0./(1.0+exp(-(a1*n+a2.*x+b)))+d;
plot(x,y,x,eval);

###  6.5 Verification

In Figure 21 we present the training over 250,000 samples, and in Figure 22 we
present the convergence of the training.

We select $n=0.333$ and we compute

$\displaystyle y(x)=ANN(n,x)=\frac{c}{1+exp(-a_{1}*n-a_{2}*x-b)}+d$ (44)

we compare with $sin(0.333\pi x)$ in Figure 23.

Figure 23: Verification of the neural network approximation of solution for
$n=0.333$ with $y(x)=ANN(n,x)=\frac{c}{1+exp(-a_{1}*n-a_{2}*x-b)}+d$

##  7 Appendix C: One-dimensional example of Physics Informed Neural Network

###  7.1 One dimensional heat-transfer problem

We focus again on the heat-transfer problem

$u^{\prime\prime}(x)+f(x)=0\quad x\in(0,0.5)$ (45)

defined over $x\in(0,0.5)$, with boundary conditions $u(0)=0$ and
$u^{\prime}(0.5)=g(0.5)$, with $g(x)=n\pi cos(n\pi x)$ and
$f(x)=n^{2}\pi^{2}sin(n\pi x)$.

###  7.2 Physics informed neural network

We define the neural network

$PINN(x)=u$ (46)

where

$PINN(x)=c\sigma\left(ax+b\right)+d=\frac{c}{1+exp(-ax-b)}+d$ (47)

We compute the derivatives

$PINN_{x}(x)=\frac{a*c*exp(-ax-b)}{\left(exp(-ax-b)+1\right)^{2}}$ (48)

and

$PINN_{xx}(x)=c\left(\frac{2a^{2}exp(-2ax-2b)}{\left(exp(-ax-b)+1\right)^{3}}-\frac{a^{2}exp(-ax-b)}{\left(exp(-ax-b)+1\right)^{2}}\right)$
(49)

###  7.3 Training

The goal of the training is to find values of the weights $a,b,c,d$

We prepare a set of samples

  * 1.

We randomly select $x\in(0,0.5)$

  * 2.

Input data $x$, output data $u=PINN(x)$

We define

$F(x)=PINN_{xx}(x)+n^{2}\pi^{2}sin(n\pi x)$ (50)

We define the error of approximation of PDE

$\displaystyle
error1(x)=0.5*F(x)^{2}=0.5*\left(PINN_{xx}(x)+n^{2}\pi^{2}sin(n\pi
x)\right)^{2}=$ $\displaystyle
0.5*\left(c\left(\frac{2a^{2}exp(-2ax-2b)}{\left(exp(-ax-b)+1\right)^{3}}-\frac{a^{2}exp(-ax-b)}{\left(exp(-ax-b)+1\right)^{2}}\right)+n^{2}\pi^{2}sin(n\pi
x)\right)^{2}$ (51)

as well as the error of approximation of the boundary condition at $x=0$

$\displaystyle
error2(0)=0.5*\left(PINN(0)-0\right)^{2}=0.5*\left(\frac{c}{1+exp(-b)}+d-0\right)^{2}$
(52)

as well as the error of approximation of the boundary condition at $x=0.5$

$\displaystyle error3(0.5)=0.5*\left(PINN_{x}(0.5)-g(0.5)\right)^{2}=$
$\displaystyle
0.5*\left(\frac{a*c*exp(-a0.5-b)}{\left(exp(-a0.5-b)+1\right)^{2}}-n\pi
cos(n\pi 0.5)\right)^{2}$ (53)

  1. 1.

Select $x$

  2. 2.

Compute $u=PINN(x)=c\sigma\left(ax+b\right)+d=\frac{c}{1+exp(-ax-b)}+d$

  3. 3.

Compute $error1(x)$, $error2(0)$, $error3(0.5)$

  4. 4.

Compute $\frac{\partial error1(x)}{\partial a},\frac{\partial
error1(x)}{\partial b},\frac{\partial error1(x)}{\partial c},\frac{\partial
error1(x)}{\partial d}$

  5. 5.

Compute $\frac{\partial error2(0)}{\partial a},\frac{\partial
error2(0)}{\partial b},\frac{\partial error2(0)}{\partial c},\frac{\partial
error2(0)}{\partial d}$

  6. 6.

Compute $\frac{\partial error3(0.5)}{\partial a},\frac{\partial
error3(0.5)}{\partial b},\frac{\partial error3(0.5)}{\partial
c},\frac{\partial error3(0.5)}{\partial d}$

  7. 7.

Correct

$\displaystyle a=a-\eta*\frac{\partial e(x)}{\partial a}$ (54) $\displaystyle
b=b-\eta*\frac{\partial e(x)}{\partial b}$ (55) $\displaystyle
c=c-\eta*\frac{\partial e(x)}{\partial c}$ (56) $\displaystyle
d=d-\eta*\frac{\partial e(x)}{\partial d}$ (57)

where $e(x)=error1(x)+error2(x)+error3(x)$

for $\eta\in(0,1)$.

We compute

$\displaystyle\frac{\partial PINN_{xx}(x)}{\partial
a}=c\left(\frac{a^{2}x*exp(-ax-b)}{(exp(-ax-b)+1)^{2}}-\frac{6a^{2}x*exp(-2ax-2b)}{(exp(-ax-b)+1)^{3}}\right.$
$\displaystyle\left.+\frac{6a^{2}x*exp(-3ax-3b)}{(exp(-ax-b)+1)^{4}}-\frac{2aexp(-ax-b)}{(exp(-ax-b)+1)^{2}}+\frac{4a*exp(-2ax-2b)}{(exp(-ax-b)+1)^{3}}\right)$
(58) $\displaystyle\frac{\partial PINN_{xx}(x)}{\partial
b}=c\left(\frac{a^{2}exp(ax+b)\left(-4exp(ax+b)+exp(2ax+2b)+1\right)}{(exp(ax+b)+1)^{4}}\right)$
(59) $\frac{\partial PINN_{xx}(x)}{\partial
c}=\left(\frac{2a^{2}exp(-2ax-2b)}{\left(exp(-ax-b)+1\right)^{3}}-\frac{a^{2}exp(-ax-b)}{\left(exp(-ax-b)+1\right)^{2}}\right)$
(60) $\frac{\partial PINN_{xx}(x)}{\partial d}=0$ (61) $\frac{\partial
PINN(0)}{\partial a}=0$ (62) $\frac{\partial PINN(0)}{\partial
b}=\frac{exp(-b)c}{\left(exp(-b)+1\right)^{2}}$ (63) $\frac{\partial
PINN(0)}{\partial c}=\frac{1}{\left(exp(-b)+1\right)}$ (64) $\frac{\partial
PINN(0)}{\partial d}=1.0$ (65) $\displaystyle\frac{\partial
PINN_{x}(0.5)}{\partial
a}=c*exp(b-a)\frac{\left((1-0.5a)*exp(2a+b)+(0.5a+1)exp(1.5a)\right)}{\left(exp(0.5a+b)+1\right)^{3}}$
(66) $\displaystyle\frac{\partial PINN_{x}(0.5)}{\partial
b}=\frac{ac*exp(b-0.5a)\left(exp(a)-exp(1.5a+b)\right)}{\left(exp(0.5a+b)+1\right)^{3}}$
(67) $\displaystyle\frac{\partial PINN_{x}(0.5)}{\partial
c}=\frac{a*exp(-0.5a-b)}{\left(exp(-0.5a-b)+1\right)^{2}}$ (68)
$\displaystyle\frac{\partial PINN_{x}(0.5)}{\partial d}=0$ (69)

Using the above formulas we have

$\displaystyle\frac{\partial error1(x)}{\partial
a}=\left(c\left(\frac{2a^{2}exp(-2ax-2b)}{\left(exp(-ax-b)+1\right)^{3}}-\frac{a^{2}exp(-ax-b)}{\left(exp(-ax-b)+1\right)^{2}}\right)+n^{2}\pi^{2}sin(n\pi
x)\right)$ $\displaystyle
c\left(\frac{a^{2}x*exp(-ax-b)}{(exp(-ax-b)+1)^{2}}-\frac{6a^{2}x*exp(-2ax-2b)}{(exp(-ax-b)+1)^{3}}\right.$
$\displaystyle\left.+\frac{6a^{2}x*exp(-3ax-3b)}{(exp(-ax-b)+1)^{4}}-\frac{2aexp(-ax-b)}{(exp(-ax-b)+1)^{2}}+\frac{4a*exp(-2ax-2b)}{(exp(-ax-b)+1)^{3}}\right)$
(70) $\displaystyle\frac{\partial error1(x)}{\partial
b}=\left(c\left(\frac{2a^{2}exp(-2ax-2b)}{\left(exp(-ax-b)+1\right)^{3}}-\frac{a^{2}exp(-ax-b)}{\left(exp(-ax-b)+1\right)^{2}}\right)+n^{2}\pi^{2}sin(n\pi
x)\right)$ $\displaystyle
c\left(\frac{a^{2}exp(ax+b)\left(-4exp(ax+b)+exp(2ax+2b)+1\right)}{(exp(ax+b)+1)^{4}}\right)$
(71) $\displaystyle\frac{\partial error1(x)}{\partial
c}=\left(c\left(\frac{2a^{2}exp(-2ax-2b)}{\left(exp(-ax-b)+1\right)^{3}}-\frac{a^{2}exp(-ax-b)}{\left(exp(-ax-b)+1\right)^{2}}\right)+n^{2}\pi^{2}sin(n\pi
x)\right)$
$\displaystyle\left(\frac{2a^{2}exp(-2ax-2b)}{\left(exp(-ax-b)+1\right)^{3}}-\frac{a^{2}exp(-ax-b)}{\left(exp(-ax-b)+1\right)^{2}}\right)$
(72) $\displaystyle\frac{\partial error1(x)}{\partial d}=0$ (73)

$\displaystyle\frac{\partial error2(x)}{\partial a}=0$ (74)
$\displaystyle\frac{\partial error2(x)}{\partial
b}=\left(\frac{c}{1+exp(-b)}+d-0\right)\frac{exp(-b)c}{\left(exp(-b)+1\right)^{2}}$
(75) $\displaystyle\frac{\partial error2(x)}{\partial
c}=\left(\frac{c}{1+exp(-b)}+d-0\right)$
$\displaystyle\left(\frac{2a^{2}exp(-2ax-2b)}{\left(exp(-ax-b)+1\right)^{3}}-\frac{a^{2}exp(-ax-b)}{\left(exp(-ax-b)+1\right)^{2}}\right)$
(76) $\displaystyle\frac{\partial error2(x)}{\partial
d}=\left(\frac{c}{1+exp(-b)}+d-0\right)$ (77)

$\displaystyle\frac{\partial error3(x)}{\partial
a}=\left(\frac{a*c*exp(-a0.5-b)}{\left(exp(-a0.5-b)+1\right)^{2}}-n\pi
cos(n\pi 0.5)\right)$ $\displaystyle
c*exp(b-a)\frac{\left((1-0.5a)*exp(2a+b)+(0.5a+1)exp(1.5a)\right)}{\left(exp(0.5a+b)+1\right)^{3}}$
(78) $\displaystyle\frac{\partial error3(x)}{\partial
b}=\left(\frac{a*c*exp(-a0.5-b)}{\left(exp(-a0.5-b)+1\right)^{2}}-n\pi
cos(n\pi 0.5)\right)$
$\displaystyle\frac{ac*exp(b-0.5a)\left(exp(a)-exp(1.5a+b)\right)}{\left(exp(0.5a+b)+1\right)^{3}}$
(79) $\displaystyle\frac{\partial error3(x)}{\partial
c}=\left(\frac{a*c*exp(-a0.5-b)}{\left(exp(-a0.5-b)+1\right)^{2}}-n\pi
cos(n\pi 0.5)\right)$
$\displaystyle\frac{a*exp(-0.5a-b)}{\left(exp(-0.5a-b)+1\right)^{2}}$ (80)
$\displaystyle\frac{\partial error3(x)}{\partial d}=0$ (81)

Figure 24: Training of the simple PINN (47) starting from $a=b=d=1.0$, and
$c=3.0$, for $\eta=0.1$. Figure 25: Convergence of error for training of PINN

###  7.4 MATLAB implementation

% Creation of dataset
i=1;
n=0.333;
for x=0.01:0.01:0.5
y=sin(n*pi*x);
dataset_in_x(i)=x;
dataset_y(i)=y;
i=i+1;
endfor
ndataset=i-1;

% Training
a1=1.0; a2=1.0; b=1.0; c=3.0; d=1.0;
eta=0.1;
r = 0 + (1-0).*rand(ndataset,1);
r=r.*ndataset;
for j=1:ndataset
i=floor(r(j));
eval = c*1.0/(1.0+exp(-(a*dataset_in_x(i)+b)))+d;
error = 0.5*(eval-dataset_y(i))2;
% Training of the PDE
x = dataset_in_x(i);
n = dataset_in_n(i);
Fx = c*(2*a*a*exp(-2*a*x-2*b) / power((exp(-a*x-b)+1),3)
\- a*a*exp(-a*x-b) / power((exp(-a*x-b)+1),2))+n*n*pi*pi*sin(n*x);
derror1da = Fx*c *( (a*a*x*exp(-a*x-b))/power((exp(-a*x-b)+1),2)
-(6*a*a*x*exp(-2*a*x-2*b))/power((exp(-a*x-b)+1),3)
+(6*a*a*x*exp(-3*a*x-3*b))/power((exp(-a*x-b)+1),4)
-(2*a*exp(-a*x-b))/power((exp(-a*x-b)+1),2)
+(4*a*exp(-2*a*x-2*b))/power((exp(-a*x-b)+1),3) );
a=a-eta* derror1da;
derror1db = Fx*c*( (a*a*exp(a*x+b))*(-4*exp(a*x+b)
+exp(2*a*x+2*b)+1)/power((exp(a*x+b)+1),4));
b=b-eta* derror1db;
derror1dc = Fx*( (2*a*a*exp(-2*a*x-2*b))/power((exp(-a*x-b)+1),3)
-(a*a*exp(-a*x-b))/power((exp(-a*x-b)+1),2));
c=c-eta* derror1dc;
derror1dd = 0;
d=d-eta* derror1dd;
% Training of the boundary condition at x=0
x=0;
derror2da = 0;
a=a-eta* derror2da;
derror2db = (c / (1+exp(-b))+d)* (exp(-b)*c)/power((exp(-b)+1),2);
b=b-eta* derror2db;
derror2dc = (c/(1+exp(-b))+d)*( (2*a*a*exp(-2*a*x-2*b))/power(exp(-a*x-b)+1,3)
\- (a*a*exp(-a*x-b))/power(exp(-a*x-b)+1,2) );
c=c-eta* derror2dc;
derror2dd = c/(1+exp(-b))+d;
d=d-eta* derror2dd;
% Training of the boundary condition at x=0.5
x=0.5;
derror3da = (a*c*exp(-a*0.5-b)/power((exp(-a*0.5-b)+1),2)
-n*pi*cos(n*pi*0.5)) *c*exp(b-a)*((1-0.5*a)*exp(2*a+b)
+(0.5*a+1)*exp(1.5*a))/power((exp(0.5*a+b)+1),3);
a=a-eta* derror3da;
derror3db = (a*c*exp(-a*0.5-b)/power((exp(-a*0.5-b)+1),2)
-n*pi*cos(n*pi*0.5))* a*c*exp(b-0.5*a)*(exp(a)+exp(1.5*a+b))/power((exp(0.5*a+b)+1),3);
b=b-eta* derror3db;
derror3dc = (a*c*exp(-a*0.5-b)/power((exp(-a*0.5-b)+1),2)
-n*pi*cos(n*pi*0.5))* a*exp(-b-0.5*a)/power((exp(-0.5*a-b)+1),2);
c=c-eta* derror3dc;
derror3dd = 0;
d=d-eta* derror3dd;
endfor

% evaluation of PINN approximation of sin(0.333*pi*x)
n=0.333;
x=0:0.01:0.5;
y=sin(n*pi.*x);
eval = c*1.0./(1.0+exp(-(a.*x+b)))+d;
plot(x,y,x,eval);

###  7.5 Verification

In Figure 24 we present the training over 50 samples, and in Figure 25 we
present the convergence of the training.

The PINN has been trained for $n=0.333$ so we compute

$\displaystyle y(x)=ANN(n,x)=\frac{c}{1+exp(a*x-b)}+d$ (82)

we compare with $sin(0.333\pi x)$ in Figure 26.

Figure 26: Verification of the PINN trained for $n=0.333$ with
$y(x)=PINN(n,x)=\frac{c}{1+exp(a*x-b)}+d$

## References

  * Brevis et al. [2021] Ignacio Brevis, Ignacio Muga, and Kristoffer G. van der Zee.  A machine-learning minimal-residual (ml-mres) framework for goal-oriented finite element discretizations.  Computers & Mathematics with Applications, 95:186–199, 2021.  Recent Advances in Least-Squares and Discontinuous Petrov–Galerkin Finite Element Methods.
  * Chen et al. [2018] Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud.  Neural ordinary differential equations.  Advances in Neural Information Processing Systems 31, 2018.
  * Collier et al. [2012] Nathan Collier, David Pardo, Lisandro Dalcin, Maciej Paszynski, and V.M. Calo.  The cost of continuity: A study of the performance of isogeometric finite elements using direct solvers.  Computer Methods in Applied Mechanics and Engineering, 213-216:353–361, 2012.
  * Dalcin et al. [2016] L. Dalcin, N. Collier, P. Vignal, A.M.A. Côrtes, and V.M. Calo.  Petiga: A framework for high-performance isogeometric analysis.  Computer Methods in Applied Mechanics and Engineering, 308:151 – 181, 2016.
  * Haghighat et al. [2020] Ehsan Haghighat, Maziar Raissi, Adrian Mourec, Hector Gomez, and Ruben Juanes.  A deep learning framework for solution and discovery in solid mechanics.  arXiv:2003.02751, 2020.
  * Hornik et al. [1989] Kurt Hornik, Maxwell Stinchcombe, and Halbert White.  Multilayer feedforward networks are universal approximators.  Neural Networks, 1989.
  * HP [2021] Deep learning driven self-adaptive hp finite element method.  Lecture Notes in Computer Science, (12742):114–121, 2021.
  * Hughes et al. [2005] Thomas J.R. Hughes, J.A. Cottrell, and Yuri Bazilevs.  Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement.  Computer Methods in Applied Mechanics and Engineering, 2005.
  * Lagaris et al. [1998] I.E. Lagaris, A. Likas, and D.I. Fotiadis.  Artificial neural networks for solving ordinary and partial differential equations.  IEEE Transactions on Neural Networks, 9(5):987–1000, 1998.
  * Michoski et al. [2020] Craig Michoski, Miloš Milosavljević, Todd Oliver, and David R. Hatch.  Solving differential equations using deep neural networks.  Neurocomputing, 399:193–212, 2020.
  * Misyris et al. [2019] George S. Misyris, Andreas Venzke, and Spyros Chatzivasileiadis.  Physics-informed neural networks for power systems.  preprint at https://www.researchgate.net/publication/337184659_Physics-Informed_Neural_Networks_for_Power_Systems, 2019\.
  * Nguyen et al. [2015] Vinh Phu Nguyen, Cosmin Anitescu, Stéphane P.A. Bordas, and Timon Rabczuk.  Isogeometric analysis: An overview and computer implementation aspects.  Mathematics and Computers in Simulation, 117:89 – 116, 2015.
  * Pang et al. [2020] Guofei Pang, Marta D’Elia, Michael Parks, and George Karniadakis.  npinns: Nonlocal physics-informed neural networks for a parametrized nonlocal universal laplacian operator. algorithms and applications.  Journal of Computational Physics, 422, 2020.
  * Paszyński [2020] Maciej Paszyński.  Classical and isogeometric finite element method.  AGH University of Science and Technology in Krakow, 2020.  https://epodreczniki.open.agh.edu.pl/handbook/1088/module/1173/reader.
  * Raissi et al. [2017] M. Raissi, P. Perdikaris, and G. Karniadakis.  Physics informed deep learning (part i): Data-driven solutions of nonlinear partial differential equations.  ArXiv, abs/1711.10561, 2017.
  * Raissi et al. [2019] Maziar Raissi, Paris Perdikaris, and George E Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics, 378:686–707, 2019.
  * Wang et al. [2020] Sifan Wang, Xinling Yu, and Paris Perdikaris.  When and why pinns fail to train: A neural tangent kernel perspective.  preprint at https://arxiv.org/abs/2007.14527, 07 2020.


# Extended Anderson Criticality in Heavy-Tailed Neural Networks

Asem Wardak  Pulin Gong<EMAIL_ADDRESS>School of Physics,
University of Sydney, New South Wales 2006, Australia

###### Abstract

We investigate the emergence of complex dynamics in networks with heavy-tailed
connectivity by developing a non-Hermitian random matrix theory. We uncover
the existence of an extended critical regime of spatially multifractal
fluctuations between the quiescent and active phases. This multifractal
critical phase combines features of localization and delocalization and
differs from the edge of chaos in classical networks by the appearance of
universal hallmarks of Anderson criticality over an extended region in phase
space. We show that the rich nonlinear response properties of the extended
critical regime can account for a variety of neural dynamics such as the
diversity of timescales, providing a computational advantage for persistent
classification in a reservoir setting.

## Introduction.—

In a diverse range of physical, biological, financial and ecological systems,
complex dynamics fluctuating across multiple scales emerge from a large number
of interacting, nonlinear units with heterogeneous properties. Understanding
the organising principles and behavior of such complex dynamics is a
longstanding topic of interest across these diverse fields [1, *Stanley1988,
*Ghashghaie1996]. In neuroscience and machine learning, neural networks with
many interacting neurons likewise exhibit complex dynamics with large
fluctuations which are critical for their information processing abilities on
real-world inputs [4, *NIPS2016_14851003]. However, the network mechanisms and
fundamental computational capabilities of complex neural dynamics remain
elusive.

The classical formulation of complex dynamics in systems with many interacting
elements is based on randomly connected networks of coupled, nonlinear units
with homogeneous connectivity [6, *PhysRevX.5.041030]. By employing
mathematical approaches such as mean-field theory and random matrix theory,
such networks have robustly predicted a phase of chaotic activity with global,
homogeneous (i.e. delocalized) fluctuations existing adjacent to an ordered,
silent regime, enabling the analysis of a wide range of systems with a
characteristic scale [8, *PhysRevE.90.062710, *PhysRevE.82.011903,
*PhysRevLett.110.118101]. The edge of the ordered and chaotic phases gives
rise to a range of critical phenomena which is thought to be necessary for
these systems to perform useful computations [12]. However, growing evidence
has shown that coupling heterogeneity is widespread in complex systems such as
biological [13, *10.7554/eLife.57443, *Buzsaki2014] and artificial neural
networks [16, 17], underscoring the need to understand the fundamental
dynamical and computational mechanisms of such heterogeneity.

Here we study the dynamics of random neural networks with heterogeneous,
heavy-tailed connectivity. After describing the fixed points of the system
using a Lévy mean-field approach, we develop a novel non-Hermitian random
matrix theory for column-structured heavy-tailed matrices to analyze the
statistical fluctuations of random neural networks around the fixed point.
This theory reveals a new regime with correlated multifractal modes which are
neither localized nor delocalized, but have aspects of both (Fig. 1).
Multifractality is characterized by the appearance of differing, nontrivial
structures appearing simultaneously over a wide variety of scales (first
proposed to explain fluctuations in turbulence by Parisi [18]), and is a
hallmark of Anderson transitions (criticality) [19]. Anderson transitions were
first described in the context of disordered electronic systems with localized
and metallic (i.e. delocalized) phases [20], and have since been analysed in a
broad sense in a wide range of systems including treelike Bethe lattices [21,
*Parisi_2019] and those exhibiting conventional second-order phase transitions
[19].

We illustrate that the heavy-tailed heterogeneity in connectivity enables
Anderson criticality to emerge in a broad parameter regime. The correlated
multifractal modes characteristic of this extended critical regime are able to
explain a range of realistic neural dynamics, including correlated
fluctuations with low-dimensional features [23, *HUANG2019337], long-range
correlations [25] and a diversity of timescales [26]. Importantly, these
correlated multifractal modes provide a profound computational advantage in
the setting of real-time reservoir computing by allowing for a persistent form
of dimensionality expansion, which is not possible in classical homogeneous
systems [27].

(a)

(b)

Figure 1: Schematic of the localization of activity fluctuations arising in
neural networks. (a) In the homogeneous neural network (blue), activity
fluctuations are delocalized and spread evenly over the spatial extent of the
system ($D_{q}=1$). (b) In the heavy-tailed neural network (red) with
heterogeneous weights (line thickness), activity fluctuations are multifractal
with a nontrivial distribution of localization over system sites ($0\leq
D_{q}\leq 1$) visible over multiple scales (inset).

## Network model and fixed points.—

We begin by extending the seminal random neural network model with interacting
nonlinear units analysed by Sompolinsky et al. and others [7, 11, 8, 9, 10],
which has the dynamics

$\displaystyle
h_{i}^{\prime}(t)=-h_{i}(t)+g\sum_{j=1}^{N}J_{ij}\phi(h_{j}(t))~{},$ (1)

where $h_{i}(t)$ is the input of the $i$-th neuron at time $t$, $J_{ij}$ is
the strength of the connection from neuron $j$ to neuron $i$, and $\phi$ is a
scalar nonlinearity that determines the neural firing rate given the input.
The theoretical results in this Letter apply for a wide range of $\phi$
explained below, while the figures use $\phi=\tanh$ for comparability with
previous models [6]. The coefficient $g$ is the gain parameter of the synaptic
input. Our aim is to investigate networks with heavy-tailed heterogeneity
which do not fall within the purview of the Gaussian large-size limit. Finite-
size heterogeneous networks can be investigated to some degree using
perturbative finite-order corrections to the homogeneous Gaussian limit [28,
*PhysRevLett.127.158302]. We instead take an approach which is exact in the
large network limit by regarding each $J_{ij}$ as an independent random
variable whose second moment is not finite, so that its probability density
has a power-law asymptotic tail,

$\displaystyle
p_{J_{ij}}(x)\stackrel{{\scriptstyle|x|\to\infty}}{{\sim}}\frac{C_{\alpha}}{2N|x|^{1+\alpha}}~{},$
(2)

where $C_{\alpha}:=\Gamma(1+\alpha)\sin(\pi\alpha/2)/\pi$ and $1<\alpha<2$.
Such heavy-tailed connectivity has been observed in Drosophila [13, 14], in
successfully trained artificial neural networks [17] and in spin-glass systems
with strong disorder [30]; theoretically elucidating the dynamical impact of
such heterogeneity has drawn increasing attention [31,
*PhysRevResearch.3.013083]. Our approach includes asymmetric Lévy matrices
$J_{ij}\sim L(\alpha,0,0,1/2N)$ with independent $\alpha$-stable [33] elements
of scale parameter $1/2N$. To prevent self-interference we set $J_{ii}=0$,
which does not change the overall statistical properties of the connectivity
matrix as $N\to\infty$. In the limit $\alpha\to 2$, we recover the classical
rate-based random neural network [6] with independent elements of variance
$1/N$.

Applying the generalized central limit theorem to Eq. (1) in the large network
limit yields a Lévy mean-field theory for the dynamics of the system at long
times, whereby each neuron behaves as

$\displaystyle h_{i}^{\prime}(t)=-h_{i}(t)+\eta_{i}(t)~{},$ (3)

where $\eta_{i}(t)\sim
L(\alpha,0,0,g^{\alpha}\langle|\phi(h_{j}(t))|^{\alpha}\rangle_{j}/2)$ is a
time-dependent $\alpha$-stable field generated by randomly weighted inputs
from the other neurons. Setting the left hand side of Eq. (3) to zero gives
the fixed points of the network, whose components are distributed as an
$\alpha$-stable random variable for all values of the gain $g$. As we discover
in our Jacobian analysis below, the zero fixed point is unstable so that
network activity is nonzero and has a macroscopic number of nonzero fixed
points at any gain $g$ for $\alpha<2$, in contrast to homogeneous networks
($\alpha=2$) for which the network exhibits a single stable fixed point when
$g<1$. Moreover, since the synaptic input $\eta_{i}(t)$ has power-law
fluctuations over neurons $i$ for heavy-tailed networks with $\alpha<2$, the
synaptic input autocorrelation vector
$\langle\eta_{i}(t)\eta_{i}(t+\tau)\rangle_{t}$ has the asymptotic tail of an
$\alpha/2$-stable distribution over neurons $i$, making the averaged synaptic
input autocorrelation $\langle\eta_{i}(t)\eta_{i}(t+\tau)\rangle_{i,t}$
infinite in general along with the local-field autocorrelation $\langle
h_{i}(t)h_{i}(t+\tau)\rangle_{i,t}$. This non-self-averaging property of
heavy-tailed network input differs greatly from classical homogeneous networks
where a Gaussian dynamical mean-field theory predicts that all neurons of a
sufficiently large population have the same autocorrelation function, and thus
the same relaxation timescale [7]. On the other hand, the heavy-tailed
distribution of local-field autocorrelations over neurons in heterogeneous
neural networks allows for a wide reservoir of timescales across neurons of
the same population even in the large network limit.

## Network stability and heavy-tailed random matrix theory.—

To determine the local stability of the network around the fixed points
predicted by the above mean-field approach, we analyze the Jacobian matrix
$-I+gJ\operatorname{diag}_{j}(\phi^{\prime}(h_{j}))$ obtained from Eq. (1)
where $\operatorname{diag}_{j}\chi_{j}$ denotes the diagonal matrix with
entries $\chi_{j}$. Shifting this Jacobian matrix yields a stability matrix
obtained by scaling the columns of $gJ$ by $\phi^{\prime}(h_{j})$, which has
the form of a column-structured non-Hermitian random matrix 111For clarity, we
refer to $gJ\operatorname{diag}_{j}(\phi^{\prime}(h_{j}))$ as the Jacobian
matrix in the remainder of the Letter as its eigenvectors remain unchanged and
its eigenvalues are shifted one unit to the right. Since the heavy-tailed
matrix $J$ (and thus $J\operatorname{diag}_{j}\chi_{j}$) has a locally
treelike structure [35], we develop a new cavity approach for column-
structured non-Hermitian heavy-tailed random matrices to obtain the spectral
density and eigenvector localization properties of the Jacobian 222See
Supplemental Material.. In recent years, cavity approaches have been used in
asymmetrically disordered contexts involving non-Hermitian random ensembles by
mapping the problem back to a symmetric, Hermitian system of twice the
dimensionality [37]. Our cavity approach to the column-structured matrix
$J\operatorname{diag}_{j}\chi_{j}$ for any $\chi_{j}$ with
$\langle|\chi_{j}|^{\alpha}\rangle_{j}<\infty$ yields the spectral density

$\displaystyle\rho(z)=\frac{y_{*}^{2}-2|z|^{2}y_{*}\partial_{|z|^{2}}y_{*}}{\pi}\left\langle\frac{|\chi_{i}|^{2}SS^{\prime}}{(|z|^{2}+|\chi_{i}|^{2}y_{*}^{2}SS^{\prime})^{2}}\right\rangle_{i}$
(4)

and inverse participation ratio (see Ref. [36] for mathematical derivations)
describing the spatial localization of the left and right eigenvectors as a
function of eigenvalue modulus $|z|$. $\langle..\rangle_{i}$ denotes averaging
over $i$ and any relevant random variables, while $S,S^{\prime}\sim
L(\alpha/2,1,0,C_{\alpha}/4C_{\alpha/2})$ are independent, skewed
$\alpha/2$-stable random samples, and $y_{*}$ is found by solving the equation

$\displaystyle
1=\left\langle\left(\frac{|\chi_{i}|^{2}S}{|z|^{2}+y_{*}^{2}|\chi_{i}|^{2}SS^{\prime}}\right)^{\alpha/2}\right\rangle_{i}~{}.$
(5)

Our random matrix theory thus unifies both classical results in column-
structured non-Hermitian random matrix theory [8, 38, *PhysRevE.85.066116] for
which $\alpha=2$ and the random variables $S,S^{\prime}$ reduce to the
constant $1$, and unstructured heavy-tailed non-Hermitian random matrices [35]
for which $\chi=1$.

| \begin{overpic}[width=216.81pt]{fig/spectrum_alpha_120.pdf}
\put(0.0,80.0){(a)} \end{overpic}
---
\begin{overpic}[width=216.81pt]{fig/relaxationtime_0.25_alpha_120.pdf}
\put(0.0,80.0){(b)} \end{overpic}
\begin{overpic}[width=325.215pt]{fig/phasediagram.pdf} \put(0.0,103.0){(c)} \end{overpic} |  | \begin{overpic}[width=216.81pt]{fig/IPR.pdf} \put(0.0,80.0){(d)} \end{overpic}
---
\begin{overpic}[width=216.81pt]{fig/D_q.pdf} \put(0.0,80.0){(e)} \end{overpic}
Figure 2: The extended Anderson critical regime in heavy-tailed random neural
networks. (a) Jacobian spectra
$J\operatorname{diag}_{j}(\tanh^{\prime}(h_{j}))$ in the stationary state
(lines) with characteristic spectral radius $r_{0.01}$ (dashed circle) and
their numerical validation (dots). (b) Distribution of the half width at
quarter-maximum of the local-field autocorrelation for $g=0.75,1.5,2.25,3$
(red, blue, yellow, purple). Log-log plot (inset) with comparison power-law
line of exponent $-2.0$. (c) Phase diagram over heavy-tailed index $\alpha$
and gain $g$. The extended Anderson critical regime emerging between the lower
and upper dashed lines predicted by Eq. (6) (labeled by $n$) corresponds with
high relaxation time variance (colored). (d) The logarithm of the IPR of the
heavy-tailed network activity fluctuation vector (red) varies widely between
the asymptotic large-$N$ delocalized ($-\log_{10}N\approx-3$) and localized
($0$) values, while the classical network’s activity fluctuations remain
delocalized (blue). (e) The multifractal dimension $D_{q}$ estimated via the
asymptotic relation $D_{q}\sim(\log_{N}\mathrm{IPR}_{q})/(1-q)$, shaded to one
standard deviation over time. $D_{q}$ estimate for a random delocalized vector
on the $N$-sphere (dashed) for comparison.

Figure 2(a) shows the eigenvalue density $\rho(z)$ of the Jacobian and its
numerical validation as a function of eigenvalue modulus $|z|$. A key
characteristic of heavy-tailed neural networks is the infinite spectral radius
$r_{0}$, given by the point at which $y_{*}=0$, so that the zero fixed point
only occurs for zero gain. However, since the eigenvalue density of the heavy-
tailed Jacobian is exponentially suppressed at large radius [34, 35], a
characteristic spectral radius $r_{p}$ can be defined by the eigenvalue
modulus $|z|$ at which the parameter $y_{*}$ drops to a fraction $p$ of its
value at $z=0$. Meanwhile, the spectral radius $r_{0}$ of the Jacobian for
classical networks ($\alpha=2$) is finite and is obtained by a spectral
analysis of finite-variance block-structured random matrices [38, 39, 8]. In
this sense, $y_{*}$ behaves as an order parameter indicating the transition
between microscopic and macroscopic numbers of eigenstates relative to system
size $N$ at a given eigenvalue modulus $|z|$. Thus, heavy-tailed networks
exhibit a quasi-ordered regime with a fixed point near zero whose magnitude is
suppressed due to the microscopic number of eigenstates above the Jacobian
stability line $\operatorname{Re}z=1$. A continuous transition parameterized
by $p\ll 1$ is then defined by the value of the gain parameter $g$ at which
the characteristic spectral radius $r_{p}=1$, distinguishing a quiescent phase
from a chaotic regime (Fig. 2(c), lower dashed line). This transition is
consistent with the point at which neural activity predicted by mean-field
theory (Eq. (3)) deviates significantly from zero (Fig. 2(c), red line).

## An extended critical phase with correlated multifractal modes.—

Using our cavity approach, we find that all of the right eigenvectors of the
heavy-tailed network Jacobian around the stationary state are multifractal for
sigmoidal $\phi$, a hallmark of Anderson transitions (see Ref. [36] for
derivation). The activity of heavy-tailed neural networks is thus dominated by
multifractal chaotic fluctuations in contrast to the spatially delocalized
chaos appearing in classical models. This theoretical prediction on heavy-
tailed network dynamics is confirmed by simulations of temporal fluctuations
of homogeneous and heavy-tailed neural activity (Fig. 2(d–e), red) which are
both chaotic at $g=1.75$. This result is remarkable from a physical standpoint
as the activity itself is delocalized due to the bounded activation function
$\phi=\tanh$, and does not visibly differ significantly from classical
networks with Gaussian dynamics (see Ref. [36]).

To investigate the behavior of the system’s multifractal modes over long
timescales, we quantify the extent to which Jacobian eigenvalues $\lambda_{i}$
corresponding to a given eigenvector change their modulus relative to unity
when new samples are chosen from the stationary distribution of neural
activity. We thus consider the Jacobian average

$\displaystyle\langle(|\lambda_{i}|-1)^{n}\rangle_{i}$
$\displaystyle=\int_{\mathbb{C}}(|z|-1)^{n}\rho(z)\,dz$ (6)

which penalises small and rewards large eigenvalues to a degree determined by
$n$, and we find a region for the gain $g$ adjacent to and above the ordered
transition line which is characterised by a greater proportion of eigenvalues
away from zero compared to the ordered transition (Fig. 2(c), upper dashed
lines). This allows us to distinguish between an active chaotic region in
which unstable fluctuations tend to be quickly suppressed in favour of new
fluctuations, and a region of temporally correlated chaotic fluctuations. The
continuous transition between these active and correlated chaotic regimes is
parameterized by the annealing strength $n$ in Eq. (6). The theoretically
predicted correlated region closes into the well-known critical point at the
ordered-chaotic phase transition ($g=1$) for classical rate-based networks
($\alpha=2$), supporting the notion that the extended region of correlated
multifractal modes is a critical regime. This extended critical phase is
characterized by a significantly nonzero stationary state and a macroscopic
proportion of unstable eigenstates relative to system dimensionality $N$,
which is fundamentally different from the classical edge of chaos occurring
around the zero fixed point when a microscopic proportion of eigenstates
crosses the stability line. Consequently, this extended critical phase remains
chaotic rather than existing solely at the edge of a chaotic phase as in
classical networks [11].

In summary, the extended critical phase of temporally correlated, spatially
multifractal fluctuations provides a demonstration of how various aspects of
realistic neural dynamics may be exhibited simultaneously, such as long-range
correlations [25] and low spatial dimensionality relative to system size [23].
This latter property arises from the localization of spatially multifractal
fluctuations onto a small number of sites relative to system size (Fig. 1).
Such behavior, along with the non-self-averaging properties characteristic of
Anderson criticality [40] which we derived for local-field autocorrelation,
suggests that the timescales across neurons in the heavy-tailed neural network
are diverse. To validate this theoretical prediction on the extended critical
region, we compute the relaxation timescales of neural autocorrelations over
random networks across heavy-tailed index $\alpha$ and network gain $g$ (Fig.
2(c)). We find that the extended critical phase is characterized by diverse
timescales (Fig. 2(b)): the relaxation timescale distribution is power-law
with index $-2.0$, consistent with that seen in cortical memory traces [26].
This behavior only exists in the critical regime [36].

Figure 3:  Persistent dimensionality expansion in heavy-tailed random neural
networks. Top, middle panels: The evolution of average signal and noise
distances in a rate-based model with $N=250$ and $g=3$, as measured by the
Hamming distance $H_{12}:=\|\phi(h_{i}^{(1)})-\phi(h_{i}^{(2)})\|_{i}^{2}$
between states $\phi(h_{i}^{(1)})$ and $\phi(h_{i}^{(2)})$. Bottom panel:
Classification accuracy of heavy-tailed (red) and classical (blue) rate-based
neural networks.

## Persistent reservoir computing.—

To explore the computational implications of the extended critical regime of
correlated multifractal chaos, we consider a reservoir computing task
described in [27] exploiting the chaotic dimensionality expansion of neural
representations $\Delta H_{12}=H_{12}^{(s)}-H_{12}^{(n)}$ between signal
($H_{12}^{(s)}$) and noise ($H_{12}^{(n)}$) distances (see Ref. [36] for setup
details). Heavy-tailed networks make use of the temporal correlations of
multifractal Jacobian eigenvectors above the stability line to enact a
persistent form of real-time computation in a reservoir computing context.
Because the extended critical regime is chaotic in heavy-tailed networks with
a significantly nonzero fixed point, this regime is able to perform
dimensionality expansion on its input in contrast to the classical edge of
chaos which resides around the zero fixed point. At the same time, the
correlated multifractal Jacobian eigenvectors work to hold off the onset of
mixing (Fig. 3, yellow) to perform persistent chaotic dimensionality expansion
on its input, allowing the computed result to remain in the system and thus
the classification performance to stay above the baseline ($0.02$) for a
longer period of time (Fig. 3, bottom). The extended critical regime of
correlated multifractal chaos is thus able to produce efficient neural
representations balancing the dimensional compression of stimuli (Fig. 3 top,
red and blue), which is useful for generalization [41], and the separation of
stimuli, in order to enact a form of persistent real-time computation.

## Discussion.—

Our theory rigorously demonstrates that heterogeneous, heavy-tailed
connectivity can endow neural circuits with Anderson criticality over an
extended parameter region, thus eliminating the fine-tuning needed in
homogeneously connected neural networks. The Anderson criticality is
characterized by correlated, low-dimensional fluctuations [25, 23] and a
diverse reservoir of timescales [26] as observed in biological neural systems;
these observations otherwise remain unexplained in a unified manner in
conventional theories. Moreover, the extended Anderson criticality provides a
unique mechanism for combining robust real-time computation with long-term
memory of the computed output. Both homogeneous and heavy-tailed networks
utilize chaos to enhance the separation of inputs for linear classification,
but after a transient period of high separability, mixing dominates and erases
the computed output from the classical homogeneous system. The recent
observations on the ubiquity of heavy-tailed coupling in pretrained deep
neural networks [17] suggest that our theory would be powerful for revealing
the shared dynamical principles for persistent computation in both biological
and artificial neural networks.

By making a novel link between Anderson criticality and the highly fluctuating
complex dynamics of neural networks, our results suggest that complex systems
operating over multiple scales should display a degree of multifractality at
some level of fluctuations of system activity, even when the activity itself
is bounded due to physical constraints. Multifractal phenomena have indeed
been seen in a wide variety of natural systems such as turbulence [18]. In
statistical and condensed matter physics, the Anderson transition appears with
multifractality around the boundary of two phases characterised by localized
and delocalized eigenvectors [19, 30]. Our heavy-tailed neural network may
thus be viewed as an inversion of these canonical models: a structurally
extended Anderson regime $0<\alpha<2$ is bounded by delocalized Jacobian modes
appearing in classical homogeneous networks with $\alpha=2$, and localized
eigenstates appearing in the limit $\alpha\to 0$ corresponding to sparse
matrices and directed random graphs with a small average degree [42]. Hence
our random matrix theory could be applied to understand how complex dynamics
emerge in physical systems.

## References

  * Dorogovtsev  _et al._ [2008] S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Critical phenomena in complex networks, Rev. Mod. Phys. 80, 1275 (2008).
  * Stanley and Meakin [1988] H. E. Stanley and P. Meakin, Multifractal phenomena in physics and chemistry, Nature 335, 405 (1988).
  * Ghashghaie  _et al._ [1996] S. Ghashghaie, W. Breymann, J. Peinke, P. Talkner, and Y. Dodge, Turbulent cascades in foreign exchange markets, Nature 381, 767 (1996).
  * Chialvo [2010] D. R. Chialvo, Emergent complex neural dynamics, Nature Physics 6, 744 (2010).
  * Poole  _et al._ [2016] B. Poole, S. Lahiri, M. Raghu, J. Sohl-Dickstein, and S. Ganguli, Exponential expressivity in deep neural networks through transient chaos, in  _Advances in Neural Information Processing Systems_, Vol. 29, edited by D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett (Curran Associates, Inc., 2016).
  * Sompolinsky  _et al._ [1988] H. Sompolinsky, A. Crisanti, and H. J. Sommers, Chaos in random neural networks, Phys. Rev. Lett. 61, 259 (1988).
  * Kadmon and Sompolinsky [2015] J. Kadmon and H. Sompolinsky, Transition to chaos in random neuronal networks, Phys. Rev. X 5, 041030 (2015).
  * Aljadeff  _et al._ [2015] J. Aljadeff, M. Stern, and T. Sharpee, Transition to chaos in random networks with cell-type-specific connectivity, Phys. Rev. Lett. 114, 088101 (2015).
  * Stern  _et al._ [2014] M. Stern, H. Sompolinsky, and L. F. Abbott, Dynamics of random neural networks with bistable units, Phys. Rev. E 90, 062710 (2014).
  * Rajan  _et al._ [2010] K. Rajan, L. F. Abbott, and H. Sompolinsky, Stimulus-dependent suppression of chaos in recurrent neural networks, Phys. Rev. E 82, 011903 (2010).
  * Wainrib and Touboul [2013] G. Wainrib and J. Touboul, Topological and dynamical complexity of random neural networks, Phys. Rev. Lett. 110, 118101 (2013).
  * Bertschinger and Natschläger [2004] N. Bertschinger and T. Natschläger, Real-Time Computation at the Edge of Chaos in Recurrent Neural Networks, Neural Computation 16, 1413 (2004).
  * Shih  _et al._ [2020] C.-T. Shih, Y.-J. Lin, C.-T. Wang, T.-Y. Wang, C.-C. Chen, T.-S. Su, C.-C. Lo, and A.-S. Chiang, Diverse community structures in the neuronal-level connectome of the drosophila brain, Neuroinformatics 18, 267 (2020).
  * Scheffer  _et al._ [2020] L. K. Scheffer, C. S. Xu, M. Januszewski, Z. Lu, S.-y. Takemura, K. J. Hayworth, G. B. Huang, K. Shinomiya, J. Maitlin-Shepard, S. Berg, _et al._ , A connectome and analysis of the adult Drosophila central brain, eLife 9, e57443 (2020).
  * Buzsáki and Mizuseki [2014] G. Buzsáki and K. Mizuseki, The log-dynamic brain: how skewed distributions affect network operations, Nature Reviews Neuroscience 15, 264 (2014).
  * Perez-Nieves  _et al._ [2021] N. Perez-Nieves, V. C. H. Leung, P. L. Dragotti, and D. F. M. Goodman, Neural heterogeneity promotes robust learning, Nature Communications 12, 5791 (2021).
  * Martin  _et al._ [2021] C. H. Martin, T. S. Peng, and M. W. Mahoney, Predicting trends in the quality of state-of-the-art neural networks without access to training or testing data, Nature Communications 12, 4122 (2021).
  * Benzi  _et al._ [1984] R. Benzi, G. Paladin, G. Parisi, and A. Vulpiani, On the multifractal nature of fully developed turbulence and chaotic systems, Journal of Physics A: Mathematical and General 17, 3521 (1984).
  * Evers and Mirlin [2008] F. Evers and A. D. Mirlin, Anderson transitions, Rev. Mod. Phys. 80, 1355 (2008).
  * Anderson [1958] P. W. Anderson, Absence of diffusion in certain random lattices, Phys. Rev. 109, 1492 (1958).
  * Abou-Chacra  _et al._ [1973] R. Abou-Chacra, D. J. Thouless, and P. W. Anderson, A selfconsistent theory of localization, Journal of Physics C: Solid State Physics 6, 1734 (1973).
  * Parisi  _et al._ [2019] G. Parisi, S. Pascazio, F. Pietracaprina, V. Ros, and A. Scardicchio, Anderson transition on the bethe lattice: an approach with real energies, Journal of Physics A: Mathematical and Theoretical 53, 014003 (2019).
  * Mastrogiuseppe and Ostojic [2018] F. Mastrogiuseppe and S. Ostojic, Linking connectivity, dynamics, and computations in low-rank recurrent neural networks, Neuron 99, 609 (2018).
  * Huang  _et al._ [2019] C. Huang, D. A. Ruff, R. Pyle, R. Rosenbaum, M. R. Cohen, and B. Doiron, Circuit models of low-dimensional shared variability in cortical networks, Neuron 101, 337 (2019).
  * He [2014] B. J. He, Scale-free brain activity: past, present, and future, Trends in Cognitive Sciences 18, 480 (2014).
  * Bernacchia  _et al._ [2011] A. Bernacchia, H. Seo, D. Lee, and X.-J. Wang, A reservoir of time constants for memory traces in cortical neurons, Nature Neuroscience 14, 366 (2011).
  * Keup  _et al._ [2021] C. Keup, T. Kühn, D. Dahmen, and M. Helias, Transient chaotic dimensionality expansion by recurrent networks, Phys. Rev. X 11, 021064 (2021).
  * Helias and Dahmen [2020] M. Helias and D. Dahmen,  _Statistical Field Theory for Neural Networks_ (Springer International Publishing, Cham, 2020).
  * van Meegen  _et al._ [2021] A. van Meegen, T. Kühn, and M. Helias, Large-deviation approach to random recurrent neuronal networks: Parameter inference and fluctuation-induced transitions, Phys. Rev. Lett. 127, 158302 (2021).
  * Tarquini  _et al._ [2016] E. Tarquini, G. Biroli, and M. Tarzia, Level statistics and localization transitions of lévy matrices, Phys. Rev. Lett. 116, 010601 (2016).
  * Kuśmierz  _et al._ [2020] L. Kuśmierz, S. Ogawa, and T. Toyoizumi, Edge of chaos and avalanches in neural networks with heavy-tailed synaptic weight distribution, Phys. Rev. Lett. 125, 028101 (2020).
  * Wardak and Gong [2021] A. Wardak and P. Gong, Fractional diffusion theory of balanced heterogeneous neural networks, Phys. Rev. Research 3, 013083 (2021).
  * Schertzer  _et al._ [2001] D. Schertzer, M. Larchevêque, J. Duan, V. V. Yanovsky, and S. Lovejoy, Fractional fokker–planck equation for nonlinear stochastic differential equations driven by non-gaussian lévy stable noises, Journal of Mathematical Physics 42, 200 (2001).
  * Note [1] For clarity, we refer to $gJ\operatorname{diag}_{j}(\phi^{\prime}(h_{j}))$ as the Jacobian matrix in the remainder of the Letter as its eigenvectors remain unchanged and its eigenvalues are shifted one unit to the right.
  * Bordenave  _et al._ [2011] C. Bordenave, P. Caputo, and D. Chafaï, Spectrum of non-hermitian heavy tailed random matrices, Communications in Mathematical Physics 307, 513 (2011).
  * Note [2] See Supplemental Material.
  * Metz  _et al._ [2019] F. L. Metz, I. Neri, and T. Rogers, Spectral theory of sparse non-hermitian random matrices, Journal of Physics A: Mathematical and Theoretical 52, 434003 (2019).
  * Rajan and Abbott [2006] K. Rajan and L. F. Abbott, Eigenvalue spectra of random matrices for neural networks, Phys. Rev. Lett. 97, 188104 (2006).
  * Wei [2012] Y. Wei, Eigenvalue spectra of asymmetric random matrices for multicomponent neural networks, Phys. Rev. E 85, 066116 (2012).
  * Solórzano  _et al._ [2021] A. Solórzano, L. F. Santos, and E. J. Torres-Herrera, Multifractality and self-averaging at the many-body localization transition, Phys. Rev. Research 3, L032030 (2021).
  * Stringer  _et al._ [2019] C. Stringer, M. Pachitariu, N. Steinmetz, M. Carandini, and K. D. Harris, High-dimensional geometry of population responses in visual cortex, Nature 571, 361 (2019).
  * Metz and Neri [2021] F. L. Metz and I. Neri, Localization and universality of eigenvectors in directed random graphs, Phys. Rev. Lett. 126, 040604 (2021).

See pages ,,1,,2,,3,,4,,5,,6,,7,,8,,9 of supp_20220211.pdf