Compute-in-memory implementation of state space models for event sequence processing

Domain: physics
Published: 2026-01-09
Editorial outcome: major_revision

Abstract
State space models have recently emerged as a powerful framework for long sequence processing, outperforming traditional methods on diverse benchmarks. Fundamentally, state space models can generalize both recurrent and convolutional networks and have been shown to even capture key functions of biological systems. Here we report an approach to implement SSMs in energy-efficient compute-in-memory hardware to achieve real-time, event-driven processing. Our work re-parameterizes the model to function with real-valued coefficients and shared decay constants, reducing the complexity of model mapping onto practical hardware systems. By leveraging device dynamics and diagonalized state transition parameters, the state evolution can be natively implemented in crossbar-based compute-in-memory systems combined with memristors exhibiting short-term memory effects. Through this algorithm and hardware co-design, we show the proposed system offers both high accuracy and high energy efficiency while supporting fully asynchronous processing for event-based vision and audio tasks.

Introduction
Sequencemodelingisafundamentalcomputationaltaskessentialfor diversefields,fromnaturallanguageprocessing1totime-seriesdata analysis2. A significant focus within the area of neuromorphic com- puting lies on sequential information processing, particularly effi- cientlyprocessingasynchronous,event-baseddatastreams.Suchdata can originate directly from sensors such as event cameras, which capture changes in luminance and generate sparse events3,4, or be derivedbyconvertingtraditionaldatasetsintospike-basedformats5to takeadvantageofsparsityincommunicationandprocessing.

SpikingNeuralNetworks(SNNs)aregenerallythemodelofchoice foreventdataprocessing,astheirfundamentaloperationisbasedon asynchronous spike events, mirroring the data format6. The asyn- chronous, event-driven nature holds considerable potential for energy-efficientdataprocessing,ascomputationisideallyperformed onlywhennewdataarrives7,8.However,despitetheirpotential,train- ingdeepSNNscanbechallengingduetothenon-differentiablenature ofspikegeneration,whichnecessitatestheuseofsurrogategradient techniquesforbackpropagation9.Whilethesetechniqueshaveproven increasinglyrobust10,11,theinherentbinarizationofactivationsmaystill limit model performance, a constraint recently identified as a key factor for scaling up SNNs12. Tradeoffs between different encoding methods and neuron modelshave also been found to affect model performance13.

Toimprovemodelaccuracy,SNNsoftenincorporatespatialpro- cessing blocks such as convolution or transformer blocks to help extractspatialfeatures14–16.However,theuseoftheseblocksrequires collectingasequenceofspikesandconvertingthesequenceofasyn- chronouseventsintoaframe-basedrepresentation,therebysacrificing theinherentsparsityandlow-latencyadvantagesoftheevent-based approach. This mismatch between asynchronous data inputs and synchronous processing algorithms motivates the search for more suitablemodels.

Furthermore, from a hardware perspective, convolution and transformer blocks used in the SNN models are not well suited for energy efficient architectures suchascompute-in-memory(CIM).In CIM,weightsarestoredstationarilyinhigh-densitymemorycellssuch asResistiveRAM(RRAM)17–19orPhase-ChangeMemory20,21,andcom- putations,primarilyvector-matrixmultiplication(VMM),aredirectly performedinsidethememoryarraythroughOhm’slawandKirchhoff’s current law22,23. This in-situ computation can dramatically minimize data movement compared to traditional von Neumann designs24–26. However,forconvolutionoperations,thesameweightsarere-usedfor many input patches thus the weight read costs can be efficiently amortizedinwell-designeddigitalsystems,diminishingthebenefitsof CIM and amplifying the challenges of the slower read in emerging devices used for CIM implementations. The need to load different inputpatchesfromthefeaturemapmemorytoperformconvolutions also requires data re-shuffling, again diminishing CIM benefits. As a result, convolution operations are generally not well matched with CIM implementations. In Transformers, the attention mechanism involvesthemultiplicationsoftheQuery,Key,andValuematrices,all ofwhichdependontheinputsotheattentionoperationisnotfriendly to stationary CIM arrays which assume one of the matrices (i.e., weights)isstatic.Consequently,practicalimplementationsofTrans- formersadopthybriddesigns:theyleveragestationaryCIMonlyfor VMM involving static weights (e.g., in Fully Connected Layers), but delegatetheinput-dependentattentionscorecomputationsortheKV cachemanagementtoseparatedigitalcircuitry27.

Toovercomeboththearchitecturalmismatchwithconvolutions andthememoryandscalingbottleneckinherentinTransformers,in thiswork,wepresentaStateSpaceModel(SSM)basedapproachthat isnativelycompatiblewithCIMhardware.Weshowthattheproposed SSMallowsasynchronousspikesequenceprocessingwithhighmodel performancewithouttheneedforexplicitspatialprocessingblocks (e.g., convolutions). SSMs have recently gained prominence for sequencemodelingespeciallyontasksrequiringthecaptureoflong- range dependencies28,29, offering strong performancecompatible or exceeding that of recurrent, convolutional, and Transformer-based methods.Grounded in control theory with theoreticalguarantee to stably approximate continuous-time convolutional kernels, SSMs modelthesequencethroughalatentstatevectorthatevolvesover timebasedonthecurrentinputandthepreviousstate30.Akeyinno- vationinrecentSSMslikeS431andS532istheuseofadiagonalizedstate transitionmatrix,whichallowsforhighlyparallelizablecomputationof thestateupdatesandefficientmodelingoflong-termdependencies throughexponentialdecaydynamics.S4setsnewstate-of-the-arton theLong-RangeArena—solvingthe16384-lengthPath-Xtaskwith91% accuracyandrunning60×fasterthanTransformer-basedmodels.Its successorS5matchesS4’sefficiencywhileaveraging87.4%onLRAand 98.5 % on Path-X. Subsequent developments such as Mamba33 and Mamba234employselectiveSSMlayerstoachievelinear-timescaling tomillion-lengthcontexts.

Beyondsuperiormodelperformance,SSMshavebeenshownto generalize recurrent, convolution, and attention-based models34,35, offering a unified theoreticalframework for sequentialdata proces- sing. SSMs do not rely on sliding-window operations such as con- volutionsorTransformerblocks,suggestingtheyarecompatiblewith asynchronousspikeprocessing36andalsopotentiallymorefriendlyto CIMhardware.However,mappingSSMsontopracticalCIMhardware presentsitsownchallenges.Firstly,manyhigh-performingSSMs(e.g., S4, S5) rely on complex-valued state representations and computa- tions.Emulatingcomplexarithmeticusingreal-valuedcircuitssuchas CIM introduces overhead in terms of complexity, area, or power. Secondly,SSMsinherentlyrelyonmaintainingandupdatingahidden statethatcapturestemporalcontext,ofteninvolvingspecificdynam- ics.WhileCIMexcelsatVMM,relyingonexternaldigitalcomponents (e.g., registers,digital logic) to storeand update the state canrein- troducesignificantoverheadbetweentheanalogCIMcoreanddigital peripheral,negatingenergyandlatencybenefits.

In this study we address these challenges through a co-design approach,bydevelopingbothahardware-friendlySSMalgorithmand ahardwarearchitectureforefficientSSMmappingandexecutionon event-baseddata.WhilerecentworkshavesuccessfullymappedSSMs ontodigitalneuromorphicprocessors37,38,ourapproachrepresentsa co-design methodology leveraging analog CIM to perform key operations, including state update, directly through device physics with potentially greater efficiency and scalability. We first re- parameterizetheSSMbyrestrictingthestatetransitiondynamicsto operatepurelywithreal-valuedcoefficients,eliminatingtheoverhead associatedwithcomplexarithmeticinhardware.Furthermore,we constrain the exponential decay parameter (λ) associated with the statedynamicstoone(orafew)fixedvalue(s)perblock,significantly simplifyinghardwareimplementation.Subsequently,weco-designthe hardware architecture to natively implement all major operations, includingthestatedynamics,inanalogdomainthroughphysics.The hardwarearchitectureconsistsofmodular,asynchronousblocksthat matchthefunctionsofSSMblocks.Withineachhardwareblockare RRAM crossbararrays thatperform the necessaryVMM operations, and a group of memristors exhibiting short-term memory (STM) effects (e.g.,TungstenOxide(WO ) based memristors39)tonatively x store the state components and implement the diagonalized state updatefunctions.Thephysicaldecayrates(λ)oftheseSTMdevicesare tunedduringfabricationtoalignwiththefewfixedvaluesrequiredby the modified SSM algorithm.Overall,our co-designed system using real-valuedSSMwithfixeddecaysconsistentlymatchesorexceedsthe performanceofleadingmodelsonbothaudioandvisionevent-stream benchmarkswithsmallnumberofparameters.

