Actually I think there's something slightly more serious/general going on, with how the discretization has been incorrectly formulated to be preemptive with the state at the cost of effectively predicting the future passthrough. This basically means that transfer functions that are not strictly proper are not possibly going to be able to behave correctly (and fortunately Lowpass/Alpha are strictly proper, so they have no passthrough). So, I need to be more sure of what this means should happen in Nengo, so I'm closing this issue for now. Sorry for the spam.

I think I've found the root of the problem. If we put an infinitesimally small passthrough on the system (which keeps the system's intended behaviour identical), then suddenly things change because Nengo inserts a delay.

sys1 = nengo.LinearFilter([0, 0, 1], [0.0001, 0.02, 1])
sys2 = nengo.LinearFilter([1e-16, 0, 1], [0.0001, 0.02, 1])

The blue line in the second graph is Nengo, which is inconsistent from the first simulation. The green lines are from lfilter (and they match for both systems). Both simulations use the General step method, and both have identical num/den after discretization within epsilon factor.

I believe the issue revolves around this line num = num[1:] if num[0] == 0 else num in Nengo, which I find worrisome. This basically says: "if there's no passthrough, then simulate one step ahead (to avoid delay). Otherwise, if there is a passthrough, then simulate it normally with a delay on the output".

Changing == 0 to np.allclose would fix this one particular case, but still give inconsistent behaviour between proper and strictly proper transfer functions. But maybe that's the price we need to pay if we want our synapses to return future (undelayed) output.

I think this is fixed in #1535. Here's the plots I get for the first systems:

Here's the plots for the second systems:

Looks fixed to me!

If we use the example from #1535 (comment) or sys2 from the first post here, with a different input signal, you can see there is still a difference (now plotting the error):

For nengo.LinearFilter([-1, 0, 1e6], [1, 2, 1]):

For nengo.LinearFilter([0.0001, 0, 0], [0.0001, 0.02, 1]):

import pylab
import numpy as np
from scipy.signal import lfilter, cont2discrete, tf2ss

import nengo

sys1 = nengo.LinearFilter([-1, 0, 1e6], [1, 2, 1])
sys2 = nengo.LinearFilter([0.0001, 0, 0], [0.0001, 0.02, 1])

for sys in (sys1, sys2):

    with nengo.Network() as model:
        x = nengo.Node(output=nengo.processes.PresentInput([1000, 0], 0.001))
        p = nengo.Probe(x, synapse=sys)
        p_x = nengo.Probe(x, synapse=None)

    sim = nengo.Simulator(model)
    sim.run(0.02)

    (num,), den, _ = cont2discrete((sys.num, sys.den), sim.dt)
    expected = lfilter(num, den, sim.data[p_x].squeeze())

    pylab.figure()
    pylab.title("len(num)=%d, len(den)=%d" % (len(sys.num), len(sys.den)))
    pylab.plot(sim.trange(), sim.data[p].squeeze() - expected.squeeze())
    #pylab.plot(sim.trange(), expected)
    pylab.xlabel("Time (s)")
    pylab.ylabel("Error")
    pylab.show()

I believe the reason is that the discrete state-space model for an LTI is defined by equation (5.5):

x[k+1] = Ax[k] + Bu[k]
y[k] = Cx[k] + Du[k]

which is also what nengolib does. However, #1535 uses the equations:

x[k+1] = Ax[k] + Bu[k]
y[k] = Cx[k+1] + Du[k]

Note the difference is in Cx[k] versus Cx[k+1].

I believe the reason we are doing that is because there is already a delay elsewhere in Nengo. Prior to #1535 (and in nengolib) this delay was compensated for by applying the z-operator to the transfer function (which you can do if and only if it is strictly proper -- that is, no passthrough). This also meant that when D != 0 on master there would be two time-step delays.

In #1535, this is instead compensated for by using the most recent value of x in the LTI equations, which, for basic cases such as lowpass filters, effectively undoes that delay and ends up doing the correct thing overall. However, in the event that there is a passthrough (D != 0) then this is no longer correct, because it is using the wrong value of u[.] across the two lines in the system. Unfortunately the example stimulus I chose in this issue was not very good, because cos is slow enough on this time-scale that it's hard to spot the difference, that is u[k] ~ [k+1] relative to the size of the passthrough. The stimulus in this post alternates each time-step which makes the difference much more clear.

A correct fix would be to use Du[k+1] in the second line of the system, but of course this requires knowing a future input value.

