Spiking neural networks for nonlinear regression — TEI-tagged view

Alexander HenkespersName, Jason K. EshraghianpersName, Henning WesselspersName (2024date). “Spiking neural networks for nonlinear regression.” Royal Society Open Sciencetitle·j.

Object of study — spiking neural networkterm·ai architectures for accurate, energy-efficient nonlinear and history-dependent regressionterm·key

Research question — how to decode binary spikesterm·ai into real numbers so spiking networks can do regressionterm·key without dense layers

Contribution — a membrane-potential decoding layer with population voting, used in LIFterm·ai, recurrent LIFterm·ai and spiking LSTMterm·key regressionterm·key networks

Full text — tagged

Spiking neural networksterm·ai for nonlinear regressionterm·key Abstract Spiking neural networksterm·ai (SNNterm·ai), also often referred to as the third generation of neural networksterm·ai, carry the potential for a massive reduction in memory and energy consumptionterm·var over traditional, second-generation neural networks. Inspired by the undisputed efficiency of the human brain, they introduce temporal and neuronal sparsityterm·ai, which can be exploited by next-generation neuromorphic hardwareterm·key. Energy efficiencyterm·key plays a crucial role in many engineering applications, for instance, in structural health monitoring. Machine learning in engineering contexts, especially in data-driven mechanics, focuses on regressionterm·key. While regressionterm·key with SNNterm·ai has already been discussed in a variety of publications, in this contribution, we provide a novel formulation for its accuracy and energy efficiencyterm·key. In particular, a network topology for decoding binary spiketerm·ai trains to real numbers is introduced, using the membrane potentialterm·key of spiking neurons. Several different spiking neural architectures, ranging from simple spiking feed-forward to complex spiking long short-term memory neural networks, are derived. Since the proposed architectures do not contain any dense layers, they exploit the full potential of SNNterm·ai in terms of energy efficiencyterm·key. At the same time, the accuracy of the proposed SNNterm·ai architectures is demonstrated by numerical examples, namely different material models. Linear and nonlinear, as well as history-dependent material models, are examined. While this contribution focuses on mechanical examples, the interested reader may regress any custom function by adapting the published source code.

Keywords: artificial neural networks, spiking neural networksterm·ai, regressionterm·key, neuromorphic computing Introduction In recent years, artificial neural networks (ANN) have gained much attention in the computational engineering sciences and applied mathematics due to their flexibility and universal approximation capabilities, both for functions [ 1 , 2 ] and operators [ 3 ]. Their outstanding but surprising generalizability capabilities are yet to be understood [ 4 ]. Their advantages have been used in a variety of applications, including fluid dynamics [ 5 – 9 ], solid mechanics [ 10 – 13 ], micromechanics [ 14 – 17 ], material parameter identification [ 18 – 22 ], constitutive modelling [ 23 – 26 ], fracture mechanics [ 27 , 28 ], microstructure generation [ 29 – 31 ], contact problems [ 32 – 35 ], heat transfer [ 36 – 39 ], uncertainty quantification [ 40 – 43 ], among numerous others. Refer to [ 44 – 46 ] for review publications.

Despite the success of ANNs, several problems arise alongside their use, such as the need for high-frequency memory accessterm·ai, which leads to high computational power demandterm·ai [ 47 , 48 ]. This results in huge costs for training and often makes it preferable to run inference in remote servers during deployment. In general, ANNs are most often trained on graphics processing units (GPUs), whose energy consumptionterm·var is problematic in embedded systems (e.g. sensor devices) as is required in automotive and aerospace applications [ 49 ] or structural health monitoring (SHM). Furthermore, high latency during prediction time can arise where acceleration or parallelization is not available. Figure 1 illustrates a potential use case of neuromorphic hardwareterm·key in SHM. Using concepts from data-driven mechanics, sensor data collected during service life may be interpreted with the aid of high-dimensional (computational) physical models, which so far have only been used to design infrastructures or products. Thanks to active research, the computational efficiency of such advanced monitoring concepts is rapidly increasing (e.g. [ 50 ], for a recent contribution). A bottleneck for practical online applications, however, remains the aforementioned hardware constraints, namely, energy demand and latency. Regressionterm·key using spiking neural networksterm·ai (SNNterm·ai) on neuromorphic hardwareterm·key may contribute to overcome persisting limitations.

SNN-based regressionterm·key offers high potential in many energy-critical applications, such as SHM. We expect that, thanks to the advent of data-driven modelling, it will become possible to interpret sensor data directly where it is recorded on neuromorphic hardwareterm·key as part of embedded devices. Current work is concerned with the solution of partial differential equations on neuromorphic hardwareterm·key, which will ultimately allow to evaluate and update computational models online.

Originally motivated by the human brain, today’s traditional ANN architectures are an oversimplification of biology, relying on dense matrix multiplication. From a numerical and computational hardware point of view, dense matrix multiplication is often suboptimal. Sparsityterm·ai is thought to be favourable as it reduces dependence on memory access and data communication [ 51 ]. In contrast, the human brain is much more efficient, where neurons are considered to be sparsely activated [ 52 ]. This stems from the fact that the brain uses sparse electronic signals for information transmission instead of dense activations. This leads to remarkable capabilities by using only approximately 10–20 W of energy. One attempt to overcome these drawbacks of ANNs is to introduce the information transmission mechanism of biological neurons into network architectures. These networks are called SNNterm·ai due to the electronic impulses or spikesterm·ai used for communication between neurons [ 53 ]. This leads to sparse activations, which can be efficiently exploited by neuromorphic hardwareterm·key, such as Loihi [ 54 ], SpiNNaker [ 55 ] and TrueNorth [ 56 ]. It has been shown that these specialized hardware chips can reduce the energy consumptionterm·var of neural network-based processes by factors of up to ×1000 [ 54 , 57 – 60 ]. Apart from their energy efficiencyterm·key in prediction, recent attempts to increase the training efficiency of SNNterm·ai can be found in [ 61 – 64 ].