By leveraging a core asynchronous state evolution mechanism based on device physics, the proposed hardware implementation offers an efficient pathway for event-driven, real-time processing. Theseresultssuggestapromisingapproachforaccurate,scalable,and highly-efficientneuromorphicsystemscapableofprocessingreal-time eventstreams.

Results
Modelarchitecture
Figure1ashowsthehigh-levelmodelarchitecture,whichisdesigned for asynchronous event processing. The input event-stream is an orderedsetofeventsE=f(t ,j )g,wheret isthetimestampandj is m m m m thechannelindexofeventm.Inaudiotasks,channelindicescorre- spondtodistinctfrequencybandswhereasinvisiontasks,channels map one-to-one with pixel locations. The channel-wise embedding layer, shown in Fig. 1b, maps an event into a D-dimensional vector based on the event’s channel index. Since the embedding matrix W 2RJ×D,whereJ isthenumberofinputchannels,isstatic embedding aftertraining,itcanbephysicallystoredinanRRAMcrossbararrayand embeddingisperformedbyreadingoutthevaluesintherowcorre- spondingtoj whenspikemarrives.Embeddingisacommonpre- m processingtechniqueusedinmodernmodels,e.g.,languagemodels. Inthiscase,bymappingeachchanneltoadistinctembeddingvector, we endow every channel with a unique feature representation, enabling the SSM blocks to process event streams within a unified featurespace.

Theembeddedeventrepresentationsarethenpassedthrougha seriesofstackedSSMblocks.Figure1cshowstheinternaloperations of a single SSM block. To simplify hardware implementation, we modify the complex-valued decompositions in the original SSM models to perform all operator factorizations in real domain (see Methods).Specifically,thecontinuous-timedynamicsofthestateare modeledbyafirst-orderdifferentialequation:

dh(t)/dt = \Lambda h(t) + B x(t), (1)

where h(t)2RH denotes the hidden state, x(t)2RHin is the input, \Lambda 2RH×Hisadiagonalmatrixencodingthestatetransitionrates,and B 2RH×Hin is the input projection matrix. An output matrix C 2 RH×Hout willbeusedtoconvertthehiddenstateintooutputvectorfor downstreamtasksortothenextlayer.

Toimplementthemodelinanevent-drivenframework,weapply a discretization method for the following state update equation:

h(t+ \Delta t) = \Lambda'(\Delta t)h(t) + B' x(t+ \Delta t), (2)

wherethestatetransitionmatrix \Lambda'(\Delta t)andtheinputprojectionmatrix B' aredefinedas

\Lambda'(\Delta t) = exp(\Lambda \Delta t), (3)

B' = \Lambda^{-1}(exp(\Lambda) - I)B (4)

whereIistheidentitymatrix.WenotethatEq.(2)canthenbeeffi- cientlyimplementedinanalogcircuitswherethe2ndtermcorresponds toaVMMfunctionandcanbeimplementedusinganRRAMcrossbar forthecurrentinputatt+ \Delta t,whilethe1sttermcanbeimplementedin parallel using a group of devices (termed state nodes, SNs here) exhibiting decay dynamics, taking advantage of the diagonalized matrix \Lambda.

AsillustratedinFig.1c,Whennoinputsareapplied(e.g.,during intervals \Delta t and \Delta t ),thestateevolvepassivelyandnativelythrough 1 2 the STM property of the SNs, implementing the exponential decay (e.g.,e^{- \Delta t1},e^{- \Delta t2})withouttheneedofclocksorsynchronization.Active computation is performed viaVMMand added to the current state, making the system operationfullyevent-driven.Theupdatedstateisthenreadout,fur- ther processed through a non-linear block, and applied as aninput eventtotriggerthenextblock’scomputation.Byleveragingphysics, including Ohm’s law, Kirchhoff’s current law, and internal device dynamics,thisapproachachieveshighlyparallelandenergy-efficient computationandiscompatiblewithasynchronousinputs.

Tofurtherfacilitatepracticalimplementation,akeydesignchoice inourarchitectureistolimitoneSSMBlocktothesamedecayratein thematrix \Lambda.Inotherwords,insteadoflearningindependentdecay ratesforeachdimensioninthestate,weset

\Lambda = \lambda I, (5)

where \lambda is a hyperparameter determined during training (SupplementaryNotes1).Atfirstglance,imposingasingleshareddecayrate acrossanentireblockmayseemexcessivelyrestrictive.However,in practice, we found that this constraint does not notably limit the model’srepresentationalcapacity.Thisislikelybecauseotherlearn- able parameters, such as the input projection matrix B and output matrix C,alongweightswithinthenonlinearblock,cancompensate the restrictions on \lambda. Furthermore, the stacking of multiple blocks provides additional mechanisms for modeling diverse temporal dynamics.Ifnecessary,theconstraintcanberelaxedbygroupingstate variablesintoafewclusterswithinoneblock,eachwithitsowndecay constant—i.e.,havingn(wheren«Hwhichisthedimensionofthestate space)perblock.Thisadjustmentcanbetterbalancetherequirements for practical hardware implementation and preserving sufficient expressivepowerforthemodel,althoughforthedatasetswetested havingasingletimeconstantperblockalreadyoffersexcellentmodel performance.

Interestingly, Eq. 2 closely resembles the dynamics of a Leaky Integrate (LI) neuron, in which the neuron’s membrane potential decays exponentially over time. This similarity suggests that SSMs couldbebiologicallyfeasible,providinganintriguingconnectionto neuralmechanismsobservedinnature.However,akeybenefitofSSM is that most operations are linear (e.g., Equation (1) only has linear terms), avoiding strong non-linearity such as spike generation and resetoperationsattheneuronlevelwhichisahallmarkforSNNs.This property makes it highly effective to train SSMs, as discussed in Orvietoetal.40.Theonlypositionalnon-linearfunctionisatthefinal outputofanSSMblock,intheformofagatedactivationfunctionsuch as:

h_e = h \odot \sigma(W \odot GELU(h) + b), (6)

whereGELU( \odot )istheGaussianErrorLinearUnitactivation,Wandbare learnedparameters, \sigma( \odot )denotesthesigmoidfunction,and \odot repre- sentselement-wisemultiplication.h_e isthencombinedwiththeoriginal stateviaaresidualconnection:

h_out = h + h_e. (7)

HardwareImplementationofSSMBlock
HavingdetailedthemathematicalframeworkoftheproposedSSM, wenowdescribeitsphysicalimplementationinaCIM-basedhard- ware.Weco-designedthehardwarearchitecturetonativelyexecute themodel’scoreoperations,followingtheprinciplesofperforming computingdirectlythroughdeviceandcircuitphysics.Specifically, we implement the state update equation by leveraging the native exponentialdecayoftheinternalstatevariableinmemristordevi- ces exhibiting STM effects39. For example, prior studies have shownthatwhenstimulatedbyaninput,theWOₓ-basedmemristor deviceconductanceincreasesaccordingtotheinputthensponta- neously relaxes due to oxygen ion diffusion, resulting in STM behavior39,41.Giventhatstateparametersevolveindependently(due totheuseofdiagonalstatedynamicsmatrix)inanSSMblock,the statecanberepresentedbyanarrayofsuchvectordecaydynamics can be parallelized and directly implemented using STM memristors(SNs).