In the context of #1535 I believe this means there won't be an ideal fix unless there is an easy way to eliminate the delay elsewhere in the graph, and move it explicitly into the synapse. Circular dependencies would occur if and only if all of the synapses along the chain have D != 0. Note that synapse=None would reduce to (A, B, C, D) = (0, 0, 0, 1) and correspondingly there would no longer be a distinction between synapse=None and synapse=0, since neither would produce a delay. To recover the old behaviour for synapse=0 (if a pure time-step delay is explicitly wanted on a connection) one could do so by setting synapse=~z where z is the discrete time-shift operator, which would be realized here as synapse=nengo.LinearFilter([1], [1, 0], analog=False). After struggling with this for the past 3 years I'm pretty confident that this is the right thing to be doing, but I'll spend a little time trying to see how easy this fix is.

I spent a few hours and figured out the proper way to implement synapses when D != 0 in the reference backend. The current delay in synapses (e.g., for synapse=0) comes from the fact that SimProcess uses mode='update' which happens at the very end of the time-step (after any other operators are able to 'read' the value of the output signal).

To follow my recommendation of "moving the delay explicitly into the synapse" we need to split up the SimProcess operator into a number of separate operators. One group of operators that uses mode='set' to set y = Cx + Du at the start of the time-step, and another that uses mode='update' to update x = Ax + Bu at the end of the time-step. I think when it's done this way the operator dependency graph properly figures out the rest. That is it will first set u at the start of the time-step, then set y, and then finally update x at the end.

The biggest issue with this approach is it requires special builder logic that only makes sense for linear transfer functions, but it makes a difference in behaviour (compared to #1535) by correcting the case of D != 0.

One possible advantage of this change (aside from correcting the behaviour) is it would give the operator optimizer the ability to group together all of the synapse multiplies (because it would be expressed in terms of only basic operators). I would expect this to give performance gains when there are lots of synapses between small ensembles.

PR #1536 provides a proof-of-concept implementation of the above suggestion.

One caveat with this is that many backends will not support this, particularly the neuromorphic ones. That said, e.g. Loihi doesn't support any synapses that can't be made out of exponential filters, so D != 0 is a non-issue there.

My impression is that D != 0 is pretty rare with biological synapse models (I'm not aware of any, but I've never used any beyond exponential, double exponential and alpha). I don't think it's worth adding significant complexity to the builder if we don't have some explicit use-cases (and even then, maybe not worth it if they're edge cases). We intentionally chose to have synapses be an "update" op because it allows feedback connections and such to work in a variety of cases (e.g. even if you have a feedback connection on a node). It would be a bit strange if this worked for some synapses but not others.

One caveat with this is that many backends will not support this, particularly the neuromorphic ones. That said, e.g. Loihi doesn't support any synapses that can't be made out of exponential filters, so D != 0 is a non-issue there.

Other backends (e.g., nengo-ocl) could keep doing the same thing that Nengo's backend did prior to #1535, and that was to have a two time-step delay per synapse whenever D != 0 and len(A) > 0. This implements the transfer function correctly, but with an extra system-level delay. Whatever the choice may be, I'd think of #1536 as being an improvement to the Nengo reference backend, and it's up to other backends to decide whether to explicitly support something or to do so with some caveats (e.g., an extra delay). At the end of the day it's inevitable that there will be differences between backends, especially with delays and synaptic discretization. The discussions / examples here can facilitate documenting them in a way that makes such caveats clear to backend designers and users alike.

My impression is that D != 0 is pretty rare with biological synapse models (I'm not aware of any, but I've never used any beyond exponential, double exponential and alpha).

One possibility is in modelling gap junctions (electrical synapses) when dt is to be kept relatively high (e.g., 1ms).

I don't think it's worth adding significant complexity to the builder if we don't have some explicit use-cases (and even then, maybe not worth it if they're edge cases).

Here are two quick neuromorphic use-cases:

1. This paper that extends the NEF for Braindrop synapses used a digital synapse:

   du = nengo.LinearFilter([1, -1], [1, 0], analog=False)
   

   (i.e., 1 - 1/z) with the Nengo reference simulator, in order to simulate a discrete-time differentiator within the network. This was done on the CPU in order to validate certain math / experiments. We have been finding Nengo quite useful for these use-cases (as a general-purpose dynamical systems solver that can swap neurons in and out for various state-variables).