What was classically in the domain of neuroscientists recently has been investigated in the context of deep learning, for example, the adoption of SNNsterm·ai to supervised learning as popularized with traditional ANNs in frameworks, such as TensorFlow [ 65 ] and PyTorch [ 66 ], resulting in similar frameworks for spiking deep learning like snnTorch [ 67 ]. Some applications of spiking deep learning include image processing using a spiking ResNet [ 68 ] and temporal data processing using spiking long short-term memory neural network (SLSTMterm·ai) variants [ 69 , 70 ]. A combination of spiking convolutional neural networks and LSTMs was proposed in Yang et al . [ 71 ]. SNNsterm·ai have been used for image segmentation [ 72 ] and localization [ 73 , 74 ].

To the best of the authors’ knowledge, the full potential of energy-efficient regressionterm·key modelling using SNNsterm·ai has up to now not yet been fully exploited. In Iannella and Back [ 75 ], an architecture using inter-spike interval temporal encoding has been proposed, where the functions learned were limited to piecewise constant functions. In Gehrig et al . [ 76 ], a SNNterm·ai was used for the regressionterm·key of angular velocities of a rotating event camera. Building on these results, Rançon et al. [ 77 ] proposed a SNNterm·ai for depth reconstruction. Both works rely on fully connected decoders for real-valued output, which increases the energy demand compared with an architecture relying only on SNNterm·ai predictions. This is due to the increased number of dense matrix multiplications in fully connected layers. In Kahana et al . [ 78 ], a DeepONet [ 79 ] using SNNterm·ai was proposed, which used a floating point decoding scheme to regress simple one-dimensional functions. This decoding scheme led to non-smooth and staircase-like predictions. In Shrestha and Orchard [ 80 ], gradient descent was applied to learn spiketerm·ai times. Here, as in the present contribution, the membrane potentialterm·key is used, but to tweak the occurrence of the first spikesterm·ai in the SNNterm·ai. In EshraghianpersName et al . [ 81 ], classification problems using the membrane potentialterm·key in the context of memristor-based hardware were explored, but no regressionterm·key tasks. Recently, a SNNterm·ai has been proposed in a computational mechanics context in Tandale and Stoffel [ 82 ], which builds up on our approach developed in the preprint [ 83 ]. However, again fully connected layers were added, which counteract the energy efficiencyterm·key.

The focus of the present work lies on network architectures that come along without any fully connected layers and are therefore optimally suited for neuromorphic hardwareterm·key, which is specifically designed for SNNsterm·ai. The proposed regressionterm·key approach works without fully connected dense encoders/decoders by using only spiking layers. The output of the proposed networks is real-valued, continuous functions and can represent highly nonlinear, history-dependent data accurately.

As regressionterm·key problems are omnipresent in engineering sciences, a flexible and broadly applicable framework would enable SNNsterm·ai to be used in a variety of engineering applications and further unfold the potential of neuromorphic hardwareterm·key. To this end, the present study aims towards the following key contributions: — Introduction of spiking neural networkterm·ai s. Concise introduction of this emerging technique. Open source benchmark code for the research community, fostering further developments of SNNterm·ai in mechanics and applied mathematics (see also Data accessibility section).

— History-dependent regressionterm·key framework . We demonstrate that SNNsterm·ai can model systems that exhibit hysteresis, thus history-dependent behaviour. We present a flexible framework to use SNNsterm·ai in complex regressionterm·key tasks, demonstrated by means of a nonlinear material model, namely isotropic hardening plasticity.

— Efficiency, sparsityterm·ai and latency . We benchmark our SNNterm·ai for neuromorphic hardwareterm·key in terms of memory and energy consumptionterm·var as compared with non-spiking equivalent networks, demonstrating that our approach is much more efficient with respect to memory and power consumption, making neural networks more sustainable. This is enabled by a membrane potential-based decoding scheme, which does not rely on fully connected decoders. While in our study an emulator is used, deployment on neuromorphic hardwareterm·key allows highly efficient usage in embedded environments.

The present work intends to introduce this important novel technique to the community of computational mechanics and applied mathematics. To concentrate on the novelties and to keep the presentation concise, we restrict ourselves to one-dimensional regressionterm·key problems. However, the framework is not restricted to single-variable regressionterm·key and is easily applicable to a multi-variable regressionterm·key.

The remainder of this article is structured as follows. In §2, the basic notations of SNNsterm·ai are derived from traditional ANNs. A simple spiking counterpart to the classical densely connected feed-forward neural network is introduced. After that, our regressionterm·key SNNterm·ai topology is proposed. Starting with a linear regressionterm·key example (linear elasticity), we point out the problems arising in SNNterm·ai regressionterm·key. This basic architecture is extended towards recurrent feedback loops in §3. The ability of these recurrent SNNsterm·ai is showcased on a nonlinear material model. To process history-dependent regressionterm·key tasks with dependencies over a large number of time steps, a spiking LSTMterm·key is introduced in §4. An application to a history-dependent plasticity model shows that SNNsterm·ai can achieve similar accuracy as their traditional counterparts while being much more efficient. The paper closes with a conclusion and an outlook towards future research directions in §5. For the code accompanying this manuscript, see the Data accessibility section at the end of this manuscript.

Spiking neural networksterm·ai for regressionterm·key SNNsterm·ai are considered to be the third generation of neural networksterm·ai. While the first generation was restricted to shallow networks, the second generation is characterized by deep architectures. A broad use of second-generation neural networks has been enabled by the availability of automatic differentiation and software frameworks such as TensorFlow [ 65 ]. To introduce SNNterm·ai, we compare them with their well-known second-generation counterparts.