Specifically,thestateupdateequationinsideanSSMblock(Eq. (2))canbeimplementedas:

G(t+ \Delta t) = exp(\Lambda \Delta t)G(t) + B' x(t+ \Delta t). (8)

Here G(t)=[g_0(t),g_1(t)...g_h(t)] is the collection of the STM memristor’s conductance at timet and represents the hidden state vector, \Lambda = \lambda I and \lambda is thelayer-wisefixeddecayconstant, \Delta tisthetime intervalsincepreviousevent,B' istheinputprojectionmatrixand x(t+ \Delta t)istheinputattimet+ \Delta twhenthecurrentspikearrives.We notethatthisequationcanbephysicallyimplementedusinganarray ofSTMmemristorsthatactasSNstostorethestateandimplementthe statedecaydynamics(1stterminEq.(8)),andanRRAMcrossbarthat implementstheVMMoperationbetweentheB' matrixandtheinput (2ndterminEq.(8)),asshowninFig.2a.Theproposedhardwarethus directlyimplementsthestateupdateequation(Eq.8),whichinturn correspondstothestateequationEq.2.

WetestedthisimplementationexperimentallyusinganRRAM CIM chip that performs the VMM operations and a group of WOₓ memristors that exhibit STM behaviors, respectively. The inset of Fig. 2B shows a photograph of the RRAM CIM chip for the VMM operations.Thischipwasfabricatedusinga65nmCMOSprocess with four integrated 64×64 RRAM arrays, which are logically combinedtoformanequivalent128×128array.EachRRAMdevice in the array is based on a one-transistor-one-resistor (1T1R) struc- turewithaphysicalsizeof200nm42.Eacharrayisconnectedwith integrated 8-bit digital-to-analog converters (DACs) for input encoding and 8-bit analog-to-digital converters (ADCs) for output readout.ThemeasuredVMMoutputsfromtheRRAMCIM chip are shown in Fig. 2b (dots), along with the ideal expected values(line).

Themeasuredoutputdistribution(bottominset,Fig.2b)dueto device variations is characterized by a standard deviation \sigma of 4.6 LSB (Least Significant Bit), corresponding to a Normalized Root- Mean-SquareError(NRMSE)of1.80%.Thehighconsistencybetween measuredandidealvaluesandlowNRMSEverifythattheRRAMCIM configurations can accurately performing the VMM operation requiredfortheB' x(t+ \Delta t)terminEq.8.MoredetailsoftheRRAM CIM chip can be found in Supplementary Fig. 1 and Supplemen- taryNote2.

The VMM output is then applied to the SNs to implement the exponentialstateevolution(1stterminEq.(8)).Figure2c,dshowsthe state evolution function implemented with WOₓ memristors that exhibitSTMbehaviors,wheretheexp(\Lambda \Delta t)G(t)functionispassively implemented based on the native conductance decay in between inputs.TheapplicationoftheVMMoutputupdatesthedecayedcon- ductancevalues,completingthestatetransitionbynaturallysumming thetwotermsofEq.8intheSTMmemristordevices.

Tomatchtherequireddecayprofiles,WOₓdeviceswithdifferent oxide thicknesses were fabricated by tuning the rapid thermal pro- cessingannealingduration(seeMethods),followingpriorstudies39,43. Longerannealingtimesyieldthickeroxidelayersandhigheroxygen vacancydensities,whichin turnslow the relaxation process due to more persistent conductivefilamentpaths.Moreinformationonthe tunabledecaybehaviorcanbefoundinSupplementaryFig.2.

Figure 2c, d show experimentally implemented Eq. 8 with two targettimeconstants.Strongagreementbetweenmeasuredoutputs (scatter points) and the expected SSM-based evolution (solid lines obtained by simulation of Eq. (2)), confirming the viability of this approachforhardwareSSMinference.

Byfixingthedecayconstanttoasinglevalueinablock,onlysix distinct decay profiles—one per layer—are required for hardware implementationoftheproposedSSMmodel.Thissignificantlyredu- cesthephysicalimplementationchallengesasprecisecontrolofthe timeconstant \tau toadifferenttargetvalueforeveryindividualmem- ristordeviceisimpracticalduringfabrication.

The updated state from each SN is then read out and used to perform the positional non-linear gating and residual connection operationsasdescribedinEqs.6and7.Theseoperations,including theGELUandsigmoidactivationfunctions,areimplementeddigitally using Look-Up Tables (LUTs) that store pre-computed values. The outputisthenpassedtothenextSSMblockasadiscreteinput.This entireprocess,fromeventreceptiontoSNstateupdate,demonstrates thetightcouplingbetweenthemathematicalformulationandphysical device dynamics, enabling a highly efficient and clockless imple- mentationoftheSSM.

Modeltrainingandperformanceanalysis
Wenowdetailthetrainingprocedureusedtodeterminethemodel’s optimalparameters.Wenotethatthetrainingprocessisconductedin simulation.Allstandardnetworkparameters—includingtheSSMinput projectionB,theoutputprojectionmatrixC,thegatingweightsand biases,and(initially)eachindividualdecayrate \lambda insideeachlayer—are trainedjointlywithbackpropagationthroughtime(BPTT)44,without anysurrogateapproximations(seeMethods).Toreconciletheneedfor per-dimensionflexibilityduringlearningwiththehardwareconstraint ofasingleshared \lambda perblock,wedetermine \lambda inthreestages.First,we traintheentirenetwork(includingeach \lambda)untilthetraininglosspla- teaus—allowingeachstatedimensiontodiscoveritsowndecayrate. Second,wecomputethearithmeticmeanofalllearned \lambda withineach layer.Finally,wefix \lambda tothislayer-wiseaverageandcontinuetraining the remaining parameters for all subsequent epochs. The detailed informationaboutthe \lambda foreachtaskisshowninSupplementaryFig.3 andSupplementaryTable1.

Wefirsttrainedandevaluatedmodelonfourevent-drivendata- sets spanning both spiking audioand event-based vision tasks. The spikingaudiodatasetsincludeSpikingHeidelbergDigits(SHD)5and SSC5.SHDcomprisesspokendigitsconvertedintospiketrainsthrough abiologicallyinspiredauditorymodel,witheachdigitrepresentedasa spatio-temporalpatternacrossmultiplefrequencychannels,whileSSC isderivedfromGoogle’sSpeechCommandscorpus,retainingword categorieslike“yes,”“no,”“up,”and“down,”butconvertingtheminto spiketrainsviaagammatone-likefilterbankfordenserandlengthier eventstreams.Theevent-basedvisiondatasetsareDVS128Gesture3 andDVS128Lips4,bothcapturedbyaDVScamera,producingasyn- chronousspikesequencesona128×128pixelgrid.DVS128Gesture encompassesdiversehandandarmmovementsthatdemandmodel- ingboththetemporalprogressionandlimitedspatialinformationof the gesture pattern, whereas DVS128 Lips highlights subtler lip movementsoftenusedforaudio-visualspeechprocessing,requiring fine-grainedtemporalencodingtodistinguishrapidmotioncues.To improverobustnessandgeneralizationunderreal-worldusecases,we appliedstandarddataaugmentationsduring training45.Thetraining losscurvesareshownSupplementaryFig.4.

Table 1 summarizes the simulation results on two prominent spikingaudiodatasets,whicharefrequentlytackledwithSNN-based architectures.AsshowninTable1,ourreal-valuedasynchronousSSM model with shared decay time constants consistently outperforms SNN-based baselines while maintaining a low parameter count. Notably, even though our method does not explicitly replicate the neuron-level spiking mechanism found in SNNs, it still excels at handling these sparse, event-driven audio signals. The results demonstratethatSSMswithdiagonalizedexponential-decaydynamics isbothrobustandeffectiveforlong-rangetemporaltasksconsistingof asynchronousspikingsequences.