2. The same also applies to synapses such as nengolib.synapses.BoxFilter which has D != 0 and is used to model the DAC (spike->analog pulse-extender) on Braindrop. IIRC the spike gets through on the time-scale of microseconds and is extended to a width of about half a millisecond, so we can accurately model it using a pure passthrough component on a time-scale of dt=1e-4 seconds.

However, neither transfer function does the correct thing with #1535, and produces an extra delay on master. #1536 solves both problems.

We intentionally chose to have synapses be an "update" op because it allows feedback connections and such to work in a variety of cases (e.g. even if you have a feedback connection on a node). It would be a bit strange if this worked for some synapses but not others.

So in #1536, the situations where a recurrent synapse would work would be D = 0. If D != 0 then you could not have a recurrent synapse. But this was already the case for synapse = None which is equivalent to (A, B, C, D) = (0, 0, 0, 1) in #1536 (i.e., assuming no extra delay is added). In other words, previously users had to be aware that setting synapse=None was the same as instantaneously passing the input through and thus could not have such a recurrent feedback loop. If they attempted to put that into the synapse instead, they would get an extra delay. Instead, #1536 has no unnecessary delays and makes it explicit that there is a circular dependency in the opgraph if and only if all synapses along some feedback loop contain a passthrough. This should come as no surprise to users versed in control theory, which is somewhat of a prerequisite to be designing synapses where D != 0.

In my opinion it's much stranger to invent a different definition of an LTI system that does not implement the discrete F(z) that the user asked for, but instead subtracts a discrete-time differentiator scaled by D:

as #1535 does (note that this achieves equality if and only if D = 0). If #1535 is the long-term fix and we're taking the view that D != 0 is an edge-case, then it would make more sense to emit a warning indicating that the passthrough is applied to the input at the wrong time-step or to drop support for D != 0 entirely.

Here's another quick example to demonstrate some of the weirdness of implementing a different transfer function from what was requested. This applies to both master and #1535, but it is different than the other examples because it does not invoke the General stepper. However, it falls into the above discussion because (num, den) = (1, 1) becomes (A, B, C, D) = ((), (), (), 1) which is handled as a special case by a different Step class.

import nengo

sys1 = nengo.LinearFilter([1], [1, 0], analog=False)
sys2 = nengo.LinearFilter([1], [1], analog=False)

for sys in (sys1, sys2):
    with nengo.Network() as model:
        x = nengo.Node(output=nengo.processes.PresentInput([1000, 0], 0.001))
        p = nengo.Probe(x, synapse=sys)
        p_x = nengo.Probe(x, synapse=None)

    with nengo.Simulator(model, progress_bar=False) as sim:
        sim.run_steps(4, progress_bar=False)
        
    print(sim.data[p].squeeze(), sys.filt(sim.data[p_x].squeeze(), y0=0))

[   0. 1000.    0. 1000.] [1000.    0. 1000.    0.]
[   0. 1000.    0. 1000.] [1000.    0. 1000.    0.]

This example demonstrates two other forms of odd behaviour that stem from the above discussion:

1. Simulator is one time-step behind filt. (I would expect both arrays to be the same on a given line.)
2. (num, den) = ((1,), (1, 0)) is equivalent to (num, den) = (1, 1), even though the former is 1/z (a single time-step delay) while the latter is an identity pasthrough. (I would expect the first line to be a delayed version of the second line.)

#1536 solves both issues (tests have been modified and added corresponding to both points).

My impression is that D != 0 is pretty rare with biological synapse models (I'm not aware of any, but I've never used any beyond exponential, double exponential and alpha). I don't think it's worth adding significant complexity to the builder if we don't have some explicit use-cases (and even then, maybe not worth it if they're edge cases).

I forgot to mention perhaps the earliest use-case for D != 0, given by Tripp and Eliasmith (2010):

A system acting as a practical differentiator should have a transfer function similar to y/u = s/(τ s + 1), where y is the output of the system and τ is a small constant. (Note that this is the derivative of y/u = 1/(τ s + 1), which is the Laplace transform of the nonamplifying impulse response (1/τ )e−t/τ , an exponential function with integral one.)

This paper argues that the analog transfer function nengo.LinearFilter([1, 0], [tau, 1]) can be achieved in a number of biologically plausible ways, thus situating it as a good model for temporal differentation. This model is normalized and generalized by nengolib.synapses.Highpass. (And to reiterate: on master this has an extra delay, #1535 gives incorrect behaviour, and #1536 produces the ideal transfer function.)