Second-generation neural networks An ANN is a parametrized, nonlinear function composition. The universal function approximation theorem [ 1 ] states that arbitrary Borel measurable functions can be approximated with ANNs. There are several different architectures for ANNs, for example, feed-forward, recurrent or convolutional networks, which can be found in standard references such as [ 84 – 88 ]. Following Hauser [ 89 ], most ANN formulations can be unified. An ANN N , more precisely, a densely connected feed-forward neural network , is a function from an input space R d x to an output space R d y , defined by a composition of nonlinear functions h ( l ) , such that Here, x denotes an input vector of dimension d x and y an output vector of dimension d y . The nonlinear functions h ( l ) are called layers and define an l -fold composition, mapping input vectors to output vectors. Consequently, the first layer h ( 0 ) is defined as the input layer and the last layer h ( n L ) as the output layer , such that The layers h ( l ) between the input and output layers, called hidden layers , are defined as where h η ( l ) is the n -th neural unit of the l -th layer h ( l ) , n u denotes the total number of neural units per layer, W η ( l ) is the weight vector of the n -th neural unit in the l -th layer h ( l ) and h ( l − 1 ) is the output of the preceding layer, where bias terms are absorbed [ 86 ]. Furthermore, ϕ ( l ) : R → R is a nonlinear activation function. All weight vectors W η ( l ) of all layers h ( l ) can be gathered in a single expression, such that where θ inherits all parameters of the ANN N ( x ) from equation (2.1) . Consequently, the notation N ( x ; θ ) emphasizes the dependency of the outcome of an ANN on the input on one hand and the current realization of the weights on the other hand. The specific combination of layers h ( l ) from equation (2.3) , neural units h η ( l ) and activation functions ϕ ( l ) from equation (2.3) is called the topology of the ANN N ( x ; θ ) . The weights θ from equation (2.4) are typically found by gradient-based optimization with respect to a task-specific loss function [ 85 ]. An illustration of a densely connected feed-forward ANN is shown in figure 2 .

Densely connected feed-forward neural network topology of an ANN N ( x ; θ ) as described in equation (2.1) .

It can be seen that the ANN described in equation (2.1) takes an input and produces an output, one at a time. If history-dependent input and output data and are considered, the formulation of the hidden layers reads where the time component is discrete. This can be understood as processing each discrete-time slice of the input vector of the preceding layer h t = 0 ( l − 1 ) → h t = 1 ( l − 1 ) → , . . . , → h t = d t ( l − 1 ) sequentially, where the weights W η ( l ) are shared over all time steps. At this stage, the formulation in equation (2.5) is purely notational, as there is no connection of the weights through different time steps.

Spiking neural networksterm·ai A SNNterm·ai can be seen as a history-dependent ANN, which introduces memory effects by means of biologically inspired processes. Standard works in theoretical neuroscience include the studies of Dayan and Abbott [ 90 ], Izhikevich [ 91 ] and Gerstner et al . [ 92 ]. Several overviews of SNNsterm·ai with respect to deep learning can be found in [ 93 – 95 ].

To this end, the activation function ϕ ( l ) in equation (2.5) can be formulated as with where U η , t ( l ) is the membrane potentialterm·key of the η -th neural unit at time t , U thr , η ( l ) denotes the membrane threshold, β η ( l ) is the membrane potentialterm·key decay rate and W η ( l ) h t ( l − 1 ) is the standard ANN weight multiplied with the preceding layer of the current time step, respectively, see equation (2.5) . Basically, the SNNterm·ai activation restricts the neural unit to output discrete pulses ( ϕ spk = 1 ) if the membrane threshold is reached by the time-evolving membrane potentialterm·key, or to remains silent ( ϕ spk = 0 ) . These pulses are called spikesterm·ai , see Figure 3 . The last summand in equation (2.6) , − ϕ spk , t − 1 ( l ) U thr , η ( l ) , is called the reset mechanism and resets the membrane potentialterm·key by the threshold potential once a spiketerm·ai is emitted. The membrane threshold and membrane potentialterm·key decay rate can be optimized during training, such that the optimization parameters of a SNNterm·ai are Spiking neuron dynamics. Input spikesterm·ai ϕ s p k i , 0 trigger changes in the membrane potentialterm·key U 1 , which when sufficiently excited beyond a threshold U t h r causes the neuron to emit an output spiketerm·ai ϕ s p k j , 1 .

The SNNterm·ai formulation in equation (2.7) is called the leaky integrate and fireterm·ai (LIFterm·ai) neuron model, and is one of the most widely used models in spike-based deep learning. It can be seen as the baseline SNNterm·ai and plays a similar role as densely connected feed-forward ANN in classical deep learning. Equation (2.7) can also be interpreted as the explicit forward Euler solution of an ordinary differential equation, describing the time variation of the membrane potentialterm·key, see EshraghianpersName et al . [ 67 ] for details.

The main difference between SNNsterm·ai and classical ANNs lies in the way information is processed and propagated through the network from neuron to neuron. In standard ANNs, inputs, hidden layers and output vectors are handled via dense matrices. In SNNterm·ai, sparsityterm·ai is introduced by using spikesterm·ai, which are single events expressed via a Dirac delta function or a discrete pulse in continuous or discrete settings, respectively. A group of spikesterm·ai over time is called a spiketerm·ai train i = [ i t , t = 0 , . . . , n t ] . To this end, a spiking neuron is subjected to a spiketerm·ai train over a time interval, consisting of spikesterm·ai (1) or zero input (0). The membrane potentialterm·key U η , t ( l ) is modulated with incoming spikesterm·ai i t . In the absence of input spikesterm·ai, the membrane voltage decays over time due to the membrane decay rate β η ( l ) . The absence of spikesterm·ai introduces sparsityterm·ai because, in every time step, the neural unit output is constrained to either 0 or 1. This fact can be exploited on neuromorphic hardwareterm·key , where memory and synaptic weights need only be accessed if a spiketerm·ai is apparent. Otherwise, no information is transmitted. In contrast, conventional ANNs do not leverage sparsely activated neurons, and most deep learning accelerators, such as GPUs or tensor processing units (TPUs), are correspondingly not optimized for it.