Results on two representative event-based vision datasets— DVS128GestureandDVS128Lips—asshowninTable2.Unlikemany existingworksonDVSdata,wedonotincludecommonlydeployed spatialmodulessuchasconvolutionaloperations,sincetheseopera- tions(1)donotsupportasynchronousprocessing(asafullfeaturemap needstobecreatedfirsttosupportthescanningofinputpatches)and (2)arenotfriendlywithCIMhardware,asdiscussedearlier.Instead, our model relies on purely temporal SSM blocks, where the spatial features are encoded through the embedding process and the attention-like state updates34. Despite this choice of avoiding con- volutionoperations,ourapproachachievesaccuracyonDVS128Ges- turethatgenerallymatchesmorecomplexarchitecturesintegrating explicit spatial processing. On DVS128 Lips, we observe similarly competitiveperformancewithaparametercountsignificantlylower thanthatofleadingmethods.

Theseresultsunderscorethatanevent-drivenSSMframework— even with simplified real-valued dynamics—can generalize across vision tasks on sparse, asynchronous data without the overhead of spatialconvolutions.Equallyimportantly,theasynchronousnatureof ourSSMimplementationcanfullytakeadvantageofthespikinginputs and the parallel computing capabilities of CIM systems for highly efficienthardwareimplementation.

Itisnoteworthythatsomesynchronousmethods,byaggregating sparseeventsintoframe-likestructures,canachievehigheraccuracy onDVSdataset.However,thisperformancecomesatthecostoffor- goingthebenefitsofevent-basedsensing.Theframe-creationstepnot onlyintroducesadditionallatencybutalsoincreasesthedatavolume, leadingtomorecomputationaldemands.In contrast,the proposed asynchronousSSMprocesseseacheventasitarrives,preservingthe data’s inherent sparsity and temporal precision. This makes our approach suited for real-time, resource-constrained applications whereminimizinglatencyandpowerconsumptionisparamount.

To provide a quantitative basis for the hardware and latency advantagesofourapproach,weperformedadetailedcomputational complexityanalysis(seeMethods).FortheDVS128Gesturedataset,our Event-SSMrequires1.68GFLOPstoprocessatypicaldatasequence.In contrast, a common alternative approach involves converting the event stream intovideoframes(e.g.,at30fps)andusingastandard CNN.AResNet−18basedmodelprocessingtheseframeswouldrequire 104.28 GFLOPs (a 62-fold increase), and a ResNet-50 would require 219.80 GFLOPs (a 131-fold increase). This significant reduction in computational load stems from our model’s ability to leverage the inherentsparsityofinputeventdata,avoidingthemassivecomputa- tionaloverheadassociatedwithframe-basedprocessing.

We performed several additional analyses in simulation to analyze the model behavior. Figure 3a–d show how the classification accuracy evolves as we vary the number of stacked SSM blocks (model depth) for different event-driven datasets. All models tested employ the same underlying architecture differing mainly in the number of stacked SSM blocks after the embedding layer. We observe that increasing the depth initially enhances performance, reflecting the model’s ability to better capture complex temporal dynamics with more layers. However, this improvement reaches a saturation point, typically around 6 blocks, indicating that adding further SSM blocks beyond a certain threshold yields minimal or no accuracy gain. The saturation effect is common in neural networks-beyond an optimal depth, adding more layers may offer diminishing returns and can increase the risk of overfitting or optimization difficulties during training. Therefore, the design of our model architecture is based on the saturation point. We note the required depth (about 6 blocks) is remarkably shallow compared to many conventional deep networks. For example, Bulzomi et al.46 and Gehrig et al.47 use ResNet48 as the model backbone for the DVS datasets which employs more than 30 layers. The fact that SSMs can effectively model long-range event sequences with just 6 blocks underscore their representational efficiency and make them particularly attractive for hardware deployment, where model mapping complexity and latency scale roughly with layer count.

To evaluate the effects of constraining the decay parameter to a single value per block, we compare the final classification accuracies across four event-driven datasets-two audio datasets and two vision- based datasets under conditions where the decay parameter \lambda is either freely learned (unfixed) or held constant (fixed), shown in Fig. 3e. For most datasets, we observe that fixing \lambda following the approach discussed in the Model Training section surprisingly results in improved accuracy by approximately 1% on average. Specifically, SHD improves from 94.7% to 95.7%, DVS128 Gesture from 96.6% to 97.3%, and DVS128 Lips from 61.1% to 63.5%. Notably, the Spiking Speech Commands (SSC) dataset exhibits a critical failure when \lambda is unfixed, leading to divergent training with invalid Nota Number (NAN) loss values, yet it achieves a stable and competitive accuracy of 84.7% when \lambda is fixed. These observations indicate that constraining \lambda may in fact help stabilize training convergence by serving as a form of regularization for event-driven datasets, including those with challenging temporal structures such as SSC.

To further investigate the instability observed in SSC, we plot in Fig. 3f the distribution of \tau = 1/\lambda across the SSM layers for SSC and other datasets. Since training with unconstrained \lambda on the SSC dataset terminated prematurely due to NAN values, the \tau values plotted here for SSC were recorded from the final stable iteration immediately before divergence occurred. Notably, on SSC \tau values are significantly larger compared to the other datasets, indicating slower decay rates in the hidden states. From a continuous-time perspective, the hidden state updates following:

h(t+ \Delta t) = exp(- \Delta t / \tau) h(t), (9)

With larger \tau, the hidden state h decays more slowly, preserving more historical information over prolonged periods. While beneficial for capturing long-term dependencies, this slow decay can cause the hidden states to accumulate excessively when the input event rate is high, which is also the case for SSC. Consequently, gradients and loss values can become unstable, resulting in divergence and NAN values. This numerical instability illustrates the challenges associated with unconstrained optimization of \lambda on datasets with demanding temporal characteristics, and the necessity of carefully regularizing \lambda to maintain numerical stability and achieve convergence.

HardwareSystemPerformanceAnalysis
We next evaluate the model implementation on CIM hardware using the SSC and DVS-Gesture datasets. To align the SSM model with analog CIM hardware constraints, all weights and activations are quantized to 8-bit integer (INT8) using quantization-aware training49,50. This step causes no accuracy degradation for the DVS-Gesture model (remaining at 97.3%) and only a minimal 0.3% drop for the SSC model (from 84.7% to 84.4%).

We then studied the effect of device non-idealities in analog hardware by introducing random hardware variations to the INT8 parameters (see Method) based on the measured device variation effect data in Fig. 2b–d. This resulted in only a modest accuracy decline for both models without any re-training or other mitigation techniques: the SSC accuracy dropped 2.39% to 82.01%, and the DVS-Gesture accuracy dropped 2.2% to 95.1%. Since the STM memristors need to implement the trained \lambda, we also analyzed the effects of variations in the device decay rate on model performance (Supplementary Fig. 5 and Supplementary Note 3). These results show that the models exhibit negligible accuracy degradation with variations up to 10% of the mean, without any re-training or variation-aware training techniques. The noise and variation analyses support the resilience of the proposed model for practical analog hardware implementations, likely due to the use of only fully-connect layers and the shallow model architecture. Future implementation of techniques such as variation-aware-training50 is expected to further improve the hardware accuracy.

To estimate power consumption, we constructed simulation for the proposed hardware system for the SSC dataset, using experimentally measured data from the CIM chip and STM memristor devices for the SSC dataset. The system consists of 6 SSM blocks, as shown in Fig. 4a.

The total power consumption for our system is estimated to be 34mW and the energy distribution across system components is shown in Fig. 4b (see Method). The 8-bit ADCs dominate the power budget, accounting for 43.7% of total power usage. Analog power consumption in the RRAM arrays follow, consuming 30.3%. SN devices for update behaviors consume 0.4%. However, we note that the estimate was based on ADCs and other digital components constantly operating at 40MHz due to the limitations of the current CIM chip. In practice, when the operations are truly operated in an event-driven fashion, the ADCs will not be active all the time and actual power consumption can be significantly lower. The power can be further improved by eliminating the analog-digital transitions to implement all operations in the analog domain, which is compatible with SSM based approaches.

