Methods for applying the Neural Engineering Framework to neuromorphic hardware
Aaron R. Voelker Chris Eliasmith

Freezeout systematics due to the hadron spectrum
    Subhasis Samanta
    December 30, 2023
================================================




We review our current software tools and theoretical methods for applying the Neural Engineering Framework to state-of-the-art neuromorphic hardware.
These methods can be used to implement linear and nonlinear dynamical systems that exploit axonal transmission time-delays, and to fully account for nonideal mixed-analog-digital synapses that exhibit higher-order dynamics with heterogeneous time-constants.
This summarizes earlier versions of these methods that have been discussed in a more biological context (Voelker & Eliasmith, 2017) or regarding a specific neuromorphic architecture (Voelker et al., 2017).




§ INTRODUCTION


This report is motivated by the recent influx of neuromorphic computing architectures including SpiNNaker <cit.>, Neurogrid <cit.>, Brainstorm <cit.>, IBM's TrueNorth <cit.>, BrainScaleS <cit.>, and ROLLS <cit.>, among others.
These are massively-parallel, low-power, analog, and/or digital systems (often a mix of the two) that are designed to simulate large-scale artificial neural networks rivalling the scale and complexity of real biological systems.
Despite growing excitement, we believe that methods to fully unlock the computational power of neuromorphic hardware are lacking. 
This is primarily due to a theoretical gap between our traditional, discrete, von Neumann-like understanding of conventional algorithms and the continuous spike-based signal processing of real brains that is often emulated in silicon <cit.>.

We use the term “neural compiler” to refer loosely to any systematic method of converting an algorithm, expressed in some high-level mathematical language, into synaptic connection weights between populations of spiking neurons <cit.>.
To fully leverage neuromorphic hardware for real-world applications, we require neural compilers that can account for the effects of spiking neuron models and mixed-analog-digital synapse models, and, perhaps more importantly, exploit these details in useful ways when possible.
There exist various approaches to neural engineering, including 
those by Denève et al. <cit.> and Memmesheimer et al. <cit.>.
However, the Neural Engineering Framework (NEF; <cit.>) stands apart in terms of software implementation (Nengo; <cit.>), large-scale cognitive modeling <cit.>, and neuromorphic applications <cit.>.


Competing methods consistently exclude important details, such as refractory periods and membrane voltage leaks, from their networks <cit.>.
The NEF, on the other hand, embraces biological complexity whenever it proves computationally useful <cit.> and/or improves contact with neuroscience literature <cit.>.
This approach to biological modeling directly mirrors a similar need to account for the details in neuromorphic hardware when building neural networks <cit.>.

Nevertheless, there are many open problems in optimizing the NEF for state-of-the-art neuromorphics.
In particular, we have been working to account for more detailed dynamics in heterogeneous models of the post-synaptic current (PSC) induced by each spike, as well as delays in spike propagation <cit.>.
The purpose of this report is to summarize our methods, both theoretical and practical, that have progressed in this direction.
There are similar challenges in extending the NEF to account for multi-compartment neuron models <cit.>, conductance-based synapses <cit.>, and to minimize the total number of spikes <cit.> – but these topics will not be addressed in this report.

The remainder of this report assumes that the reader is already comfortable with Nengo and the NEF in the context of engineering dynamical systems (i.e., “Principle III”). 
A technical yet accessible overview of the NEF can be found in <cit.>, and a tutorial for Nengo (version 2) can be found in <cit.>.



§ ACCOUNTING FOR SYNAPTIC DYNAMICS



This section provides a theoretical account of the effect of higher-order linear synapse models on the dynamics of the network, by summarizing the extensions from <cit.> and <cit.>.
This yields two novel proofs of Principle III from the NEF, and generalizes the principle to include more detailed synapses, including those modeling axonal transmission delays.



 §.§ Linear systems



Here we focus our attention on linear time-invariant (LTI) systems:

    ẋ⃗̇(t)    = Ax⃗(t) + Bu⃗(t) 
    y⃗(t)    = Cx⃗(t) + Du⃗(t)

where the time-varying vector x⃗(t) represents the system state, y⃗(t) the output, u⃗(t) the input, and the time-invariant “state-space” matrices (A,  B,  C,  D) fully determine the system's dynamics.
We will omit the variable t when not needed.