Spiking neural networksterm·ai training Unfortunately, the spiking activation ϕ spk , t ( l ) in equation (2.6) is non-differentiable. To use the backpropagation algorithm from standard ANNs, the activation is replaced using a surrogate gradientterm·mth during the backward pass. Several different formulations have been proposed (e.g. [ 51 , 96 , 97 ]). In this work, the arcus tangent surrogate activation from Fang et al . [ 98 ] is used, for some input x . The surrogate ϕ surr ( x ) is continuously differentiable and preserves the gradient dynamics of the network. Thus, for training using backpropagation and its variants, ϕ surr ′ is employed. Illustrations can be found in figures 4 and 5 .

Computational graph of an unrolled SNNterm·ai. Forward pass.

Computational graph of an unrolled SNNterm·ai. Backward pass.

Network topology The key question for using SNNterm·ai in regressionterm·key is how to transform real input values into binary spikesterm·ai and binary spiketerm·ai information at the output layer back to real numbers. The former task is called spike encodingterm·mth , whereas the latter is called spike decodingterm·mth . Popular forms of encoding include rate encoding , latency encoding, delta modulation , among others. Similarly, different decoding strategies exist, such as rate decoding and latency decoding . An illustration of various encoding and decoding strategies is shown in figure 6 . Refer to EshraghianpersName et al . [ 67 ] for an overview and detailed description. In this work, a constant current injection is chosen for the encoding part, whereas a novel population voting on membrane potentialterm·key approach is chosen for the decoding part. It will be demonstrated throughout the numerical examples that, with the proposed spike decodingterm·mth scheme, accurate and efficient SNNterm·ai regressionterm·key is made possible.

A sample of spike-based encoding and decoding strategies. ( a ): Real-valued inputs are encoded into spikesterm·ai by means of different strategies [ 67 ], for example: (i) high intensities or large values result in a large number of spikesterm·ai (top left); (ii) high intensities or large values result in early spiketerm·ai firing (centre left); and (iii) delta modulation where spikesterm·ai are produced for positive gradients of the input function (bottom left). ( b ): In classification, the predicted class is determined via (i) the number of spikesterm·ai (top right) or (ii) the first occurrence of spikesterm·ai (bottom right). Regressionterm·key is introduced in §2.2.

All network topologies used in the upcoming numerical examples follow a general scheme, which is flexible and suited for regressionterm·key tasks. First, the real input x t is provided as a constant input to the first layer h t ( 0 ) , for all time steps t , such that Then, several SNNterm·ai layers h t ( l ) follow, where the exact formulation will be given for every numerical example. The output of the last spiking layer h t ( n L ) is transformed into a decoding layer h t dec , which takes the membrane potentialterm·key of every time step as input and outputs real numbers which is essentially the formulation of equation (2.7) , where no spikesterm·ai and reset mechanisms are used. The transformed values are then transferred to a ‘population voting layer’, where the output of all neurons of the decoding layer are averaged to give real numbers. This results in where n o denotes the number of neurons in the population voting layer and again, no spikesterm·ai or reset mechanisms are used.

The final spiking regressionterm·key topology network S can be written as To summarize, information flows in the following form: — The input layer h t const = x t equation (2.10) takes constant current (real numbers x t ).

— The input is then transformed into binary spikesterm·ai in the spiking layers h t ( l ) ( equation (2.6) ) — and transformed back into real numbers in the translation layer h t dec Equation (2.11) ().

— The output is formed in the population layer h t pop Equation (2.12) .

A graphical interpretation is given in figure 7 .

Topology of the spiking regressionterm·key network is introduced in §2.1.

For all the following numerical examples, the AdamW optimizer from Loshchilov and Hutter [ 99 ] is used. The parameters are set as follows: learning rate α = 1 × 10 − 3 , exponential decay rates for the first- and second-moment estimates β 1 = 0.9 and β 2 = 0.999 , respectively, weight decay λ = 0.01 . The training was carried out on a Nvidia GeForce RTX 3090 GPU using snnTorch [ 67 ] and PyTorch [ 66 ]. In this work, the mean relative error E is used, which is defined as for some input ∙ and baseline ∙ ^ . If the error is reported for all time steps, ∙ is a vector containing the values of all time steps. If the error is reported for the last time step, ∙ is equal to the last component of the corresponding history-dependent vector.

Numerical experiment: linear regressionterm·key We first study the performance of the proposed LIFterm·ai topology in a simple linear regressionterm·key problem. To this end, the general model described in §2.1 with LIFterm·ai defined in equation (2.7) is used, resulting in the following network topology.

To begin with, a simple linear elastic material model with strains in the range of ε = [ 0 , 0.001 ] and fixed Young’s modulus E = 2.1 × 10 5 MPa is considered, such that the resulting stress σ is The training data consists of strain input, uniformly sampled in the interval ε = [ 0 , 0.001 ] , and stress output calculated according to equation (2.16) . Three datasets are generated from equation (2.16) , namely, training set, validation set and test set consisting each of n train = n val = n test = 1024 samples. All three datasets are standardized using the mean and standard deviation from the training set. The batch size is chosen as n batch = 1024 . The number of neurons n u is chosen as n u = 128 and is kept constant over all layers. The training is carried out for 2 × 10 3 epochs. The model performing best on the validation set is chosen for subsequent evaluations. The mean relative error accumulated over all time steps and the mean relative error of the last time step with respect to the test set are reported.