Methods
Modelinitialization
WefirstassembleamodifiedHiPPO–LegS57operatorbycomputinga lower-triangularmatrixwhose(i,j)entry(forj≤i)isgivenbythepro- ductof 1+2iand 1+2j,withitsdiagonalentriesshiftedbysub- tractingtherowindex.Torendertheoperatordiagonalizablebyan orthogonaltransform,wepauffiffigffiffimffiffiffiffiffieffiffinffiffititwitharank-onecorrectioncon- structedfromavectorof i+1=2termsandthenenforcestrictsym- metry by averaging the result with its transpose. The resulting real symmetricmatrixSadmitsanorthogonaleigen-decomposition

S=Vdiag(\Lambda)VT

inwhich \Lambda isguaranteedtoconsistsolelyofrealvalues.

Allsubsequentprojections—includingtheinputmappingBand any output weights—are initialized using conventional real-valued schemes (e.g., LeCun normal58) without splitting into real and ima- ginary parts. By maintaining this wholly real-valued pipeline—from operator construction through eigen-decomposition and parameter initialization—ourmethodproducesanSSMblockwhosestateupdates and nonlinear projections can be implemented directly on analog crossbar arrays and digital lookup tables, avoiding any reliance on complexhardwareunits.

Devicefabrication
TheCIMchipconsistsoffour64×64crossbararraysbasedon1T1R cellswithaphysicaldiameterof200nm42.Thesecrossbarswere monolithicallyintegratedaboveCMOScircuitryusinga65nmprocess. Eachtileincorporates8-bitDACsforinputencodingandshared8-bit outputdigitization.Adifferential(dual-column)schemeis employedtorepresentsignedweights,enablingsymmetricmapping andimprovedlinearity.Theperipheralcircuitsandcontrollogicare managedbyanon-chipRISC-Vprocessortoenablefullyintegratedand autonomousCIMoperation.TheCIMchipismeasuredusingacustom- designedtestboard.

TheWO memristordeviceswithasizeof1μm×1μmwerefab- ricated onSiO/Sisubstrate,following the methodology outlined in our previous work39,43. The bottom electrode (BE) patterns were definedbyopticalphotolithography(GCAAS200AutoStep),followed by 60nmW metal deposition by DC reactive sputtering using a W metaltargetatroomtemperature,andsubsequentlyformedusinga lift-off process. A 250nm thick SiO film was deposited by plasma- enhancedchemicalvapordeposition(PECVD),followedbydirectional etchingusingRIE(LAM9400)toexposetheBEtopsurfaceandcreate aspacerstructure.Rapidthermalannealing(RTA)wasperformedto formtheWO layer,byoxidizingtheexposedpartoftheWBEinO atmosphereat400°Cfor20–80s.Afterwards,theAu/Titopelectro- des(TEs)wereformedbyopticallithography,liftoffandE-beammetal deposition.AdditionalRIEprocesswasappliedtoselectetchoutthe WO exposethe BEs for electricalcontacts.The devices werechar- acterizedusingaKeithley4200Ssemiconductoranalyzerwith4225- RPMmodule.

ParallelizedBPTT
WetrainthemodelusingBPTT,computinggradientswithrespectto alltemporaldependenciesviaautomaticdifferentiation.Theforward pass implements the discrete-time SSM recurrence h(t+ \Delta t)= \Lambda'(\Delta t)h(t)+B' x(t+ \Delta t)forsequencesoflengthL.Ratherthan computing this recurrencesequentially,we leverage aparallelasso- ciative scanalgorithm thatexploits the associativity of linear recur- rencesthroughaoperator A,b : (A,b) \odot (A,b) = (A \odot A, A \odot b + b). This formulationenablessimultaneouscomputationofallhiddenstatesvia a balanced binary tree reduction, achieving O(logL) parallel time complexitycomparedtoO(L)forsequentialevaluationgivensufficient parallel processors59. During the backward pass, automatic differ- entiationthroughthescanoperationcomputesgradientsforallSSM parameters—includingstatematrix \Lambda,inputprojectionB,outputpro- jectionC,andlearnablediscretizationtimesteps—whilepreservingthe samelogarithmiccomplexity.ThisparallelizedBPTTimplementation allowsefficientgradient-basedoptimizationonlongeventsequences whilemaintainingfulltemporalcreditassignmentacrossalltimesteps.

TheoreticalFLOPsanalysis
WeconductedatheoreticalFLOPsanalysistocomparethecomputa- tionalcomplexityofEvent-SSMandResNetmodelsforDVS-Gesture dataprocessing.TheEvent-SSMmodelhas2stageseachcomposingof 3SSMblocks.Thestateinthefirststagehasthedimensionof128while the second has 256. We assume the input has sequences of 65,536 eventswhichistypicaltoDVS-Gesturedataset.

For models based on synchronous processing, event data are converted to 30fps video (approximately 180 frames for 6s), with each frame having a resolution of 128×128×2 channels (ON/OFF events). We analyzed both ResNet−18 and ResNet-50 as standard architectures.

Ouranalysisassumes:(1)multiply-accumulate(MAC)operations count as 2 FLOPs, (2) matrix multiplication H ×H×H counts as 2×H ×H×H FLOPs,and(3)normalizationoperationscountas5 FLOPsperelement.(4)ResNetprocesseseachframeindependently.

For Event-SSM, we calculate FLOPs for each layer including embeddingprojection(L×D×2),andforeachSSMlayer:inputpro- jection(L×H ×H×2),stateevolution(L×H ×2),outputprojection (L×H×H ×2), feedforward network (L×H2 ×2), and normal- ization(L×H×5),whereD,L,H ,H,H denoteembeddingdimen- sion,sequencelength,inputfeaturedimension,andstatedimension, outputfeaturedimensionrespectively.Poolingwithstride16reduces thesequencelengthfrom65,536to4,096afterthefirststage,andto 256afterthesecondstage.ForResNet,weusethestandardformula FLOPs=H ×W ×C ×(K2×C ×2)foreachconvolutionallayer, whereH ,W ,C ,K,C denotetheoutputheight,outputwidth, numberofoutputchannels,kernelsize,andnumberofinputchannels, respectively. Total FLOPs equal per-frame FLOPs multiplied by the numberofframes(180).

ResultsshowthatEvent-SSMrequires1.68GFLOPs,whileResNet- 18andResNet-50require104.28and219.80GFLOPs,respectively.This corresponds to 62:12× and 130:93× computational efficiency improvementsforEvent-SSM,primarilyattributedto(1)thesparsityof eventrepresentation(65,536eventsvs.5,898,240pixelvalues,a90× density difference), (2) elimination of fixed frame-rate processing overhead,(3)compactmodeldesign,and(4)real-valuedoperations throughout.

Noisemodeling
To evaluate the model’s robustness againsthardwarenon-idealities, we conducted a noise- injection simulation. In this simulation, we injected noise after each computationalstep, namelyVMMand the stateupdateintheSNs,tothephysicalhardware.Theinjectednoise wasmodeledusingaGaussiandistribution(N( \mu, \sigma^2))basedonvalues empirically derived directly from the experimental characterization data of the fabricated CIM chip and SN devices. This methodology ensuresthatoursimulationreflectstheanalognoiseprofileinherentin thephysicalsystem.

Poweranalysis
ThepowerfortheCIMchip,whichperformstheVMM,wasmeasured directlyonthePCBboard.Wemeasuredthepowerconsumptionwhile thechipwasperformingatypicalVMMoperationatitsoperational clock frequency of 40MHz. Energy consumption is calculated as Power Consumptionper array × number of arrays needed and the measured power consumption per array is about 0.028mW. It is importanttonotethat,duetothelimitationsofthecurrenttestchip, thismeasurementwasconductedunderaglobalclock,andtheclock frequencyismuchhigherthantherequirementofthedataset(After pooling, SSC dataset will have roughly 5000 events/second, corre- spondingtooperationatthekilohertzlevel).Wethereforebelievethat a future implementation using a fully event-driven CIM chip would exhibitevenlowerpowerconsumption,asitwouldtheoreticallyonly performVMMcomputationswhenaninputeventarrives,thusfully leveragingdatasparsity.