Principle III from the NEF states that in order to train the recurrent connection weights to implement (<ref>)—using a continuous-time lowpass filter  to model the PSC—we use Principle II to train the decoders for the recurrent transformation (τ A + I)x⃗, input transformation τ B u⃗, output transformation Cx⃗, and passthrough transformation Du⃗ <cit.>.
This drives the recurrent synapses with the signal τẋ⃗̇ + x⃗ so that their output is the signal x⃗, in effect transforming the synapses into perfect integrators with respect to ẋ⃗̇.
The vector x⃗ is then represented by the population of neurons via Principle I.
Thus, this provides a systematic approach for training a recurrent neural network to implement any linear dynamical system.
Now we show that this approach generalizes to other synaptic models.

For these purposes, the transfer function is a more useful description of the LTI system than (<ref>).  
The transfer function is defined as the ratio of Y(s) to U(s), given by the Laplace transforms of y⃗(t) and u⃗(t) respectively. 
The variable s denotes a complex value in the frequency domain, while t is non-negative in the time domain. 
The transfer function is related to (<ref>) by the following:

    F(s) = Y(s)/U(s) = C(sI - A)^-1B + D .


The transfer function F(s) can be converted into the state-space model (A,  B,  C,  D) if and only if it can be written as a proper ratio of finite polynomials in s. 
The ratio is proper when the degree of the numerator does not exceed that of the denominator. 
In this case, the output will not depend on future input, and so the system is `causal'. 
The order of the denominator corresponds to the dimensionality of x⃗, and therefore must be finite. 
Both of these conditions can be interpreted as physically realistic constraints where time may only progress forward, and neural resources are finite.



  
Continuous Identity

In order to account for the introduction of a synaptic filter h(t), we replace the integrator s^-1 in (<ref>) with H(s), where H(s) = ℒ{ h(t) }.
This new system has the transfer function C(H(s)^-1I - A)^-1B + D = F(H(s)^-1).
To compensate for this change in dynamics, we must invert the change-of-variables s ↔ H(s)^-1.
This means finding the required F^H(s) such that F^H(H(s)^-1) is equal to the desired transfer function, F(s).
We highlight this as the following identity:

    F^H ( 1/H(s)) = F(s) .

Then the state-space model (A^H,  B^H,  C^H,  D^H) satisfying (<ref>) with respect to F^H(s) will implement the desired dynamics (<ref>) given h(t).



  
Discrete Identity

For the discrete (i.e., digital synapse) case, we begin with F(z) and H(z) expressed as digital systems.
The form of H(z) is usually determined by the hardware, and F(z) is usually found by a zero-order hold (ZOH) discretization of F(s) using the simulation time-step (dt), resulting in the discrete LTI system:

    x⃗[t+dt]    = A̅x⃗[t] + B̅u⃗[t] 
    y⃗[t]    = C̅x⃗[t] + D̅u⃗[t] .

Here, we have the same relationship as (<ref>),

    F(z) = Y(z)/U(z) = C̅(zI - A̅)^-1B̅ + D̅.

Therefore, the previous discussion applies, and we must find an F^H(z) that satisfies:

    F^H ( 1/H(z)) = F(z) .

Then the state-space model (A̅^H, B̅^H, C̅^H, D̅^H) satisfying (<ref>) with respect to F^H(z) will implement the desired dynamics (<ref>) given h[t].
In either case, the general problem reduces to solving this change-of-variables problem for various synaptic models.
We now provide a number of results.
More detailed derivations are available in <cit.>.



  
Continuous Lowpass Synapse

Replacing the integrator s^-1 with the standard continuous-time lowpass filter, so that H(s) =  1/τ s + 1:

    F^H(τ s + 1) = F(s)  F^H(s) = C(sI - (τ A + I))^-1(τ B) + D

which rederives the standard form of Principle III from the NEF <cit.>.



  
Discrete Lowpass Synapse

Replacing the integrator with a discrete-time lowpass filter H(z) = 1 - a/z - a in the z-domain with time-step dt, where a = e^-dt/τ:

    F^H ( z-a/1-a) = F(z)  F^H(z) = C̅( zI - 1/1 - a(A̅ - aI) )^-1( 1/1-aB̅) + D̅.

Therefore, A̅^H = 1/1 - a(A̅ - aI), B̅^H = 1/1-aB̅, C̅^H = C̅, and D̅^H = D̅.
This mapping can dramatically improve the accuracy of Principle III in digital simulations (e.g., when using a desktop computer) <cit.>.






  
Delayed Continuous Lowpass Synapse

Replacing the integrator with a continuous lowpass filter containing a pure time-delay of length λ, so that H(s) = e^-λ s/τ s + 1:

    F^H ( τ s + 1/e^-λ s) = F(s)  F^H(s) = F ( 1/λ W_0(ds) - 1/τ) ,

where d = λ/τe^λ/τ and W_0(xe^x) = x is the principal branch of the Lambert-W function <cit.>.[This assumes that |η| < π and λRe[ s ] + λ/τ > - ηη, where η := λIm[ s ] <cit.>.]
This synapse model can be used to model axonal transmission time-delays due to the finite-velocity propagation of action potentials, or to model feedback delays within a broader control-theoretic context. 
To demonstrate the case where a pure time-delay of length θ is the desired transfer function (F(s) = e^-θ s), we let c = e^θ/τ and r = θ/λ to obtain the required transfer function:

    F^H(s) = c ( W_0(ds)/ds)^r = cr ∑_i=0^∞(i+r)^i-1/i! (-ds)^i .

We then numerically find the Padé approximants of the latter Taylor series.
More details and validation may be found in <cit.> and <cit.>.

To show the utility of these derivations, Fig. <ref> implements a bandpass filter, 
    F(s) = 1/1/w_0^2 s^2 + 1/w_0Q s + 1
 where w_0 = 2π· 20 rad/s is the angular frequency band and Q = 5 is a constant that is inversely proportional to the cutoff width. This compares the method of (<ref>), which is the standard version of Principle III used in NEF simulations, to the method of (<ref>), which is the discretized version of Principle III that takes into account the simulation time-step. This indicates that it is important to consider the discrete nature of the synapse when using digital hardware. We remark that the standard version is still correct for analog hardware (i.e., as dt → 0).






  
General Case

Finally, we consider any linear synapse model of the form:

    H(s) = 1/∑_i=0^k c_i s^i,

for some polynomial coefficients ( c_i ) of arbitrary degree k.
To the best of our knowledge, this class includes the majority of linear synapse models used in the literature.
For synapses containing a polynomial numerator with degree q > 1 (e.g., considering the Taylor series expansion of the box filter ϵ^-1 (1 - e^-ϵ s) s^-1), we take its [0, q]-Padé approximants to transform the synapse into this form within some radius of convergence.[This is equivalent to the approach taken in <cit.>.]
To map F(s) onto (<ref>), we begin by defining our solution to F^H(s) in the form of its state-space model:

    A^H    = ∑_i=0^k c_i A^i ,      C^H    = C ,
    
    B^H    = ( ∑_j=0^k-1 s^j ∑_i=j+1^k c_i A^i-j-1) B ,      D^H    = D ,

which we claim satisfies (<ref>).
To validate this claim, we assert the following algebraic relationship (proven in <cit.>):

    H(s)^-1I - A^H = ∑_i=0^k c_i ( s^i I - A^i ) = ( ∑_j=0^k-1 s^j ∑_i=j+1^k c_i A^i-j-1) (sI - A).

We now verify that (<ref>) satisfies (<ref>):

    F^H(H(s)^-1)    = C^H(H(s)^-1I - A^H)^-1 B^H + D^H 
       = C(H(s)^-1I - A^H)^-1( ∑_j=0^k-1 s^j ∑_i=j+1^k c_i A^i-j-1) B + D 
       = C(sI - A)^-1 B + D = F(s) .*

Since s^j is the j^th-order differential operator, this form of B^H states that we must supply the j^th-order input derivatives u⃗^(j), for all j = 1 … k - 1.
To be more precise, let us first define B^H_j:= ( ∑_i=j+1^k c_i A^i-j-1) B.
Then (<ref>) states that the ideal state-space model must implement the input transformation as a linear combination of input derivatives, ∑_j=0^k-1 B^H_j u⃗^(j).

However, if the required derivatives are not included in the neural representation, then it is natural to use a ZOH method by assuming u⃗^(j) = 0, for all j = 1 … k - 1:

    B^H = ( ∑_i=1^k c_i A^i-1) B ,

with A^H, C^H, and D^H as in (<ref>).
This is now an equivalent model to (<ref>) assuming ZOH, and in the form of the standard state-space model (<ref>).

The same derivation also applies to the discrete-time domain, with respect to the discrete synapse (corresponding to some implementation in digital hardware):

    H(z) = 1/∑_i=0^k c̅_i z^i.

Here, the only real difference (apart from notation) is the discrete version of (<ref>):

    B̅^H = ( ∑_j=0^k-1∑_i=j+1^k c̅_i A̅^i-j-1) B̅.

These mappings are made available by  and  in nengolib 0.4.0 <cit.>.
Additional details may again be found in <cit.>.



 §.§ Nonlinear systems



Here we derive two theorems for nonlinear systems, by taking a different perspective that is consistent with <ref>.
This generalizes the approach taken in <cit.>, which considered the special case of a pulse-extended (i.e., time-delayed) double-exponential.


We wish to implement some desired nonlinear dynamical system,

    ẋ⃗̇(t) = f(x⃗, u⃗) ,

using (<ref>) as the synaptic filter h(t).
Letting w⃗(t) = f^h(x⃗, u⃗) for some recurrent function f^h and observing that x⃗(t) = (w⃗∗ h)(t), we may express these dynamics in the Laplace domain:

    X⃗(s)/W⃗(s)   = 1/∑_i=0^k c_i s^i
          W⃗(s)    = X⃗(s) ∑_i=0^k c_i s^i = ∑_i=0^k c_i [ s^i X⃗(s) ] 
          w⃗(t)    = ∑_i=0^k c_i x⃗^(i)

since s is the differential operator. This proves the following theorem:
 
Let the function computed along the recurrent connection be:

    f^h(x⃗, u⃗) = ∑_i=0^k c_i x⃗^(i)

where x⃗^(i) denotes the i^th time-derivative of x⃗(t), and c_i are given by (<ref>). Then the resulting dynamical system is precisely (<ref>).


For the discrete case, we begin with some desired nonlinear dynamics expressed over discrete time-steps:

    x⃗[t+dt] = f̅(x⃗, u⃗) ,

using (<ref>) as the synaptic filter h[t], followed by an analogous theorem:
 
Let the function computed along the recurrent connection be:

    f̅^h(x⃗, u⃗) = ∑_i=0^k c̅_i x⃗^[i]

where x⃗^[i] denotes the i^th discrete forwards time-shift of x⃗, and c̅_i are given by (<ref>).
Then the resulting dynamical system is precisely (<ref>).

The proof for the discrete case is nearly identical.
For sake of completeness, let w⃗[t] = f̅^h(x⃗, u⃗) for some recurrent function f̅^h and observe that x⃗[t] = (w⃗∗ h)[t]:

    X⃗(z)/W⃗(z)   = 1/∑_i=0^k c̅_i z^i
          W⃗(z)    = X⃗(z) ∑_i=0^k c̅_i z^i = ∑_i=0^k c̅_i [ z^i X⃗(z) ] 
          w⃗[t]    = ∑_i=0^k c̅_i x⃗^[i]

since z is the forwards time-shift operator. 



  
Continuous Lowpass Synapse

For standard Principle III, we have H(s) = 1/τ s + 1  k = 1, c_0 = 1 and c_1 = τ, 

    f^h(x⃗, u⃗) = c_0 x⃗^(0) + c_1 x⃗^(1) = x⃗ + τẋ⃗̇ = τ f(x⃗, u⃗) + x⃗.

Note that (<ref>) is consistent with (<ref>) and with Principle III from the NEF.



  
Discrete Lowpass Synapse

For the discrete case of Principle III, we have H(z) = 1 - a/z - a, where a = e^-dt / τ  k = 1, c̅_0 = -a(1 - a)^-1, c̅_1 = (1 - a)^-1,

    f̅^h(x⃗, u⃗) = c̅_0 x⃗^[0] + c̅_1 x⃗^[1] =  (1 - a)^-1(f̅(x⃗, u⃗) - ax⃗) .

Note that (<ref>) is consistent with (<ref>).



  
Continuous Double Exponential Synapse

For the double exponential synapse:

    H(s) = 1/(τ_1 s + 1)(τ_2 s + 1) = 1/τ_1 τ_2 s^2 + (τ_1 + τ_2)s + 1


    f^h(x⃗, u⃗)    = x⃗ + (τ_1 + τ_2) ẋ⃗̇ + τ_1 τ_2 ẍ⃗̈
       = x⃗ + (τ_1 + τ_2) f(x⃗, u⃗) + τ_1 τ_2 ( ∂ f(x⃗, u⃗)/∂x⃗· f(x⃗, u⃗) + ∂ f(x⃗, u⃗)/∂u⃗·u̇⃗̇) .

In the linear case, this simplifies to:

    f^h(x⃗, u⃗) = ( τ_1 τ_2 A^2 + (τ_1 + τ_2) A + I ) x⃗ + ( τ_1 + τ_2 + τ_1 τ_2 A ) B u⃗ + τ_1 τ_2 B u̇⃗̇.




  
Linear Systems

As in <ref>, Theorems <ref> and <ref> require that we differentiate the desired dynamical system.
For the case of nonlinear systems, this means determining the (possibly higher-order) Jacobian(s) of f, as shown in (<ref>).
For the special case of LTI systems, we can determine this analytically to obtain a closed-form expression.
By induction it can be shown that:

    x⃗^(i) = A^ix⃗ + ∑_j=0^i-1 A^i-j-1 B u⃗^(j).

Then by expanding and rewriting the summations:

    f^h(x⃗, u⃗)    = ∑_i=0^k c_i x⃗^(i)
       = ∑_i=0^k c_i [ A^i x⃗ + ∑_j=0^i-1 A^i-j-1 B u⃗^(j)] 
       = ( ∑_i=0^k c_i A^i )_4emRecurrent
    Matrixx⃗ + ∑_j=0^k-1( ∑_i=j+1^k c_i A^i-j-1)B_7emInput Matricesu⃗^(j).

The discrete case is identical:

    f̅^h(x⃗, u⃗) = ( ∑_i=0^k c̅_i A̅^i )_4emRecurrent
    Matrixx⃗ + ∑_j=0^k-1( ∑_i=j+1^k c̅_i A̅^i-j-1)B̅_7emInput Matricesu⃗^[j].

This gives a matrix form for any LTI system with a k^th order synapse, provided we can determine u⃗^(j) or u⃗^[j] for 0 ≤ j ≤ k-1.





Again, (<ref>) is consistent with (<ref>) and (<ref>), as is (<ref>) with (<ref>).





§ ACCOUNTING FOR SYNAPTIC HETEROGENEITY



We now show how <ref> can be applied to train efficient networks where the i^th neuron has a distinct synaptic filter h_i(t), given by:

    H_i(s) = 1/∑_j=0^k_i c_ij s^j.

This network architecture can be modeled in Nengo using nengolib 0.4.0 <cit.>. 
This is particularly useful for applications to neuromorphic hardware, where transistor mismatch can change the effective time-constant(s) of each synapse.
To this end, we abstract the approach taken in <cit.>.
We show this specifically for Theorem <ref>, but this naturally applies to all methods in this report.

Recalling the intuition behind Principle III, our approach is to separately drive each synapse h_i with the required signal f^h_i(x⃗, u⃗) such that each PSC becomes the desired representation x⃗.

Thus, the connection weights to the i^th neuron should be determined by solving the decoder optimization problems for f^h_i(x⃗, u⃗) using the methods of <ref> with respect to the synapse model h_i. 
This can be repeated for each synapse to obtain a full set of connection weights.
While correct in theory, this approach displays two shortcomings in practice: (1) we must solve n optimization problems, where n is the number of post-synaptic neurons, and (2) there are 𝒪(n^2) weights, which eliminates the space and time efficiency of using factorized weight matrices <cit.>.

We can solve both issues simultaneously by taking advantage of the linear structure within f^h_i that is shared between all h_i.
Considering Theorem <ref>, we need to drive the i^th synapse with the function:

    f^h_i(x⃗, u⃗) = ∑_j=0^k_i c_ijx⃗^(j).

Let d⃗^j be the set of decoders optimized to approximate x⃗^(j), for all j = 0 … k, where k = max_i k_i.
By linearity, the optimal decoders used to represent each f^h_i(x⃗, u⃗) may be decomposed as:

    d⃗^f^h_i = ∑_j=0^k_i c_ijd⃗^j .

Next, we express our estimate of each variable x⃗^(j) using the same activity vector a⃗:

    x⃗̂⃗^(j) = a⃗d⃗^j.

Now, putting this all together, we obtain:

    a⃗d⃗^f^h_i = ∑_j=0^k_i c_ija⃗d⃗^j = ∑_j=0^k_i c_ijx⃗̂⃗^(j)≈ f^h_i(x⃗, u⃗) .


Therefore, we only need to solve k+1 optimization problems, decode the “matrix representation” [ x⃗̂⃗^(0), x⃗̂⃗^(1), …, x⃗̂⃗^(k)], and then linearly combine these k+1 different decodings as shown in (<ref>)—using the matrix of coefficients c_ij—to determine the input to each synapse.
This approach reclaims the advantages of using factorized connection weight matrices, at the expense of a factor 𝒪(k) increase in space and time efficiency.



§ DISCUSSION


We have reviewed three major extensions to the NEF that appear in recent publications.
These methods can be used to implement linear and nonlinear systems in spiking neurons recurrently coupled with heterogeneous higher-order mixed-analog-digital synapses.
This provides us with the ability to implement NEF networks in state-of-the-art neuromorphics while accounting for, and sometimes even exploiting, their nonideal nature.

While the linear and nonlinear methods can both be used to harness pure spike time-delays (due to axonal transmission) by modeling them in the synapse, the linear approach provides greater flexibility.
Both extensions can first transform the time-delay into the standard form of (<ref>) via Padé approximants, which maintains the same internal representation x⃗ as the desired dynamics (within some radius of convergence).
But the linear extension also allows the representation to change, since it is only concerned with maintaining the overall input-output transfer function relation.
In particular, we derived an analytic solution using the Lambert-W function, which allows the neural representation x⃗, and even its dimensionality, to change according to the expansion of some Taylor series.
The linear case is also much simpler to analyze in terms of the network-level transfer function that results from substituting one synapse model for another.
For all other results that we have shown, the linear extension is consistent with the nonlinear extension, as they both maintain the desired representation by fully accounting for the dynamics in the synapse.



§ ACKNOWLEDGEMENTS


We thank Wilten Nicola for inspiring our derivation in <ref> with his own phase-space derivation of Principle III using double exponential synapses for autonomous systems (unpublished).
We also thank Kwabena Boahen and Terrence C. Stewart for providing the idea used in <ref> to separately drive each h_i, and for improving this report through many helpful discussions.



ieeetr






§ APPENDIX




 §.§ The Neural Engineering Framework



For convenience we provide a self-contained treatment of the Neural Engineering Framework (NEF) first proposed by Eliasmith and Anderson <cit.>. A simplified technical summary can be found in <cit.>.

The NEF is a set of three mathematical principles for describing neural computation.
This has the aim of relating low-level aspects of biology to high-level cognition.
Practically speaking, it is a theoretical framework that may be used to implement algorithms using neural components as the primitive operations.
This latter perspective is most relevant when engineering systems using neuromorphic hardware.


The NEF is most commonly applied to building dynamic spiking neural networks, but also applies to non-spiking and feed-forward networks. 
Recently, the NEF has been used to build what remains the world's largest functioning brain model <cit.>, capable of performing perceptual, motor, and cognitive tasks.
This model relies fundamentally on the principles described in part below.



  §.§.§ Principle I – Representation


We use x⃗(t) ∈ℝ^k to denote a continuous time-varying vector that is represented by a population of neurons.  
To describe this representation, we define a nonlinear (spiking or non-spiking) encoding and a linear decoding which together determine how neural activity relates to the represented vector.

Specifically, we choose encoders E ∈ℝ^n × k, gains α_i > 0, and biases β_i as parameters for the encoding, which map x⃗(t) to neural activities.
These parameters can be chosen based on prior knowledge (e.g., tuning curves, firing rates, sparsity, etc.), or randomly.
Most often, we select distributions from which to sample these parameters, incorporating known constraints on the system under consideration.
For a spiking neuron, the encoding is:

    δ_i^x⃗(t) = G_i [ α_i ⟨ E_i, x⃗(t') ⟩ + β_i : 0 ≤ t' ≤ t ]

where δ_i^x⃗ is the neural spike-train generated by the i^th neuron and G_i[·] is the model for a single neuron (e.g., a leaky integrate-and-fire (LIF) neuron, a conductance neuron, etc.).
We then analytically determine the firing rates of each neuron under constant inputs: 

    r_i(v⃗) = 𝔼_t ≥ 0[ δ_i^x⃗  ;  x⃗(t) = v⃗]


In the case of non-spiking neurons, we use its rate model in place of (<ref>). 

In biological systems, any transfer of information between layers of a network are affected by a post-synaptic response.  
Consequently we introduce a post-synaptic filter h(t), referred to as the synapse model, that accounts for the temporal dynamics of a receiving neuron's synapse.

To determine the vector represented by the activities, we define a linear decoding: 

    (x̂⃗̂∗ h)(t) = ∑_i=1^n (δ_i^x⃗∗ h)(t) D_i


To complete our characterization of a representation, we need to determine the linear decoders D_i.  
This determination is the same for the first and second principles, so we consider it next.



  §.§.§ Principle II – Transformation


The second principle of the NEF addresses the issue of computing nonlinear functions of the vector space.
The encoding remains as defined in (<ref>).
However, we may decode nonlinear functions of x⃗(t) by the same linear combination of temporal basis functions:

    (f̂(x⃗) ∗ h)(t) = ∑_i=1^n (δ_i^x⃗∗ h)(t) D^f_i
 

In general, we solve for the decoders D^f ∈ℝ^n × l to compute a desired nonlinear function f : S →ℝ^l, where S = {x⃗(t) : t ≥ 0} is the domain of the signal –
typically the unit k-ball {v⃗∈ℝ^k : v⃗ = 1} or the k-cube [-1, 1]^k.
To account for fluctuations introduced by neural spiking, or other sources of uncertainty, we introduce a noise term η∼𝒩(0, σ^2) to regularize the solution:

    D^f = argmin_D ∈ℝ^n, l∫_S f(v⃗) - ∑_i=1^n (r_i(v⃗) + η) D_i ^2   d^kv⃗


Note this training only depends on r_i(v⃗) for v⃗∈ S.
It does not depend on x⃗(t).
Performing this regularized least-squares optimization provides the decoders needed by equations (<ref>) and (<ref>) to represent the vector or compute transformations, respectively.
We can use singular value decomposition (SVD) to characterize the class of functions that can be accurately approximated by the chosen basis functions [: pp. 185–217].






The accuracy of this approach for spiking neurons relies on r_i(v⃗) being a suitable proxy for δ_i^x⃗ whenever x⃗(t) = v⃗.
This is clearly true for constant x⃗, and is also a suitable model in practice for low-frequency x⃗(t).


In general, we can use (<ref>) and (<ref>) to derive the connection weight matrix between layers that computes the desired function f. 

Specifically, the weight matrix 𝒲∈ℝ^m, n which decodes from neuron i and encodes into neuron j is given by:

    𝒲_ji = ⟨ E_j, D^f_i ⟩ = (E(D^f)^T)_ji


That is, the matrices D and E are simply a low-dimensional factorization of 𝒲.
This observation implies that the process of decoding and then encoding is equivalent to the more familiar expression for mapping activities between layers of neurons in a conventional neural network: 


    G [ α⃗𝒲δ⃗^x⃗(t) + β⃗]


The crucial difference is that a constant dimensionality k means O(n) operations per time-step to implement (<ref>) as compared to O(n^2) operations for typical weight matrices.

In essence, structure is imposed by the choice of representing x⃗(t) in neural activities, which in turn becomes latent within the factored weight matrix 𝒲.
This additional structure gives us a way to constrain the model by describing its state, and therefore understand in detail what the network is computing.



  §.§.§ Principle III – Dynamics


Given our ability to compute nonlinear functions with Principle II, the methods described here apply to both linear and nonlinear dynamical systems.  
However, for simplicity of exposition, we focus our discussion on linear time-invariant (LTI) systems:

    ẋ⃗̇(t)    = Ax⃗(t) + Bu⃗(t) 
    y⃗(t)    = Cx⃗(t) + Du⃗(t)

where the time-varying vector x⃗(t) represents the system state, y⃗(t) the output, u⃗(t) the input, and the time-invariant matrices (A, B, C, D) fully define the system dynamics. 

Principle III from the NEF states that in order to train the recurrent connection weights to implement (<ref>) using a continuous-time lowpass filter h(t) = 1/τ e^-t/τ, we use Principle II to train the decoders for:


  * the recurrent transformation f(x⃗(t)) = (τ A + I)x⃗(t);

  * the input transformation f(x⃗(t)) = τ B x⃗(t);

  * the output transformation f(x⃗(t)) = Cx⃗(t); and

  * the passthrough transformation f(x⃗(t) ) = Dx⃗(t) [; pp. ].




Fig. <ref> depicts this mapping in the Laplace domain. In particular, X(s) is the Laplace transform of the vector represented in neurons via Principle I.
The functions on each connection implement the desired matrix transformations as described above using Principle II.
In <ref> we extend this mapping to other models of the synaptic filter H(s), and provide two novel proofs of Principle III.