Although the linear regressionterm·key problem is a time-independent problem, SNNterm·ai has an inherent time dependency. Instead of time, the SNNterm·ai is thus trained on a sequence of strain steps Δ ε . We investigate the effect of the number of strain (time) steps on the prediction accuracy, with results illustrated in figure 8 .

Elasticity—error of all strain (time) steps. Mean relative error for all strain steps with respect to the different total number of strain steps. The error is rising for a larger number of strain steps.

It can be seen that the mean relative error is lowest for d t = 5 strain steps. For d t = 2 , the error is larger. This could be caused by a lack of a sufficient number of strain steps for the neuron dynamics to effectively be calculated. It can be understood as a failure due to too large strain steps in the explicit stepping scheme in equation (2.7) . Clearly, the highest error can be observed for d t = 100 strain steps. In contrast, as depicted in figure 9 , the error at the last strain step is lowest for d t = 100 strain steps. To illustrate the cause, the prediction of the network for two different samples, one for d t = 5 and one for d t = 100 strain steps are shown in figures 10 and 11 , respectively. While good agreement on the endpoints is apparent, fluctuation during the rest of the strain steps causes the error to rise. Seemingly, the LIFterm·ai has difficulties regressing a large number of strain steps. This could be caused by the lack of recurrent connections in the LIFterm·ai formulation from equation (2.15) , where history dependency is only weakly included in the form of the membrane potentialterm·key. Therefore, recurrent LIFsterm·ai will be introduced in §3.

Elasticity—error at last strain (time) step. Mean relative error for the last strain step with respect to the different total number of strain steps. The error is converging for a larger number of strain steps.

Elasticity—prediction in five strain (time) steps. Prediction of the LIFterm·ai from equation (2.15) for d t = 5 . Deviations from the true solution can be observed in the middle part.

Elasticity—prediction in 100 strain (time) steps. Prediction of the LIFterm·ai from equation (2.15) for d t = 100 . Fluctuations around the true solution can be observed.

Remark : The seemingly simple linear regressionterm·key task provides a challenge for SNNterm·ai, as effectively an ordinary differential equation has to be fitted to a linear function while relying on binary information transmission and inexact gradients.

Nonlinear regressionterm·key using recurrent leaky integrate and fireterm·ai In order to counter the problems of vanishing information for a large number of strain (time) steps encountered in the preceding section, a recurrent SNNterm·ai architecture is proposed (§3.1). Its performance is demonstrated by means of a numerical example in §3.2.

Recurrent leaky integrate and fireterm·ai The standard LIFterm·ai is a feed-forward neuron, such that information is flowing unidirectionally in the form of spikesterm·ai. By adding a feedback loop, a recurrent LIFterm·ai (RLIFterm·ai) can be formulated, which builds on the standard recurrent neural network (RNN) formulation. This enables the network to use relationships along several time steps for the prediction of the current time step. It was shown in Pascanu et al . [ 100 ] that recurrent loops can retain information for a relatively longer number of time steps when compared with their non-recurrent counterparts.

Here, the formulation of the hidden layer in equation (2.5) includes additional recurrent weights V η ( l ) , such that In this RNN, the influence of the preceding time step is explicitly included by means of additional recurrent weights V η ( l ) . The resulting set of trainable parameters reads The RNN formulation can be included in the LIFterm·ai formulation from equation (2.7) to obtain an RLIFterm·ai, such that where U η , t ( l ) is again the membrane potentialterm·key of the η th neural unit at time t , U thr , η ( l ) denotes the membrane threshold, β η ( l ) is the membrane potentialterm·key decay rate and W η ( l ) h t ( l − 1 ) is the standard ANN weight multiplied with the preceding layer at the current time step. Additionally, V η ( l ) denotes the recurrent weights from equation (3.1) . This leads to the following set of trainable parameters Numerical experiment: Ramberg–Osgood The performance of the RLIFterm·ai is investigated in nonlinear function regressionterm·key. As a first test case, the nonlinear Ramberg–Osgood power law for modelling history-independent plasticity is considered, in which stress σ and strain ε are related via Here, ε is the infinitesimal, one-dimensional elastic strain, σ denotes the one-dimensional Cauchy stress, E is Young’s modulus, n is a material constant describing the hardening behaviour of plastic deformation and σ Y is the yield strength of the material. Note that this plasticity model is only suited for a single loading direction and does not incorporate accumulation of plastic strain.

To this end, the general model described in §2.1 using RLIFterm·ai defined in equation (3.3) is used. For a given strain history ε t , the stress response for different values of yield stress σ Y is predicted with the following parametric architecture The training data consists of yield strength σ Y as input for fixed strains in the interval ε = [ 0 , 0.01 ] for d t = 20 strain steps. The yield strength is uniformly sampled in the interval σ Y = [ 100 , 500 ] MPa, and the stress output is calculated according to equation (3.5) . The Young’s modulus is chosen as E = 2.1 × 10 5 MPa and n = 10 . Three datasets are generated, namely, training set, validation and test set with n train = n val = n test = 1024 samples, respectively. All three sets are standardized using the mean and standard deviation from the training set. The batch size is chosen as n batch = 1024 . The number of neurons n u is chosen as n u = 128 and is kept constant over all layers. The training is carried out for 5 × 10 3 epochs. The model performing best on the validation set is chosen for subsequent evaluations. The mean relative error and the mean relative error of the last strain step with respect to the test set are reported.