PowerconsumptionoftheSTMmemristorarrayisdynamicand event-driven. Therefore, we estimated the total energy E bysum- mingtheenergyofeachindividualinput.Weusedthemethodology fromYooetal.60,asshowninthefollowingequation:

E = \sum_{n=1}^N V_{pulse}^2 \cdot G_{x,n} \cdot \gamma_{pulse} \cdot t

where G_{x,n} represents the conductance of the specific SN devicebeingaddressed.ForeachSNupdates,whichcorrespondstoan inputspikex(t+ \Delta t)inEq.8,aburstofN1.4Vper50µspulseswere appliedonthememristor,wherethenumberofpulsesNispropor- tionaltotheVMMoutput.

Datafromtheremainingcomponentsareobtainedfromsynthe- sizedcircuitsusingTSMC28nmnode,followinganapproachbyWang et al.24 Inside each block, positive and negative weight values are mappedontoseparateRRAMarrays.The8-bitinputactivationsare applied in a bit-serial fashion using two voltage pulses, each repre- sentinga4-bitvalue(16discretelevels).Theresultinganalogcurrents fromRRAMarraysaredigitizedby8-bitADCsandsubsequentlyshift- addeddigitallytoobtainan8-bitpartialsum,representingtheoutput activationvaluesthatareappliedtothecorrespondingSNdevices.In ourmodel,thelayeratthefrontusesasmallerstatedimensionandthe layers at the back use a larger state dimension to allow for richer informationcompressedacrosstime.EachSSMblockalsorequiretwo LUTstoimplementthenon-linearfunctions,namely,GELUandSig- moid.Thepowerofthesenonlinearfunctionsaresimilarlysimulated. Giventhatthesefunctionsoperateonthe8-bitquantizedvaluesfrom ADCs, each non-linear operation can be efficiently realized using a compact LUT containing only 28 (256) pre-computed entries. The powerforoneLUTisestimatedtobe0.25mWat40MHz.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript entitled ‘Compute-in-Memory Implementation of State Space Models for Event Sequence Processing’ presents a SSM & CIM co-design method to implement efficient asynchronous processing for event sequences. This work is noteworthy in replacing the numeric iteration in the SSM with the short term memory effect of WOx devices, while using RRAM for VMM operations. The authors claimed that the proposed method outperforms existing SNN-based methods with respect to benchmark performance and power. However, the manuscript exhibits several weaknesses that require further improvement, as detailed below.\n1) The authors present their method as a real-time, event-driven hardware system. Therefore, it is expected to demonstrate a complete system with sufficient design details and measurement results. Critical parts should include the WOx-based SSM block, RRAM-based embedding layer and VMM modules, as well as how they work (e.g., how to implement line 338-341 physically?). Fig. 1 and 3 show several highly abstracted figures without details, making a big gap between the proposed concept and the hardware implementation. Moreover, the CIM chip of Fig. 3(A) is the only hardware image I can find in this paper, and again, without any details. It is suggested to provide a more thorough documentation of the hardware implementation to allow readers to follow.\n2) It is not clear whether the results in Table 1 and 2 were obtained through simulation or hardware measurement. Please clarify in both caption and main text. It is hard to “imagine” the deployment of millions of parameters onto the proposed system without knowing the implementation details. If it was done by simulation, please provide full details about the simulation model, e.g. how was the noise modelled?\n3) Following the last comment, please clearly present the data source throughout the manuscript. It is known the figure 3(C), (D) and S2 were experimentally measured, while fig 4(B) was simulated by TSMC PDK. What about other data, e.g. Fig 2? For the above reasons, it is hard to tell weather the results are technically sound.\n4) The authors are mixing multiple concepts, which may hinder reader comprehension. The term Reservoir Node should originate from reservoir computing. The authors previously published reservoir computing works with the same device WOx and the same decay property (ref. 58). As the decay factor (cid:1) in eq. 2 and 9 are determined by the WOx devices, it is fixed and less likely to be trained. Therefore the definition appears to be in line with reservoir computing. In fact, reservoir computing can be considered a SSM as the reservoir layer produces a state space. The authors should clarify the necessity of this new terminology and discuss how their SSM-based approach different from reservoir computing frameworks.\n5) The training process could be described in more detail. Please provide more details about the implementation of BPTT and how it handles the noisy dynamical behavior of WOx and RRAM? Are they differentiable? Also, it would be helpful to show the loss curve for the tasks.\n6) The authors claimed fully asynchronous processing. But, could the authors clarify how the claimed asynchronicity is maintained given that the ADC/DAC likely require clocks? The extent of asynchronicity needs to be clarified.\n7) In Line 374, the authors didn’t explain how to measure the power of RN devices and CIM chip. This is an important part that would benefit from more details.\n8) Both figures and text are poorly organized. The main problems are:\n(i) fig 1 (B) and (C) are somewhat unclear. Why the embedding vector is pointed to a straight line? It is not understood what exactly is the embedding. Both captions mention ‘internal structure’. The current figures may illustrate some functionalities but the structure is not clearly depicted.\n(ii) Insufficient labeling in figures obscures essential experimental/design details. The authors should label all the terms. For instance, what are delta t1 and delta t2 in Fig 1(c)? What does every arrow means in fig 4(A)? In fig. 4(B), the power distribution percentages quoted in the text (lines 387-389) do not match those shown in Figure 4(B). And What is LUT? It has not been mentioned throughout the text. In table 2, the acc. for DVS128 is 63.5%, that is, different from the 62.8% in line 280. Please double check the paper.\n(iii) It is suggested to totally reorganize the Results and Methods section. The majority of line 121-145 looks like a figure caption. And there are too many method descriptions in the Result section. Please first clear up what information should be in the Results, Methods, figure caption and SI, and then rewrite accordingly.\n(iv) The error of Fig 3(A) should be quantified using e.g. NRMSE.\n(v) The authors put the comparison table after fig 1. It is suggested to restructure the paper to follow a logical progression: problem -> algorithm -> hardware -> validation -> comparison.\nMinor:\nLine 142: simply -> simplified\nLine 326: double check the ref. 30. RN wasn’t found here.\nLine 374: system performance or power consumption?\nLine 40: Transformer block wasn’t found in ref. 11,12 and 13.", "reviewer_2": "This paper presents a methodology for implementing SSMs on CIM hardware, and holds significant potential by crossing over two emerging fields. In my view, this paper can enable significant impact. Some minor notes that could enhance the quality and technical accuracy of the manuscript are provided below. Note there are several references used to back up several points in this review. I do not expect the authors to include these citations.\nIntroduction\n- The authors state that the use of surrogate gradients is what compromises model accuracy. While this is broadly the accepted view and commonly cited reason, I'd contend that surrogate gradients are generally quite robust (especially in the era of QAT and extreme low-precision models that have been shown to scale up). I expect the main limiter is binarization of activations. This drives excessive information compression. This constraint is highlighted in Zhu et al.'s work in SpikeGPT (TMLR, 2024).\n- The paragraph that highlights the mismatch for convolutions and transformers for CIM makes a lot of sense. Buchel & Sebastian et al.'s work in \"Efficient scaling of large language models with mixture of experts and 3D analog in-memory computing\", Nat Comp Sci, 2025 is a good reference for this, where all Transformer/self-attention operations are delegated to digital circuitry. However, an additional challenge is that applying CIM to dense-layers alone does not address the primary bottleneck of Transformer ops for long sequence length. I.e., Transformers will dominate the computational cost relative to dense layers. This paragraph (or the start of the following paragraph) could be made more compelling by more clearly stating that SSMs overcome the memory-bottleneck of self-attention/KV-caching, making them a natural fit for CIM architectures.\n- There is some tightly related work where SSMs are deployed on neuromorphic/event-based hardware by Abreu et al. (Neuromorphic Principles for Efficient Large Language Models on Intel Loihi 2, ICLR SCOPE Workshop, 2025) & Zhu et al. (Scalable MatMul-Free Language Modeling on Neuromorphic Hardware, arXiv 2025). Neither uses CIM, and instead relies on digital, async accelerators.\n- I would encourage citing the original SSM+DL work by Voelker, Kajic & Eliasmith, \"Legendre memory units: Continuous-time representation in recurrent neural networks\", neurips, 2018, where referring to \"grounded in control theory\".\n- The use of real-valued coefficients in the transition matrix is acceptable, and the use of non-dynamic/input-dependent transition matrices for the update is a fair fit for CIM, though the dynamic transition matrix is a big part of what boosts performance in recurrent architectures. I think it is fair to state frankly that this is a limitation, although, if it was delegated to digital computation, it's a tiny fraction of the overall compute.\nResults\n- The graph in Figure 1(B) is a unclear to me (the square).\n- The discussion on nonlinearities could use the phrase \"positional nonlinearity\" which 1) preserves universal function approximation, 2) without breaking the ability to represent training as 1D-conv.\n- In model training, it's said that BPTT is used. Training an SSM is indeed functionally equivalent to BPTT, though some folk take issue with it as BPTT is defined as having linear time complexity while training SSMs follows log-complexity with time. I don't think there is concrete guidance here - I have no fixation on this if you prefer not to change this.\n- It is unclear to me where the asynchronousity in the SSM comes from. What is sparsifying the activations?\n- While not essential, an ablation study to identify what features of the RRAM/ADC non-idealities impact performance.\n- There is limited information on how sparsity provides an advantage in this particular set up.", "reviewer_3": "This manuscript presents an algorithm-hardware co-design approach for efficient, real-time, and asynchronous processing using state space models (SSMs). The key contributions can be summarized as follows: (1) a re-parameterized 6-block SSM model in which the decay parameter (λ) is shared within each block, thereby reducing hardware complexity while improving accuracy; (2) a memristor-based compute-in-memory (CIM) architecture to accelerate input vector-matrix multiplications (VMMs); and (3) an exploitation of the intrinsic short-term conductance decay behavior of memristors to perform hidden state (h) updates. However, several parts of the manuscript lack clarity and require further elaboration and justification.\n1. Recent References on RRAM CIM\nThe reference list regarding RRAM CIM could be strengthened by including recent works published in Nature and Nature Electronics. Incorporating these high-impact studies would better position the manuscript within the current state of the art and highlight its relevance to ongoing developments.\n2. Crossbar vs. 1T1R Structures\nThe manuscript should address the impact of sneak current in crossbar arrays and provide a discussion comparing the advantages and trade-offs of crossbar versus 1T1R structures. Such an analysis would help readers better appreciate the design choices and their implications for scalability and reliability.\n3. Figure 1(A) & (B): Clarity of Asynchronous Event Representation\nThe presentation of asynchronous events through the channel-wise embedding layer is unclear. For instance, the y-axis values are not specified, and it is difficult to determine whether the spikes in (A) correspond directly to the square wave in (B). Providing clarification or additional explanation would greatly improve readability.\n4. Synchronous vs. Asynchronous Processing\nThe manuscript emphasizes the importance of asynchronous processing, but the distinction from synchronous processing is not sufficiently detailed. The authors are encouraged to elaborate on the differences, including their respective advantages and disadvantages. For example, in Table 2, Peng et al. [50] report comparable performance despite not being asynchronous. Enhancing this comparison table to highlight the significance of the proposed approach would strengthen the manuscript.\n5. Figure 2(A–D): Caption and Interpretation\nThe caption states “Analysis of model architecture and decay parameter (λ) choices on classification accuracy,” yet the figure appears to show the number of blocks versus accuracy. Does the caption misstate the content, or is there a specific relationship between λ and the number of blocks? Furthermore, the authors claim that the proposed model requires only six blocks, whereas counterparts such as ResNet require more than 30 layers. Given the distinction between the terms “block” and “layer,” is it valid to conclude that the proposed model achieves lower latency? A direct comparison in terms of latency or other performance metrics would make this claim more convincing.\n6. RRAM Mapping Strategy\nThe authors map positive and negative weights onto separate RRAM arrays, while prior work often employs 2’s complement mapping, requiring only a single array. What motivated the use of separate arrays in this work? Discussion of the resulting area overhead and associated trade-offs would add valuable context.\n7. Fixed Decay Parameter (λ) and Accuracy\nLine 276 states that inference accuracy improves when fixing the decay parameter (λ). The authors should provide a more detailed mathematical derivation or theoretical explanation to clarify why and how fixing λ enhances accuracy.\n8. Flexibility Limitations of Memristor Conductance Decay\nUsing memristor conductance decay as RN is an intriguing approach. However, because the decay rate is largely determined by oxide thickness during fabrication, flexibility may be limited for supporting different SSM models or enabling weight updates. The authors should elaborate on potential solutions to overcome this limitation, particularly given the importance of reprogramming capabilities at the edge.\n9. Voltage Levels and Endurance in Conductance Updates\nIn Equation (9), the horizontal line voltages of the input projection memristor array are used as input spikes x(t+Δt) to update conductance. The authors should clarify the voltage levels applied in this process. Moreover, since each conductance update effectively involves multiple write operations, a discussion of endurance concerns under repeated updates would be valuable."}, "author_response": ""}, {"round": 2, "reviewer_comments": {"reviewer_1": "I thank the authors for having carefully taken my comments into account. The authors made substantial improvements on the manuscript. The claimed points are now endorsed by enough details and results. My only two minor comments are: 1) please avoid using some strange wordings like ‘software-level’. Do you mean simulation? If so please simply use ‘simulation’. 2) Now I understand that the main difference between the proposed method and reservoir computing is that the decay factor lambda of the devices were also trained, providing extra flexibilities for the system to better adapt to the tasks. My only concern is that was the lambda trained offline and then somehow configured the devices to the values of target lambda? Please clarify. If so, how to overcome the device variabilities and sim2real problems?\nApart from these minor comments, I can recommend publication in Nature communications.", "reviewer_2": "The authors have fully addressed all of my concerns. I expect this paper will help advance the field.", "reviewer_3": "The previous comments have been addressed by the authors.\n1. However, I am wondering why the authors claim that the transformer is not compatible with CIM, given that numerous CIM studies have been proposed for transformer-based architectures.\n2. Additionally, regarding the decay parameter (λ), the authors mention that a fixed value is used. Could the authors discuss whether this parameter may vary due to temperature fluctuations or temporal decay, and how such variations might affect the application performance or computational accuracy?"}, "author_response": "We thank all reviewers for their valuable and constructive comments on this paper. We are pleased to learn that all reviewers agreed our revisions addressed their previous comments. In the revised manuscript we addressed the additional valuable suggestions and comments by the reviewers. The detailed point-by-point response can be found below.\nReviewer #1:\nComments to the Author:\nI thank the authors for having carefully taken my comments into account. The\nauthors made substantial improvements on the manuscript. The claimed points are\nnow endorsed by enough details and results.\nMy only two minor comments are:\n1. Please avoid using some strange wordings like ‘software-level’. Do you mean\nsimulation? If so please simply use ‘simulation’.\nResponse:\nWe thank reviewer for the positive feedback and the valuable comments.\nFollowing the reviewer’s suggestion, we have changed the text correspondingly:\nIn p8\n‘We note that the training process is conducted in simulation.’\nIn p10:\n‘Table 1 summarizes the simulation results on two prominent spiking audio datasets,\nwhich are frequently tackled with SNN-based architectures.’\nIn p13:\n‘We performed several additional analyses in simulation to analyze the model behavior.’\nAnd the caption of the Table 1 and Table 2 are also changed correspondingly:\n‘Table 1: Comparison of our SSM model to the state-of-the-art SNN based models on Spike Audio Datasets\nin simulation’\n‘Table 2: Comparison of our SSM model to other leading models on the DVS128 Gesture and DVS128 Lips\ndatasets in simulation.’\n2. Now I understand that the main difference between the proposed method and\nreservoir computing is that the decay factor lambda of the devices were also\ntrained, providing extra flexibilities for the system to better adapt to the tasks. My\nonly concern is that was the lambda trained offline and then somehow configured\nthe devices to the values of target lambda? Please clarify. If so, how to overcome\nthe device variabilities and sim2real problems?\nResponse:\nThank you for this question. Your understanding is correct: is determined during\noffline training, and we then implement these target values through physical STM\n𝜆\nmemristor devices.\nTo address the concern about variations in the STM memristor devices, we conducted\na new simulation study to evaluate the model's robustness against variations in .\nThis result is included in the revised manuscript as Supplementary Information S5.\n𝜆\nCopied below:\n‘S5. Model Performance against Decay Constant (λ) Variations.\nFig S5. Model accuracy as a function of device variations in the decay constant . The decay constant of the STM\nmemristor devices are modeled to follow a Gaussian distribution centered on the target mean value, and the standard\n𝜆 𝜆\ndeviation of the distribution is increased from 5% to 20% of the mean. Error bars represent results from 40 independent\ntest runs at each variation level. The dotted line represents the ideal baseline accuracy without variation from the\nquantized-weight model.\nPhysical device fabrication is subject to inherent process variability, meaning the actual\ndecay constant of a fabricated STM memristor device will likely deviate from its\nintended target value. To analyze the model's sensitivity to device-to-device variations,\n𝜆\nwe modeled the device decay constant λ as following a Gaussian distribution centered on\nthe layer's target mean . The standard deviation of the distribution was increased from\n5% to 20% of the target mean . The model’s accuracy was then evaluated by randomly\n𝜆\n𝜆\nassigning values to the STM memristor devices following the distribution, and\nrepeated 40 times for each distribution to obtain statistically significant results.\n𝜆\nThe results, plotted in Fig. S5, show that the model maintains high accuracy even with\nsubstantial variations up to 10% of the target mean . This result suggests that the\nproposed network can tolerate inherent variability of physical devices, confirming the\n𝜆 𝜆\nfeasibility of our hardware co-design approach.’\nWe also added relevant discussions to the “Hardware System Performance Analysis” Section\nin p14:\n“Since the STM memristors need to implement the pre-trained λ, we also analyzed the effects\nof variations in the device decay rate on model performance (Supplementary\nInformation S5). These results show that the models exhibit negligible accuracy\ndegradation with variations up to 10% of the mean, without any re-training or variation-\naware training techniques.”\nReviewer #2:\nThe authors have fully addressed all of my concerns. I expect this paper will help\nadvance the field.\nResponse:\nThanks for your positive comments.\nReviewer #3:\nThe previous comments have been addressed by the authors.\nResponse:\nThank you for your positive comments.\n1. However, I am wondering why the authors claim that the transformer is not\ncompatible with CIM, given that numerous CIM studies have been proposed for\ntransformer-based architectures.\nResponse:\nThank you for raising this question which helps us further clarify the challenges for\nimplementing Transformer models in CIM. In the revised manuscript, we specifically\nexplained the mismatch between the attention mechanism in Transformers and the CIM\nsystem:\nIn p2:\n‘In Transformers, the attention mechanism involves the multiplications of the Query, Key,\nand Value matrices, all of which depend on the input so the attention operation is not\nfriendly to stationary CIM arrays which assume one of the matrix (i.e. weights) is static.\nConsequently, practical implementations of Transformers adopt hybrid designs: they\nleverage stationary CIM only for VMM involving static weights (e.g., in Fully Connected\nLayers), but delegate the input-dependent attention score computations or the KV cache\nmanagement to separate digital circuitry27.’\n27. Büchel, J. et al. Efficient scaling of large language models with mixture of experts and 3D analog\nin-memory computing. Nat Comput Sci 1–14 (2025) doi:10.1038/s43588-024-00753-x.\n2. Regarding the decay parameter (λ), the authors mention that a fixed value is used.\nCould the authors discuss whether this parameter may vary due to temperature\nfluctuations or temporal decay, and how such variations might affect the application\nperformance or computational accuracy?\nResponse:\nWe thank the reviewer for this valuable comment. Since in our implementation λ is\ndetermined during offline training and then implemented through physical STM\nmemristor devices, deviations from the target value can potentially affect the physical\nsystem’s performance.\nTo address the concern about the variation in the STM memristor devices, we\nconducted a new simulation study to evaluate the model's robustness against\nvariations in λ. This result is included in the revised manuscript as Supplementary\nInformation S5. Copied below:\n‘S5. Model Performance against Decay Constant (λ) Variations.\nFig S5. Model accuracy as a function of device variations in the decay constant . The decay constant of the STM\nmemristor devices are modeled to follow a Gaussian distribution centered on the target mean value, and the standard\n𝜆 𝜆\ndeviation of the distribution is increased from 5% to 20% of the mean. Error bars represent results from 40 independent\ntest runs at each variation level. The dotted line represents the ideal baseline accuracy without variation from the\nquantized-weight model.\nPhysical device fabrication is subject to inherent process variability, meaning the actual\ndecay constant of a fabricated STM memristor device will likely deviate from its\nintended target value. To analyze the model's sensitivity to device-to-device variations,\n𝜆\nwe modeled the device decay constant λ as following a Gaussian distribution centered on\nthe layer's target mean . The standard deviation of the distribution was increased from\n5% to 20% of the target mean . The model’s accuracy was then evaluated by randomly\n𝜆\nassigning values to the STM memristor devices following the distribution, and\n𝜆\nrepeated 40 times for each distribution to obtain statistically significant results.\n𝜆\nThe results, plotted in Fig. S5, show that the model maintains high accuracy even with\nsubstantial variations up to 10% of the target mean . This result suggests that the\nproposed network can tolerate inherent variability of physical devices, confirming the\n𝜆 𝜆\nfeasibility of our hardware co-design approach.’\nWe also added relevant discussions to the “Hardware System Performance Analysis” Section\nin p14:\n“Since the STM memristors need to implement the pre-trained λ, we also analyzed the\neffects of variations in the device decay rate on model performance (Supplementary\nInformation S5). These results show that the models exhibit negligible accuracy\ndegradation with variations up to 10% of the mean, without any re-training or variation-\naware training techniques.”\nThe reviewer also mentioned the effects of temperature, which can indeed also affect\nthe memristor behavior (and affect analog computing devices and circuits in general).\nUnlike device variations that are random, temperature effects cause systematic errors\nand can potentially be modeled and compensated. In the revised manuscript, we\nacknowledged potential temperature effects and proposed how these systematic\neffects can potentially be compensated.\nIn p16:\n‘… recent device-level research demonstrates, the time constant of the SN can be actively\nmodulated in-situ by an external control signal, such as a bias voltage66. Such\napproaches make it possible to dynamically adjust the state decay rates on-the-fly and\neffectively implement modern selective SSMs in the proposed architecture. Additionally,\nsystematical shifts, such as temperature effects on the device’s decay constant, can\npotentially be compensated by techniques such as bias voltage control, enhancing the\nsystem’s real-world robustness.’\n66. Li, Y. et al. Homogeneous Memristors with Tunable Decay Dynamics for Self-\nAdaptive Reservoir Computing. ACS Nano 19, 28186–28196 (2025)."}, {"round": 3, "reviewer_comments": {"reviewer_1": "Thank you for your efforts. I have no further comments and recommend publication.", "reviewer_2": "The authors have addressed all of my previous points. The manuscript is in strong shape and should have a considerable degree of impact.", "reviewer_3": "I thank the author for addressing my comments. I have no further questions."}, "author_response": ""}]

Decision letter