In figure 12 , different stress–strain curves are depicted for different yield strength values, obtained by solving equation (3.5) for the stress σ with a classical Newton–Raphson method. The results of five different samples, randomly chosen from the n t e s t = 1024 test samples, can be seen in figure 13 . For the test set, a mean relative error for all strain steps of 8.7934 × 10 − 2 and a mean relative error for the last strain step of 8.0200 × 10 − 2 are obtained. The predictions on these five samples are more accurate than would be suspected from the mean relative error. The cause can be found in figure 14 , where the mean relative error for all strain steps is plotted for every sample of the test set. It can be observed that a small number of samples has a much higher error than the rest, which impacts the error measure. This is caused by the purely data-driven nature of the experiment and can be tackled with approaches introduced in, for example, [ 101 – 103 ] .

Ramberg–Osgood—reference solutions. Stress–strain curves of the Ramberg–Osgood material model for five different values of the yield stress σ Y obtained with Newton–Raphson algorithm.

Ramberg–Osgood—RLIFterm·ai prediction. Prediction of the RLIFterm·ai from equation (3.6) for five different yield strength σ Y sampled from the test set for the nonlinear Ramberg–Osgood plasticity law.

Ramberg–Osgood—RLIFterm·ai test error. The mean relative error over all strain steps for 1024 samples from the test set for the numerical experiment is described in §3.2. It can be seen that some outliers have large error values, resulting in a mean error over all samples of 8.7934 × 10 − 2 . Most samples have a significantly lower error.

Nevertheless, the RLIFterm·ai is able to regress on the varying yield strength σ Y and can predict the resulting nonlinear stress–strain behaviour, as can be seen in the predictions ( figure 13 ). Deviations can be observed around the yield point as well as the endpoints of the curves. To be able to take into account long-term history-dependent behaviour, the RLIFterm·ai formulation will be expanded towards the incorporation of explicit long-term memory in the next section, where a more complex plasticity model is investigated.

History-dependent regressionterm·key using spiking long short-term memory network To extend the limited memory of the RNN in §3, a SLSTMterm·ai (§4.1) is proposed. Previously, SLSTMterm·ai hs been considered for classification problems [ 71 ]. Herein, a novel regressionterm·key approach with a population decoding layer is discussed by means of a history-dependent plasticity model (§4.2).

Spiking long short-term memory network A SLSTMterm·ai is the spiking version of the standard LSTM [ 104 ], where the latter is defined as with where f t denotes the forget gate with sigmoid activation ϕ sigmoid or tangent hyperbolicus activation ϕ tanh and corresponding weights W f , V f with absorbed biases. The same nomenclature holds for the input gate i t , the output gate o t , the cell input c ~ t and the cell state c t with their respective activations and weights. The new cell state c t and the output of the LSTM h t are formed using the Hadamard or point-wise product ⊙ . The parameters of the LSTM are its weights, such that For detailed derivations and explanations of standard LSTM, see, for example, [ 85 , 86 ]. The SLSTMterm·ai can be obtained from the LSTM by using spiketerm·ai activations within the LSTM formulation from equation (4.1) , such that where the output h η , t ( l ) is used to determine if a spiketerm·ai is produced In other words, the output of h η , t ( l ) can be interpreted as the membrane potentialterm·key of the SLSTMterm·ai, such that h η , t ( l ) = U η , t ( l ) . A decay parameter β is not used in this formulation. Rather than using decay to remove information from the cell state c η , t ( l ) , this is achieved by carefully regulated gates. The corresponding optimization parameters of the SLSTMterm·ai are Basically, the cell state c η , t ( l ) acts as long-term memory, just like in the standard LSTM formulation. The communication between layers is handled via spiketerm·ai trains that depend on the membrane potentialterm·key h η , t ( l ) = U η , t ( l ) in equation(4.3) and the activation function ϕ spk , t ( l ) from equation (4.4) .

Numerical experiment: Isotropic hardening using spiking long short-term memory network The following numerical experiments aim to investigate the performance of the proposed SLSTMterm·ai on nonlinear, history-dependent problems. Therefore, a one-dimensional plasticity model with isotropic hardening is investigated. Following Simo and Hughes [ 105 ], the model is defined by where is the additive elastoplastic split of the small-strain tensor ε into a purely elastic part ε el and a purely plastic part ε el ; denotes the elastic stress–strain relationship for the Cauchy stress tensor σ and elastic modulus E ; describes the flow rule and isotropic hardening law with consistency parameter γ and equivalent plastic strain α ; gives the yield condition f ( σ , α ) with hardening modulus K ; denotes the Kuhn–Tucker complementarity conditions; and describes the consistency condition.

In figure 15 , different stress–strain paths are shown for varying strains. Especially long-time dependencies are of interest. To this end, the predictive capabilities of the SNNterm·ai are investigated for inference over d t = 100 strain steps, where the elastoplastic model is evaluated using a classical explicit return-mapping algorithm, see Simo and Hughes [ 105 ].

Isotropic hardening—reference solutions. Five stress–strain curves sampled from the isotropic hardening material model for different maximum strains obtained from equation (4.6) .

The training data consists of strain as input, uniformly sampled in the interval ε = [ 0 , 0.01 ] , and stress as output calculated according to equation (4.6) . The yield stress is chosen as σ Y = 300 MPa, the elastic modulus E = 2.1 × 10 5 MPa and the hardening modulus as 2.1 × 10 4 MPa. Three datasets are generated, namely, a training set with n train = 10 240 samples and a validation set and test set with n val = n test = 1024 samples. All three sets are standardized using the mean and standard deviation from the training set. The batch size is chosen as n batch = 1024 . The training is carried out for 500 epochs. The model performing best on the validation set is chosen for subsequent evaluations. The mean relative error accumulated over all time steps and the mean relative error of the last strain step with respect to the test set are reported.

Influence of hyperparameter The first study investigates the prediction accuracy as a function of (i) the number of output neurons, which participate in the population regressionterm·key outlined in §2.1 and (ii) different capacities of the SLSTMterm·ai in the sense of layer width. To this end, the SLSTMterm·ai defined in equation (4.3) is used, resulting in the following architecture: Multiple simulationsterm·mth with output neurons and hidden layers drawn from the grid n u × n o = [ 16 , 32 , 64 , 128 , 256 ] × [ 16 , 32 , 64 , 128 , 256 ] are carried out. The resulting mean relative error for all strain steps with respect to the test set is shown in figure 16 , whereas the resulting mean relative error of the last strain step with respect to the test set is depicted in figure 17 . A clear convergence behaviour can be observed for the number of hidden neurons n u , where larger numbers of neurons lead to lower errors. For the number of output neurons n o , a tendency can be observed upon convergence with respect to n u . For the largest number of hidden neurons n u = 256 , the mean relative error over all strain steps and the mean relative error of the last strain step get larger for n o = [ 128 , 256 ] output neurons, whereas for n o = [ 16 , 32 , 64 ] the errors are almost the same. The lowest mean relative error for all strain steps is 5.2445 × 10 − 2 for n u = 256 hidden neurons per layer and n o = 64 output neurons. The lowest mean relative error for the last strain steps is 2.8729 × 10 − 3 n u = 256 hidden neurons per layer and n o = 32 output neurons. Again, the seemingly high errors are caused by outliers polluting the average, as also observed in §3.2. The same counter-measures discussed in §3.2 can be applied to prohibit outliers, for example, by enforcing thermodynamic consistency.

Isotropic hardening—error versus width. The mean relative error of the last strain step versus the number of hidden neurons per layer is shown for different numbers of output neurons in the isotropic hardening experiment from §4 using the SLSTMterm·ai from equation (4.7) .

Isotropic hardening—error versus width. The mean relative error of all strain steps versus the number of hidden neurons per layer is shown for different numbers of output neurons in the isotropic hardening experiment from §4 using the SLSTMterm·ai from equation (4.7) .

Energy and memory efficiency For the second experiment, the SLSTMterm·ai using n o = 64 output neurons and n u = 256 hidden neurons per layer are compared with a standard LSTM with an equal number of optimization parameters. The aim of this study is the comparison of the prediction accuracy, but also the difference in memory and energy consumptionterm·var on neuromorphic hardwareterm·key. For both ANN variants to be comparable, the same topology is chosen for the LSTM as for the SLSTMterm·ai, such that where the last two layers are replaced by densely connected conventional feed-forward neural networks. Again, the training was carried out for 5 × 10 3 epochs and the same datasets from the previous experiments are used. The standard LSTM from equation (4.8) reached a mean relative error of 4.8611 × 10 − 2 over all strain steps and a mean relative error of 4.7569 × 10 − 3 for the last strain step. The SLSTMterm·ai from equation (4.7) reached a mean relative error of 9.3832 × 10 − 2 over all strain steps and a mean relative error of 4.0497 × 10 − 3 for the last strain step. The resulting prediction for one strain path is illustrated in figure 18 . Clearly, both networks are able to accurately predict the history-dependent, nonlinear stress–strain behaviour.

Isotropic hardening—LSTM versus SLSTMterm·ai. Prediction of a single load path using the return-mapping algorithm as a reference, the standard LSTM and the spiking LSTMterm·key formulation.

Some deviations from the SLSTMterm·ai can be seen at the beginning of the curve. The dynamics of the spiking formulation result in a higher mean relative error over all strain steps with respect to the LSTM. However, the endpoint has a better fit than the LSTM. This is seen in the lower error at the last strain step.

To assess the potential of interfacing our model in embedded, resource-constrained sensors in the wild, we performed a series of power profiling experiments for our SNNsterm·ai (both using LIFterm·ai neurons and SLSTMsterm·ai) when processed on the Loihi neuromorphic chip [ 54 ]. These results are compared against their non-spiking equivalents on an NVIDIA V100 GPU. Data were extracted using the energy profiler in KerasSpiking v. 0.3.0. The emulation tool takes all non-zero activations, accumulates the operation count of each activation (i.e. the fan-out of the activated neuron) and scales this by the estimated energy cost per operation. Specifically, this operation is a ‘read-and-accumulate’, where the presence of a spiketerm·ai requires all connected synapses to be ‘read’ from memory and accumulated with other terms. Other literature refers to this as ‘SynOps’ or ‘Synaptic Operations’ [ 106 – 110 ]. This is a coarse, but reasonable, approximation for small-scale models as neuromorphic hardwareterm·key skips processing non-zero activations. For large-scale models that require inter-chip data communication, the energy cost of overhead is not accounted for in such a model. For our purposes, with lightweight, low-power models, this is an acceptable approximation. It is also the de facto metric for measuring energy efficiencyterm·key in silico from the NeuroBench [ 110 ] initiative, which aims to find representative benchmarking of neuromorphic models.

The first difference in energy usage is that the spiking implementation is measured in an ‘event-based’ manner, where processing only occurs when a neuron emits a spiketerm·ai. In contrast, a non-spiking network processed on a GPU continuously computes with all activations. The second difference is that SNNsterm·ai require multiple strain steps of a forward pass, whereas their non-spiking counterparts do not (unless the input to the network varies over time). Note that the cost of overhead did not need to be accounted for (i.e. transferring data between devices) because all models fit on a single device.

Each network has been broken up into its constituent layers to measure how much they contribute to energy usage on each device. The total energy consumptionterm·var per forward pass of the non-spiking network on the V100 is 512 nJ, whereas the equivalent SNNterm·ai is 4.25 nJ. This represents a 120× reduction in energy consumptionterm·var. The non-spiking LSTM network consumed 5.7 µJ while the proposed SLSTMterm·ai architecture required 24 nJ, a 238× reduction. Detailed results are summarized in table 1 .

Comparison between spiking and non-spiking forward-pass energy consumptionterm·var and memory usage.

Conclusion and outlook In the present study, a framework for regressionterm·key using SNNsterm·ai was proposed based on a membrane potentialterm·key spiking decoder and a population voting layer. Several numerical examples using different spiking neural architectures investigated the performance of the introduced topology towards linear, nonlinear and history-dependent regressionterm·key problems.

First, a simple feed-forward SNNterm·ai, the LIFterm·ai, was derived from the classical densely connected feed-forward ANN. It was shown that the SNNterm·ai can be seen as a special kind of activation function, which produces binary outputs, so-called spikesterm·ai. These spikesterm·ai are used to propagate information through a possibly deep SNNterm·ai. The spikesterm·ai occur due to the dynamic behaviour of the membrane potentialterm·key inside the neuron, which rises when spikesterm·ai appear at the input and decays over time if no spikesterm·ai appear. If a certain threshold value is reached, the membrane potentialterm·key is reset and the neuron emits a spiketerm·ai itself. This formulation introduces more hyperparameters, which fortunately can be learned during training. The spikesterm·ai introduce sparsityterm·ai in the network, which can be effectively exploited by neuromorphic hardwareterm·key to improve latency, power and memory efficiency. The non-differentiability of the binary spikesterm·ai is circumvented by surrogate gradientsterm·mth during backpropagation.

Next, a network topology was proposed, which decodes binary spikesterm·ai into real numbers, which is essential for all kinds of regressionterm·key problems. A decoding layer takes the membrane potentialsterm·key of all neurons in the last spiking layer and propagates them to a population voting layer, which provides its mean potential resulting in a real number. The proposed topology can be used for arbitrary temporal input and output dimensions. A simple experiment on a linear elastic material model using LIFsterm·ai showed that the proposed topology is able to regress the problem. It was shown that errors are introduced for a large number of strain steps. This problem was overcome by introducing RLIFterm·ai, which extends the LIFterm·ai by recurrent feedback loops. An experiment using a nonlinear Ramberg–Osgood plasticity model showed that the proposed topology using RLIFterm·ai is able to regress varying yield limits accurately. The final extension was concerned with the introduction of explicit long-term memories inspired by the classical LSTM formulation, resulting in a spiking LSTMterm·key. The performance of this SLSTMterm·ai was investigated on a history-dependent isotropic hardening model, where different load paths were accurately regressed. During prediction, the SLSTMterm·ai was able to generalize even better than the LSTM for the final load step. Furthermore, the convergence of the proposed method was shown. Note that an extension towards two- and three-dimensional mechanical problems is equally possible. This will involve learning a functional relationship for each stress component, which can be either achieved by a multi-input/output architecture or by using individual networks for each stress component. This treatment is analogous to second-generation ANN, see for instance [ 111 ], where stress data has been decoupled using Proper orthogonal decomposition (POD).

Power profiling and memory analysis were conducted on the LIFterm·ai and SLSTMterm·ai networks to compare efficiency on neuromorphic hardwareterm·key as against a GPU. The Loihi neuromorphic processor was able to achieve a 120× reduction in energy consumptionterm·var when processing the dense LIFterm·ai network, and the SLSTMterm·ai offered a 238× reduction in energy during inference.

The range of possible future application scenarios enabled by regressionterm·key with SNNterm·ai are manifold. For instance, today’s sensing systems cannot capture all quantities that are relevant for structural health monitoring. In the context of mechanics, displacement and strain are quite easy to assess, but the mechanical stress, which reflects the actual response of structures and materials to deformation, remains a so-called hidden quantity . Physics-informed machine learning offers the potential to reconstruct hidden quantities from data by leveraging information from physical models, given in the form of partial differential equations. It is expected that the developments in the field of neuromorphic hardwareterm·key will foster the development of a new generation of embedded systems, which will ultimately enable control of structures and processes based on partial differential equations.

Acknowledgements A.H. was supported by an ETH Zurich Postdoctoral Fellowship. H.W. and A.H. acknowledge the support of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) in the project DFG 255042459/GRK2075-2: Modelling the constitutional evolution of building materials and structures with respect to aging. H.W. additionally acknowledges funding in the project DFG 501798687: Monitoring data-driven life cycle management with AR based on adaptive, AI-supported corrosion prediction for reinforced concrete structures under combined impacts, which is a subproject of the DFG Priority Program 2388: Hundred plus — Extending the Lifetime of Complex Engineering Structures through Intelligent Digitalization .

Contributor Information Alexander HenkespersName, Email: ahenke@ethz.ch; alexanderhenkes1@gmail.com.

Jason K. EshraghianpersName, Email: jeshragh@ucsc.edu.

Ethics This work did not require ethical approval from a human subject or animal welfare committee.

Data accessibility The code is available at: https://github.com/ahenkes1/HENKESpersName_SNNterm·ai and Zenodo [ 112 ].

Declaration of AI use We have not used AI-assisted technologies in creating this article.

Authors’ contributions A.H.: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, project administration, resources, software, supervision, validation, visualization, writing—original draft, writing—review and editing; J.K.E.: conceptualization, funding acquisition, software, supervision, writing—original draft, writing—review and editing; H.W.: conceptualization, funding acquisition, resources, supervision, visualization, writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Conflict of interest declaration The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Funding No funding has been received for this article.

References Associated Data Data Availability Statement The code is available at: https://github.com/ahenkes1/HENKESpersName_SNNterm·ai and Zenodo [ 112 ].