WO2024243136A2

WO2024243136A2 - Systems and methods of scalable convex recurrent neural network training for interpreting brain dynamics

Info

Publication number: WO2024243136A2
Application number: PCT/US2024/030221
Authority: WO
Inventors: Fatih DINC; Hidenori Tanaka
Original assignee: Harvard University; Leland Stanford Junior University; NTT Research Inc
Current assignee: Harvard University; Leland Stanford Junior University; NTT Research Inc
Priority date: 2023-05-19
Filing date: 2024-05-20
Publication date: 2024-11-28
Anticipated expiration: 2025-11-19
Also published as: WO2024243136A3

Abstract

Computer-implemented methods and systems for training a neural network to interpret neural dynamics, the method comprising obtaining functional connectivity between brain regions and artificial neurons simulated by the neural network, capturing neural activities between biological neurons, processing the captured neural activities to obtain training data for the neural network, and constructing the neural network using the training data, to reproduce neural dynamics between the artificial neurons.

Description

PATENT COOPERATION TREATY APPLICATION FOR SYSTEMS AND METHODS OF SCALABLE CONVEX RECURRENT NEURAL NETWORK TRAINING FOR INTERPRETING BRAIN DYNAMICS INVENTOR: Hidenori Tanaka, Belmont, MA, USA Faith Dinc, Stanford, CA, USA Attorney Docket No.: NTTR-2023003pct

Attorney Docket No.: NTTR-2023003pct [0001] SYSTEMS AND METHODS OF SCALABLE CONVEX RECURRENT NEURAL NETWORK TRAINING FOR INTERPRETING BRAIN DYNAMICS [0002] CROSS-REFERENCE TO RELATED APPLICATIONS/PRIORITY [0003] The present invention claims priority to United States Provisional Patent Application Number 63/503,324 filed May 19, 2023, which is incorporated by reference into the present disclosure as if fully restated herein. Any conflict between the incorporated material and the specific teachings of this disclosure shall be resolved in favor of the latter. Likewise, any conflict between an art-understood definition of a word or phrase and a definition of the word or phrase as specifically taught in this disclosure shall be resolved in favor of the latter. [0004] BACKGROUND Understanding how the brain works is a complex problem that researchers so far have attacked from various angles. Despite the many beneficial applications of such understanding, current technologies fail for many reasons. There remains a real and pressing a need for the networks to scale to 10000s of neurons and be able to learn and predict in real-time support of brain-machine interfaces. [0005] SUMMARY [0006] Wherefore, it is an object of some embodiments of the presently disclosed invention to overcome some or all of above-mentioned shortcomings and drawbacks associated with the current technology. [0007] The presently disclosed invention relates to computer-implemented methods for training a neural network to interpret neural dynamics, the method comprising obtaining functional connectivity between brain regions and artificial neurons simulated by the neural network, capturing neural activities between biological neurons, processing the captured neural activities to obtain training data for the neural network, and constructing the neural network using the training data to reproduce neural dynamics between the artificial neurons. In some embodiments the neural Attorney Docket No.: NTTR-2023003pct network includes a convex solver. In some embodiments the convex solver comprises a convex weighted loss function. In some embodiments, computer-implemented method further comprises optimizing a loss of the convex solver via a fixed-point estimation. In some embodiments, computer-implemented method further comprises predicting a spatial and temporal activity of neural population dynamics. In some embodiments the neural network is configured to extract large-scale functional connectivity between brain regions and/or simulated networks. In some embodiments the neural network is trained to predict the temporal activity of a corresponding biological neuron in the brain. In some embodiments the neural network is configured to uncover attractor structures within the brain dynamics. In some embodiments the neural network is further configured to design external interventions to control neural population dynamics in the brain. In some embodiments the neural network is trained to interpret neural dynamics in real-time to support brain-machine interfaces. In some embodiments wherein the convex solver is further configured to perform optimization updates accelerated on a graphics processing unit (GPU). In some embodiments the convex solver is part of an alternating direction method of multipliers (ADMM) framework. In some embodiments the convex solver is configured to enforce constraints on the neural network during training. In some embodiments the constraints include preventing self-excitation of neurons within the neural network. In some embodiments the convex solver is configured to utilize biological priors for inductive biases during training. In some embodiments the convex solver is configured to handle input noise and conversion noise during the training of the neural network. In some embodiments the convex solver is configured to perform outlier detection during the training process. In some embodiments the convex solver is configured to scale to datasets collected from more than 10,000 neurons for offline processing. [0008] The presently disclosed invention further relates to non-transitory computer-readable mediums comprising instructions that, when executed by one or more processors, cause a computer system to obtain functional connectivity between brain regions and artificial neurons simulated by the neural network, capture neural activities between biological neurons, process the captured neural activities to obtain training data for the neural Attorney Docket No.: NTTR-2023003pct network, and construct the neural network using the training data to reproduce neural dynamics between the artificial neurons. [0009] The presently disclosed invention further relates to a computer system for training a neural network to interpret neural dynamics, the system comprising one or more processors and a non-transitory computer-readable medium including instructions that, when executed by the one or more processors, cause the computer system to obtain functional connectivity between brain regions and artificial neurons simulated by the neural network, capture neural activities between biological neurons, process the captured neural activities to obtain training data for the neural network, and construct the neural network using the training data to reproduce neural dynamics between the artificial neurons. In some embodiments the neural network is a data constrained recurrent neural network. [0010] The presently disclosed invention is directed to methods and apparatuses systems and methods that relate to a scalable convex recurrent neural network trained for interpreting brain dynamics are described herein. These systems and methods allow extraction of large-scale functional connectivity between brain regions and/or simulated networks. Embodiments of the present disclosure may include a computer- implemented method for training a neural network to interpret neural dynamics, including obtaining functional connectivity between brain regions and artificial neurons simulated by the neural network, including capturing neural activities between biological neurons. Embodiments may also include processing the neural activities between biological neurons to obtain training data for the neural network. [0011] Embodiments may also include constructing the neural network using the training data. In some embodiments, the neural network may include a convex solver. Embodiments may also include reproducing neural dynamics between the artificial neurons simulated from the neural network. In some embodiments, the neural network may be a data constrained recurrent neural network. [0012] In some embodiments, the convex solver includes a convex weighted loss function. In some embodiments, the method may include optimizing a loss of the convex solver via Attorney Docket No.: NTTR-2023003pct a fixed-point estimation. In some embodiments, the method may include predicting a spatial and temporal activity of neural population dynamics. [0013] In some embodiments, the convex solver includes a convex weighted loss function. In some embodiments, the instructions, when executed, cause the one or more processors to further perform optimizing a loss of the convex solver via a fixed-point estimation. In some embodiments, the instructions, when executed, cause the one or more processors to further perform including predicting a spatial and temporal activity of neural population dynamics. [0014] Embodiments of the present disclosure may also include a non-transitory computer readable medium including instructions to cause the one or more processors to perform obtaining a functional connectivity between brain regions and artificial neurons simulated by the neural network, including capturing neural activities between biological neurons. Embodiments may also include processing the neural activities between biological neurons to obtain training data for the neural network. [0015] Embodiments may also include constructing the neural network using the training data. In some embodiments, the neural network may include a convex solver. Embodiments may also include reproducing neural dynamics between the artificial neurons simulated from the neural network. In some embodiments, the neural network may be a data constrained recurrent neural network. [0016] Embodiments of the present disclosure may also include a computer system for training a neural network to interpret neural dynamics, including one or more processors. Embodiments may also include a non-transitory computer readable medium including instructions, when executed by the one or more processors, causing the computer system, to obtain a functional connectivity between brain regions and artificial neurons simulated by the neural network, including capturing neural activities between biological neurons. [0017] Embodiments may also include processing the neural activities between biological neurons to obtain training data for the neural network. Embodiments may also include constructing the neural network using the training data. In some embodiments, the neural network may include a convex solver. Embodiments may also include Attorney Docket No.: NTTR-2023003pct reproducing neural dynamics between the artificial neurons simulated from the neural network. In some embodiments, the neural network may be a data constrained recurrent neural network. [0018] In some embodiments, the convex solver includes a convex weighted loss function. In some embodiments, the instructions, when executed, cause the one or more processors to further perform optimizing a loss of the convex solver via a fixed-point estimation. [0019] In some embodiments, the instructions, when executed, cause the one or more processors to further perform the step of predicting a spatial and temporal activity of neural population dynamics. [0020] The foregoing summary is illustrative only and is not intended to be in any way limiting. Features from any of the disclosed embodiments can be used in combination with one another, without limitation. In addition to the illustrative aspects, embodiments, and features described above, further aspects, embodiments, and features will become apparent by reference to the drawings and the following detailed description. [0021] Various objects, features, aspects, and advantages of the present invention will become more apparent from the following detailed description of preferred embodiments of the invention, along with the accompanying drawings in which like numerals represent like components. The present invention may address one or more of the problems and deficiencies of the current technology discussed above. However, it is contemplated that the invention may prove useful in addressing other problems and deficiencies in a number of technical areas. Therefore, the claimed invention should not necessarily be construed as limited to addressing any of the particular problems or deficiencies discussed herein. [0022] BRIEF DESCRIPTION OF THE DRAWINGS [0023] The accompanying drawings, which are incorporated in and constitute a part of the specification, illustrate various embodiments of the invention and together with the general description of the invention given above and the detailed description of the drawings given below, serve to explain the principles of the invention. It is to be Attorney Docket No.: NTTR-2023003pct appreciated that the accompanying drawings are not necessarily to scale since the emphasis is instead placed on illustrating the principles of the invention. The invention will now be described, by way of example, with reference to the accompanying drawings in which: [0024] Fig. 1 is a schematic representation showing interpreting and controlling brain dynamics with data constrained recurrent neural networks (dRNNs). The inventors present a large-scale functional connectivity task for artificial recurrent neural networks (RNNs), where each artificial neuron is trained to predict the temporal activity of a corresponding biological neuron in the brain. By extracting RNNs with functional connectivity, the inventors aim to achieve three goals in future work: (i) predict neural population dynamics, (ii) uncover attractor structures, and (iii) design external interventions to control neural population dynamics in the brain. [0025] Fig. 2 is a schematic representation showing CoRNN: Convex and scalable solver for dRNNs via ADMM framework. The flowchart illustrates their CoRNN algorithm, which optimizes the parameters of a recurrent neural network to match the activity of a hidden unit with that of a generator (brain). The loss landscape is convex due to their careful choice of objective function. The algorithm starts by finding a fixed point for initialization and computing the inverse of the Hessian. It then starts a loop that goes through a certain number of iterations. At each iteration, the algorithm predicts the objectives, calculates the prediction error, updates the parameter values, and checks for convergence. The final output is the optimized parameters, which represent the functional connections between the neurons in the recurrent neural network. [0026] Figs. 3A and 3B are two graphs showing CoRNN achieves the best combination of speed and accuracy, dramatically accelerating convergence. The plot illustrates the relationship between accuracy (correlation between ground truth and inferred weights) and training time, measured in seconds on a log scale. The blue line represents the performance of their newly developed method, which combines fixed- point initialization and CoRNN, and achieves high accuracy while significantly accelerating convergence compared to other methods. Parameters: Fig.3A: α = 0.1, nrec Attorney Docket No.: NTTR-2023003pct = 200, T = 3000, ^^conv ! N (0, 10⁻⁴), ^^input ! N (0, 10⁻²) and Fig.3B: α = 0.1, nrec = 200, T = 3000, ^^conv ! N (0, 10⁻⁵), ^^input ! N (0, 10⁻¹). For both, error bars are standard error of the mean (s.e.m.) over 10 randomly initialized networks. [0027] Figs. 4A and 4B are two graphs showing CoRNN is much more stable and tolerant to variations in hyperparameter choices than FORCE. This plot describes the relationship between accuracy (correlation between ground truth and inferred weights) and the strength of regularization applied to the weights during training. The results show that while the FORCE algorithm, as previously applied to a similar task, is highly sensitive to the choice of hyperparameters, CoRNN maintains high accuracy across scales regardless of the strength of regularization. Fig.4A with fixed point initialization for all and Fig.4B with random initialization for Pytorch models and FORCE, whereas CoRNN is fixed-point initialized to define comparable scales. Parameters are the same as Fig. 3A. Error bars are sem. [0028] Fig.5 is a graph that shows ablating the fixed-point initialization degrades the convergence speed of the resulting models. The inventors conducted an ablation study to investigate the effectiveness of their proposed fixed-point initialization method. The results show that the CoRNN algorithm with fixed-point initialization (blue) outperforms all other methods, including the vanilla CoRNN without this fixed- point initialization technique. Parameters: α = 0.1, n_rec = 500, T = 10000, ^^conv ! N (0, 10⁻⁵), ^^input ! N (0, 10⁻¹). Error bars are sem. [0029] Figs. 6A – 6C are experimental results showing CoRNN can reproduce an RNN performing 3-bit flip flop task through dRNN training. The inventors use CoRNN to reproduce the behavior of a synthetic RNN trained on a k-bit flip-flop task, where the network is given a short input pulse at a random time. If the input pulse is +1, the output of the network must transition to +1 or stay at +1. If the input pulse is -1, the network’s output must transition to -1 or stay at -1, ignoring input pulses that do not match. Through dRNN training, CoRNN is able to accurately reproduce (Fig. 6A) the inputs and outputs of the task-trained RNN, (Fig. 6B) spatiotemporal pattern of the neural population dynamics, and (Fig. 6C) the inputs and recurrent weights of the synthetic RNN. Parameters: α = 0.9, n_rec = 1000, T = 200 trials (100 data points per Attorney Docket No.: NTTR-2023003pct trial), ^^conv ! N (0, 10⁻³), ^^input ! N (0, 10⁻²). [0030] Figs.7A and 7B are two graphs showing Increasing Trial Data is Necessary for Learning from Larger Networks. Fig.7A shows reconstruction accuracy of network outputs, measured as the correlation between the output units, and Fig. 7B shows reconstruction accuracy of network weights, measured as the correlation between recurrent weights, are plotted as a function of dataset size (number of trials) averaged over three networks. The different colors in the plot correspond to different number of neurons in the generator. Parameters: α = 0.9, ^^conv ! N (0, 10⁻³), ^^input ! N (0, 10⁻²). [0031] Fig. 8 is two plots showing CoRNN can scale to the largescale datasets. Both CPU and GPU implementation of CoRNN are fast and scale faster than data generation process itself. Parameters: α = 0.9, ^^conv !N (0, 10⁻³), ^^input ! N (0, 10⁻²). [0032] Fig.9 is a block diagram illustrating a computing device, according to some embodiments of the present disclosure. [0033] Fig. 10 is a flowchart illustrating a computer-implemented method for training a neural network, according to some embodiments of the present disclosure. [0034] DETAILED DESCRIPTION [0035] The present invention will be understood by reference to the following detailed description, which should be read in conjunction with the appended drawings. It is to be appreciated that the following detailed description of various embodiments is by way of example only and is not meant to limit, in any way, the scope of the present invention. In the summary above, in the following detailed description, in the claims below, and in the accompanying drawings, reference is made to particular features (including method steps) of the present invention. It is to be understood that the disclosure of the invention in this specification includes all possible combinations of such particular features, not just those explicitly described. For example, where a particular feature is disclosed in the context of a particular aspect or embodiment of the invention or a particular claim, that feature can also be used, to the extent possible, in combination with and/or in the context of other particular aspects and embodiments of the invention, and in the invention generally. The terms “comprise(s),” “include(s),” Attorney Docket No.: NTTR-2023003pct “having,” “has,” “can,” “contain(s),” and grammatical equivalents and variants thereof, as used herein, are intended to be open-ended transitional phrases, terms, or words that do not preclude the possibility of additional acts or structures. are used herein to mean that other components, ingredients, steps, etc. are optionally present. For example, an article “comprising” (or “which comprises”) components A, B, and C can consist of (i.e., contain only) components A, B, and C, or can contain not only components A, B, and C but also one or more other components. The singular forms “a,” “and” and “the” include plural references unless the context clearly dictates otherwise. Where reference is made herein to a method comprising two or more defined steps, the defined steps can be carried out in any order or simultaneously (except where the context excludes that possibility), and the method can include one or more other steps which are carried out before any of the defined steps, between two of the defined steps, or after all the defined steps (except where the context excludes that possibility). [0036] The term “at least” followed by a number is used herein to denote the start of a range beginning with that number (which may be a range having an upper limit or no upper limit, depending on the variable being defined). For example, “at least 1” means 1 or more than 1. The term “at most” followed by a number is used herein to denote the end of a range ending with that number (which may be a range having 1 or 0 as its lower limit, or a range having no lower limit, depending upon the variable being defined). For example, “at most 4” means 4 or less than 4, and “at most 40% means 40% or less than 40%. When, in this specification, a range is given as “(a first number) to (a second number)” or “(a first number)-(a second number),” this means a range whose lower limit is the first number and whose upper limit is the second number. For example, 25 to 100 mm means a range whose lower limit is 25 mm, and whose upper limit is 100 mm. [0037] The embodiments set forth the below represent the necessary information to enable those skilled in the art to practice the invention and illustrate the best mode of practicing the invention. For the measurements listed, embodiments including measurements plus or minus the measurement times 5%, 10%, 20%, 50% and 75% are also contemplated. For the recitation of numeric ranges herein, each intervening number there between Attorney Docket No.: NTTR-2023003pct with the same degree of precision is explicitly contemplated. For example, for the range of 6-9, the numbers 7 and 8 are contemplated in addition to 6 and 9, and for the range 6.0-7.0, the number 6.0, 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, and 7.0 are explicitly contemplated. [0038] The term “substantially” means that the property is within 80% of its desired value. In other embodiments, “substantially” means that the property is within 90% of its desired value. In other embodiments, “substantially” means that the property is within 95% of its desired value. In other embodiments, “substantially” means that the property is within 99% of its desired value. For example, the term “substantially complete” means that a process is at least 80% complete, for example. In other embodiments, the term “substantially complete” means that a process is at least 90% complete, for example. In other embodiments, the term “substantially complete” means that a process is at least 95% complete, for example. In other embodiments, the term “substantially complete” means that a process is at least 99% complete, for example. [0039] The term “substantially” includes a value that is within 10% less than or greater than the indicated value. In certain embodiments, the value is within 5% less than or greater than of the indicated value. In certain embodiments, the value is within 2.5% less than or greater than of the indicated value. In certain embodiments, the value is within 1% less than or greater than of the indicated value. In certain embodiments, the value is within 0.5% less than or greater than of the indicated value. [0040] The term “about” includes when value is within 10% of the indicated value. In certain embodiments, the value is within 5% of the indicated value. In certain embodiments, the value is within 2.5% of the indicated value. In certain embodiments, the value is within 1% of the indicated value. In certain embodiments, the value is within 0.5% of the indicated value. [0041] In addition, the invention does not require that all the advantageous features and all the advantages of any of the embodiments need to be incorporated into every embodiment of the invention. [0042] Turning now to Figs.1 - 12, a brief description concerning the various components of the present invention will now be briefly discussed. Attorney Docket No.: NTTR-2023003pct [0043] Experimental neuroscience is advancing rapidly. As this progress continues, it becomes increasingly important to update theoretical and computational tools to keep pace. With the recent advances in targeted stimulation of individual neurons, biologically interpretable computational models with generative power are needed to supply these experiments with real-time computational control. [0044] Among those biologically interpretable models are the data-constrained recurrent neural networks (dRNNs), where the number of units in the artificial network is matched to the number of neurons recorded, and the RNNs are trained to match the observed neural activities. The previous work in this paradigm uses FORCE learning, which was originally developed as an alternative to the reservoir computing paradigm for training read-out weights to produce a target function and then later applied to training recurrent weights to create dynamic attractors. [0045] Data-constrained RNNs have immense potential to generalize the decades old approaches of estimating functional connectivity between brain regions through co- fluctuations that are still being utilized in systems neuroscience to allow reverse engineering underlying computations, to provide real-time computational control for causal experimentation, and to serve as a computational complement to the widely used optogenetics. Influencing biological experiments computationally on the fly requires models that can be trained and provide predictions in real-time. Additionally, biological interpretability of the underlying architecture is necessary for reverse-engineering the underlying computation. These requirements render many traditional machine learning approaches trained in large clusters with billions of parameters impractical for a direct application to experiments performed with limited computational resources, a large portion of which would be allocated to other computational procedures such as real- time extraction of neural signals. [0046] When developing optimization procedures for data-constrained RNNs, another area of particular interest is the scalability of the training to data collected from several days and from more than 10,000 neurons for offline processing. With a model having scalable and fast learning capabilities, the inventors could potentially process multi- day imaging datasets and/or the brain-wide large datasets of the future and test Attorney Docket No.: NTTR-2023003pct biological hypothesis regarding large-scale computations at an instance. Although seven years passed since the original attempt at reproducing neural activities with data- constrained interpretable RNNs, the ability of the data constrained RNN training paradigm to influence biological/experimental progress has been limited if any. Perhaps one major contribution to this is the incompatibility of the FORCE approach with the realities of the experimental data, as well as its extremely low speed of convergence. In fact, the dRNN training has never been truly validated via a comparative study of potential optimizers. Moreover, previous work attempted to solve two problems at once, e.g. designing an optimizer for the data constrained RNN training paradigm and transforming the experimental neural activities to a format trainable by a vanilla recurrent neural network with sigmoid nonlinearity. Yet, both problems are interdisciplinary in nature and require full attention of their own. Thus, for the scope of this work, their focus is explicitly on the former, i.e. the validation of the data-constrained RNN training. The inventors leave the latter problem, e.g. the transformation of the experimentally observed neural activities into a suitable data format, to future work. In this disclosure, the inventors first present a task called “large- scale functional connectivity extraction,” in which the inventors attempt to learn the functional connectivity between “recorded neurons” simulated from a generator model. Then, the inventors introduce a convex and scalable solver, called CoRNN, and validate its accuracy and speed with respect to alternatives including FORCE employed by the seminal works. Contributions include the following: • The inventors introduce the large-scale functional connectivity extraction task and validate reproduction ability of data constrained RNNs on ground truth simulations • The inventors introduce a new convex solver called CoRNN: Scalable, stable, and interpretable convex recurrent neural networks for copying brain dynamics. • The inventors introduce a fixed-point initialization strategy to speed up convergence not only for CoRNN but for FORCE as well. [0047] Unlike the traditional setting of dynamical system reconstruction task, where powerful algorithms (with potentially hidden units and deep layers) are trained to reproduce the underlying dynamical systems of an observation, the data-constrained RNN training Attorney Docket No.: NTTR-2023003pct aims to match the artificial units with the number of observed neurons to interpret currents shared between brain regions and/or population of neurons as proxy for functional connectivity. While ignoring the hidden units may seem too restrictive and confounding when reproducing observed neural activities, several critical biological observations motivate the use of dRNNs: 1) the observed neural activities tend to lie on low dimensional manifolds, 2) the noise correlations limit the maximum amount of information in the population, which can be recovered with few thousands of neurons, and 3) the currents shared between brain regions and/or simulated networks of networks are robust dRNNs as task, a more

of functional connectivity estimation, which was originally performed via considering co-fluctuations between brain regions. [0048] 2. Theory 2.1. Data constrained RNN training for large scale functional connectivity task the inventors start by introducing the large-scale functional connectivity extraction task, as shown in Fig.1, which aims to reproduce interactions between brain regions and/or groups of neurons (e.g. functional connectivity). The inventors assume that the observations of real data can be captured by a generator RNN model. Then, the goal of the dRNNs is to reproduce the internal

the generator RNN. In most general case, observed neurons can be subsampled, or the genera tor might follow a different statistical model than dRNN, in which case the connectivity matrix learned by dRNN would not necessarily coincide with the underlying synaptic connectivity; rather be [0049] For the scope

evidence for the optimization problem; thus, focus on the case

matching generator and learners. In this case, the terms functional and synaptical connectivity coincide, allowing a simple interpretation/scoring of the reproduction quality. The inventors test not only the and control possibility,

In its essence, large scale functional connectivity task is a compromise between co-fluctuations Attorney Docket No.: NTTR-2023003pct looked traditionally through noise and/or canonical correlations between brain regions; and the dynamical system reconstruction task, which traditionally uses more accurate/powerful approaches.

[0050] 2.2. Generator model and estimator family of interest Following the trend in neuroscience, the inventors choose to focus on a recurrent neural network-based generator model. the are assumed to be of the form:

._* = 4⁹ 6_> ^:; ?^< 5_* ⁷ 6⁸¹ )₆ + 4^{9 D3EFG} 6_> ^AB ? 5_* ^* 6^@ C₆ + 0_* H (1) where z_i is the total input to the recurrent unit i, r_i is the firing rate of the unit i, and f (.) is a pre-defined non-linearity, 5_* ⁷ 6⁸¹ is the weight matrix for the recurrent connections, 5_* ^D 6³ is the weight matrix from the input to the recurrent units, u is the input vector, and 0^D3EFG * is the input noise to the recurrent units. Unlike previous literature, given their goal of reconstruction, the inventors also consider the existence of a conversion noise 0_* ¹²³ that accounts for errors and/or miss-matches in the non-linearity. Additionally, any self- excitation of

[0051] Suitable a can (recursively) by discretizing the

[0052] where the time scale parameter is defined as I' = ^J& K . Here, the model parameters influence .L lin

i,t early such that

Attorney Docket No.: NTTR-2023003pct where we define the short-hand notations x = [r; u] and M^N = [W^rec; Wⁱⁿ]^T as concatenated matrices. The goal of the estimation problem is to find M^N such that M^N ≈ θ. [0053] 2.3 Convexification problem, a reasonable choice

However, time recurrent definition of the general estimator )L t,i, this

minimization problem would require backpropagation through time. This approach would miss one of the inventors’ main able to scale well to tens of [0054] during training

turns the potentially infinite time problem into a single time-step one, which are called “single step prediction error paradigm”. Under this condition, the L2 error transforms to a simpler form:

$&O?H*>-?>P/$&H* where we define dt,i = _P and correspondingly Q^R t,i = f (.L) becomes an However, this global loss in Fig. 2. ?OS [0055] that ^V _THU would

_W follow a linear + sigmoid form, the inventors replace the L2 loss with a cross-entropy loss that is well known to be convex under this estimator, see logistic regression. On the other of a naive replacement, the inventors use a weighted loss of the

form Attorney Docket No.: NTTR-2023003pct Here, the

factors chosen as c_t,i = [X ( Q_& ^W H_* ]^-1 for the CoRNN solver for scalability as discussed below. For now, the regularized general loss function is written as [0056] loss function, such that any given i, the (vector, not

matrix) parameters := (i) = M^N _ji and c_t := c_t,i can be considered separately (See Appendix A) The gradient and Hessian of this smaller problem can be found as:

[0057] In this naive form, the Hessian would need to be computed separately for each subproblem and after each parameter update as it depends on Q^R t. This motivates the choice of c_t = [1−Q_& ^W]⁻¹, since near the solution Q^R _t ≈ d_t, the Hessian approximates to:

Attorney Docket No.: NTTR-2023003pct which is both parameter Z^R (and thus M^N) and target d_t independent; this is commonly shared across all subproblems as long as the convex program starts and remains nearby the optimum. Owing to its fixed nature, this Hessian is defined as the fixed-point Hessian and use this globally to perform optimization updates. Finally, the gradient of the full problem can be obtained by promoting Eq. 8a to a matrix: where we define the prediction error matrix E. It is worth noting that each column of the gradient matrix contains the gradient for the scalar subproblem corresponding to the ith output. Then, defining the inverse (fixed-point) Hessian matrix as we obtain the following update rule (See Appendix A) [0058] where all quantities that can be pre-computing the prediction error followed by matrix multiplications are underlined; this is well suited to be accelerated on a graphics processing unit (GPU). [0059] 2.5 Choosing a suitable fixed point. The approximation performed at the local Hessian ?>S_T ^{\^} stage, e.g. _ X ^{(0) *} ?>S_T^ , requires that we initialize the parameters θ := θ_fp close to _ M , where M^* is the solution to the optimization problem. If this condition is met, the Attorney Docket No.: NTTR-2023003pct approximate local Hessian derived in Eq. 9 would describe the local topology near the solution, the fixed point, and the steps in-between. As the inventors show in Appendix A, a suitable fixed-point candidate is the solution to the least squares problem on the currents Z:

This step is computationally negligible, as the underlined matrix is already pre- computed for the iterations that follow. [0060] 3. Results. 3.1 Convex recurrent neural networks (CoRNN) via ADMM framework. When the inventors derived the update equations of the fast solver in the previous section, the inventors omitted the constraints 5_* ⁷ *⁸¹ = 0. However, these constraints could be biologically relevant, and perhaps be vital to eliminate overfitting in the reproduced network. In fact, for the generator model given in Eq.1, r_t,i ≈ r_t+1,i for small α; making the unconstrained problem ill-posed since self-excitation would explain a good amount of the activity of each neuron. Thus, to incorporate the equality constraint, the inventors developed a solver using the alternating direction method of multipliers (ADMM), summarized in Fig.2. [0061] In the ADMM solver, the primal problem is divided into two subproblems that are solved individually during separate iterations, yet linked through a consensus variable that ensures agreement between the solutions of both subproblems upon convergence. In this framework, first subproblem solves the unconstrained optimization problem whereas the second one enforces the constraints (See Appendix A for derivations). While the inventors primarily consider the equality constraint in this work, ADMM framework can be used to add additional L1 and/or nuclear norm regularizations to the learned connectivity matrix, and utilize biological priors for further inductive biases. Finally, the inventors note that given the highly irregular nature of the prediction error when dt,i ≈ ±1, the inventors supply the fast solver with a simple automated outlier detection step explained in Appendix A. [0062] 3.2 Chaotic recurrent neural network reproduction. To validate CoRNN’s reproduction abilities, the inventors utilize randomly connected chaotic recurrent neural networks Attorney Docket No.: NTTR-2023003pct without any external input as generators. Specifically, the inventors pick the recurrent weights randomly from a normal distribution 5_* ⁷ *⁸¹! N `aH ^b c_@:;<d with zero mean and b standard deviation of _c@:;<. For e f X and sufficiently large g₇₈₁, these networks can sustain their activities without any external stimulus and exhibit chaotic behavior. For the chaotic generators, the inventors pick e = h throughout this disclosure to ensure that network is deep into the sustainable spontaneous activity regime. The inventors pick the chaotic regime specifically, since previous work established that teacher forcing when training RNN to output a function can have potentially harmful consequences to the stability of the network training. The inventors aim to test whether this generalizes to the reproduction task considered here. [0063] To validate that even chaotic dynamics can be reconstructed via a single step prediction error paradigm, the inventors initialized randomly connected recurrent neural networks as generators, extracted neural activities for a fixed duration of length T, and trained the dRNNs using CoRNN, Newton’s solver, gradient descent, and FORCE on varying number of iterations to benchmark their speed and accuracy. All algorithms were able to train, some achieving near perfect accuracies in reproducing the internal connectivity of the generator network, as shown in Fig. 3. [0064] Whether minimizing weighted or the naive logistic loss leads to more accurate reconstruction depended on the underlying noise distribution of the data as shown in Fig. 3, yet CoRNN outperformed all other algorithms in speed in both cases with several orders of magnitude and converged in sub-second times even with the CPU implementation, whereas others took up to several seconds for training. In fact, FORCE was not able to converge in a pre-allocated maximum training time, showing the inefficiency of FORCE training for the reproduction task. To compare, CoRNN converged faster at a given run than the computer could import the Pytorch module at the start of the code.! [0065] 3.3. Robustness to regularization Having tested the speed of the algorithms, the inventors next turn their attention to robustness to regularization. If the instability claims regarding network training in generalizes to the network reproduction problem, Attorney Docket No.: NTTR-2023003pct then the single-step prediction error paradigm, regardless of the optimizer, would only be successful for a narrow range of regularization parameters. [0066] To test this scenario, the inventors first simulated a set of networks with parameters matching those in Fig. 3A and trained using various optimizers for a given number of iterations and varying (normalized) regularization levels. The results are shown in Fig. 4A, with regularization level = 1 corresponding to the default regularization parameter of the respective algorithm. Unlike expectations, all single-step prediction error solvers showed high robustness to changes in the regularization level for several orders of magnitude, whereas FORCE either only maintained the accuracy achieved already the fixed-point initialization or outperformed the initialization only for a hyper-tuned regularization value and led to a completely incorrect reproduction for others. Then, the inventors tested whether the fixed-point initialization can account for the high accuracies of the single-step prediction error solvers. Specifically, the inventors initialized the networks to a random initialization instead of the fixed point and re-did the analysis of Fig. 4A in Fig. 4B. The results indicate that fixed point initialization is not necessary to obtain high accuracies, although it has two distinct advantages: i) fixed point initialization leads to a slightly higher accuracy, which the inventors test more explicitly next; and ii) it stabilizes the FORCE solver for the dRNN training. In Fig. 4B, the inventors still initialize CoRNN to the fixed point to have a reference-points between two subfigures, although CoRNN can learn without fixed-point initialization as discussed next. [0067] 3.4. Fixed point initialization increases convergence speed. Inspired by their observation in Fig. 4B that fixed point initialization increases the accuracies slightly, the inventors next tested whether the asymptotic accuracy is increased with the better initialization, or whether the convergence speed is slowed down with the random initialization compared to fixed point initialization. To test this scenario, the inventors simulate a slightly larger network of 500 neurons, which would amplify any difference in convergence speeds compared to a smaller network, for T = 10000 data points and re-ran the algorithms without fixed-point initialization. As shown in Fig. 5, CoRNN with or without fixed point initialization eventually converges to the same asymptotic Attorney Docket No.: NTTR-2023003pct accuracy, while the fixed-point initialization increases the speed of convergence at least an order of magnitude for CoRNN. [0068] 3.5. 3-bit flip flop task Having provided validity evidence for the speed and accuracy of CoRNN, or single-step prediction error solvers in general, via the randomly connected chaotic RNN generators; the inventors next focus on an RNN with a more structured connectivity. Specifically, the inventors train a model RNN to perform a 3- bit flip flop task following as shown in Fig. 6. [0069] The 3-bit flip flop task had 3 dimensions of input and 3 of output. Each trial lasted 100 timepoints. Input dimensions could take on the values of 0, -1, or 1. At every timepoint, and for each dimension independently, there was a p = 0.05 probability of being in a - 1 or 1 input state with equal probability, otherwise the state was 0. The output dimensions were initialized to 0, and would “flip” to -1 or 1 depending on the input in the same dimension, and stay in that state until the input dimension took on the other value. For instance, if the first input dimension was -1 at a certain time, then the output in that dimension would be -1 from that time forward, until the first input dimension took on a value of 1. At that point the first output dimension would take on a value of 1, until the first input dimension took on a value of -1. Thus, values in the input dimensions “flip” the state of the output dimensions, like a switch. [0070] Fig. 6A shows the output of an example test trial, not used for the generator RNN training or the CoRNN reproduction, with 200 time points followed by 100 time points of intertrial interval. Not only CoRNN reproduced network can perform the task, it gives the same output as the original network; potentially making similar mistakes as the network. Thus, as a first step, the inventors confirm that CoRNN can learn to output the same outputs as the original network. [0071] Next, the inventors look at underlying neural activities of an example (original) generator and reproduced (CoRNN learned) network in Fig. 6B. The reproduction matches not only the output of the network, but also the internal firing dynamics of individual units. Yet, as shown in Fig. 6B, the synaptic connections are not perfectly learned even in such a simple scenario with matching models, which requires more trials and potentially interventional data as explored next. Although surprisingly, the Attorney Docket No.: NTTR-2023003pct input synaptic connections are perfectly learned in this scenario as shown in Fig.6C, a potentially useful observation for causal biological experiments. 3.6. More trials may be necessary for larger networks Inspired by their observation that firing rates may be perfectly predicted even in test trials for large networks without exactly reproducing the underlying synaptic connectivity, the inventors next tested whether increasing the number of trials could lead to better reproduction of the original network. [0072] Fig.7 shows that the output of the network can be reproduced with far smaller numbers of trials, whereas accurately reproducing the internal connectivity requires significantly more trials. The data hungriness of the estimation process increases with the number of units in the generator network, which again emphasizes the importance of an optimization algorithm that can scale well to multiple day recordings. 3.7. CoRNN scales well to larger networks and recordings So far, the inventors showed that CoRNN works faster than alternative algorithms with run times at most in few minutes. Now the inventors show that, owing to the first order nature of the updates, CoRNN also scales well to larger networks and data points as shown in Fig.8. Remarkably, CoRNN remains sub-minute for 200 trials (approximating an hour-long recording) and 1000 neurons, meaning that CoRNN can indeed be trained real-time even in batches in biological experiments. [0073] 4. Discussion. In this disclosure, the inventors introduced a convex and scalable solver (CoRNN) for training data-constrained recurrent neural networks to reproduce brain activity. The inventors validated CoRNN on simulated data via the large-scale functional connectivity extraction task, where the inventors used both randomly connected chaotic networks and RNNs trained to perform a 3-bit flip flop task as generators. The inventors showed that CoRNN is as accurate yet faster than existing approaches by orders of magnitude and can easily scale to the large datasets of today, training in seconds to minutes depending on the number of trials and neurons. For this work, the inventors focused on developing a scalable, fast, and accurate solver for the large-scale functional connectivity extraction task. Thus, the inventors conceive several practical aspects such as subsampling of neuronal populations, miss-matches between non-linearities, miss-matches in time scales, and bringing experimental data in a Attorney Docket No.: NTTR-2023003pct suitable format for the dRNN training. Further development of data constrained RNNs can support systems neuroscience research aiming to understand interarea interactions and replace the canonical/pairwise correlations as a proxy for functional connectivity estimation. Moreover, reverse-engineering the learned network can provide deeper understanding of how neuronal populations compute. Finally, with the advent of interventional data collection methods such as single-cell targeted optogenetics, dRNNs can be trained to test causal connections computationally and supply theoretical predictions at a scale previously unachievable for experimentation. [0074] Fig. 9 is a functional block diagram illustrating an example computing device 100 configured to train a neural network to interpret brain dynamics, incorporated with teachings of the present disclosure, according to some embodiments. Fig. 9 illustrates only one particular example of computing device 100, and many other examples of computing device 100 may be used in other instances and may include a subset of the components included in the example computing device 100. Computing device 100 may include additional components not shown in FIG. 9. [0075] As shown in the example of FIG. 9, computing device 100 includes one or more processing devices 102, one or more input devices 104, one or more output devices 106, one or more communication devices 108, one or more storage devices 110, network interface 112, and memory devices 114. Processing devices 102, input devices 104, output devices 106, communication devices 108, storage devices 110, network interface 112, and memory devices 114 are interconnected via bus 190. [0076] Processing devices 102 may include one or more execution cores (CPUs). For example, computing device 100 may also include a peripheral controller hub (PCH) (not shown). In another example, computing device 100 may also include a sensors hub (not shown). Input devices 104 and output devices 106 may include, for example, user interface device(s) including a display, a touch-screen display, printer, keypad, keyboard, etc., sensor(s) including accelerometer, global positioning system (GPS), gyroscope, etc., communication logic, wired and/or wireless, storage device(s) including hard disk drives, solid-state drives, removable storage media, etc. I/O ports for input devices 104 and output devices 106 may be configured to transmit and/or receive commands and/or Attorney Docket No.: NTTR-2023003pct data according to one or more communications protocols. For example, one or more of the I/O ports may comply and/or be compatible with a universal serial bus (USB) protocol, peripheral component interconnect (PCI) protocol (e.g., PCI express (PCIe)), or the like. [0077] Processing devices 102 may implement functionality and/or execute instructions within computing device 100. The instructions may be stored completely or partially within storage devices 110 and memory devices 114. For example, processing devices 102 on computing device 100 may receive and execute instructions stored by storage devices 110 that provide the functionality of training a neural network to interpret brain dynamics. These instructions executed by processing devices 102 may cause computing device 100 to store and/or modify information within storage devices 48 during program execution. [0078] Storage devices 110 may generally comprise a random-access memory (“RAM”), a read-only memory (“ROM”), and a permanent mass storage device, such as a disk drive or SDRAM (synchronous dynamic random-access memory). Computing device 100 may store program code for modules and/or software routines. [0079] Storage devices 110 may also store operating system, kernel, kernel space, user space, and user protected address space. Storage devices 110 may store one or more processes (i.e., executing software application(s)). One or more processes may generally execute in parallel, i.e., as a plurality of processes and/or a plurality of threads. These software components may be loaded from a non-transient computer-readable storage medium into storage devices 110 using a drive mechanism associated with a non-transient computer-readable storage medium, such as a floppy disc, tape, DVD/CD-ROM drive, memory card, or other like a storage medium. In some embodiments, software components may also or instead be loaded via a mechanism other than a drive mechanism and computer-readable storage medium (e.g., via network interface 112). [0080] Embodiments of the presently disclosed invention are a computer software system and a solution to an important optimization problem. Some embodiments include convex optimization method that fits a recurrent neural network to reproduce the neuronal activities of mice or other animals observed from experimental data through calcium Attorney Docket No.: NTTR-2023003pct and voltage imaging. Some embodiments could fall within the domain of “listening to” the neural activities, but beyond listening, embodiments model the activities. Moreover, some embodiments have practical implications for brain machine interface research. [0081] A distinguishing factor of some embodiments of the presently disclosed invention to previous technology is the formulation of the problem in a language that matches the biological reality of experiments. FORCE, the previous approach, has several short- comings that is described herein. Moreover, the inventors develop the convex solver that is 1-2 order of magnitude faster than the closest competition. [0082] Among many possible uses or applications of some embodiments of the presently disclosed invention are brain-machine interfaces using this approach, research regarding large-scale connectivity/discovery of functionally connected brain circuits from calcium imaging recordings etc. Embodiments of the presently disclosed invention may be used with efferent BMIs, passive BMIs, sensory BMIs, two-photon (2P) calcium imaging, genetically encoded calcium indicators (GECOs), for example. [0083] The invention illustratively disclosed herein suitably may explicitly be practiced in the absence of any element which is not specifically disclosed herein. While various embodiments of the present invention have been described in detail, it is apparent that various modifications and alterations of those embodiments will occur to and be readily apparent those skilled in the art. However, it is to be expressly understood that such modifications and alterations are within the scope and spirit of the present invention, as set forth in the appended claims. Further, the invention(s) described herein is capable of other embodiments and of being practiced or of being carried out in various other related ways. The present disclosure also contemplates other embodiments “comprising,” “consisting of” and “consisting essentially of,” the embodiments or elements presented herein, whether explicitly set forth or not. In addition, it is to be understood that the phraseology and terminology used herein is for the purpose of description and should not be regarded as limiting. The use of “including,” “comprising,” or “having” and variations thereof herein is meant to encompass the Attorney Docket No.: NTTR-2023003pct items listed thereafter and equivalents thereof as well as additional items, while only the terms “consisting of” and “consisting only of” are to be construed in the limitative sense.

APPENDIX A CORNN: Convex optimization of recurrent neural networks for rapid inference of neural dynamics Fatih Dinc^∗ Adam Shai^∗ Mark J. Schnitzer^† Department of Applied Physics CNC Program Howard Hughes Medical Institute Stanford University Stanford University CNC Program Stanford, CA 94305 Stanford, CA 94305 Stanford University Stanford, CA 94305 Hidenori Tanaka^† Physics & Informatics Laboratories, NTT Research, Inc. Sunnyvale, CA 94085 Center for Brain Science, Harvard University Cambridge, MA 02138 Abstract Advances in optical and electrophysiological recording technologies have made it possible to record the dynamics of thousands of neurons, opening up new possibili- ties for interpreting and controlling large neural populations in behaving animals. A promising way to extract computational principles from these large datasets is to train data-constrained recurrent neural networks (dRNNs). Performing this training in real-time could open doors for research techniques and medical applications to model and control interventions at single-cell resolution and drive desired forms of animal behavior. However, existing training algorithms for dRNNs are inefficient and have limited scalability, making it a challenge to analyze large neural record- ings even in offline scenarios. To address these issues, we introduce a training method termed Convex Optimization of Recurrent Neural Networks (CORNN)¹. In studies of simulated recordings, CORNN attained training speeds ∼100-fold faster than traditional optimization approaches while maintaining or enhancing modeling accuracy. We further validated CORNN on simulations with thousands of cells that performed simple computations such as those of a 3-bit flip-flop or the execution of a timed response. Finally, we showed that CORNN can robustly reproduce network dynamics and underlying attractor structures despite mismatches between generator and inference models, severe subsampling of observed neurons, or mismatches in neural time-scales. Overall, by training dRNNs with millions of parameters in subminute processing times on a standard computer, CORNN constitutes a first step towards real-time network reproduction constrained on large-scale neural recordings and a powerful computational tool for advancing the understanding of neural computation. 1 Introduction Understanding the relationship between neural dynamics and computational function is fundamental to neuroscience research [1]. To infer computational structure from neural population dynamics, ∗These authors contributed equally to this work. †These authors co-supervised this work. 1CORNN software and reproduction code are available at https://github.com/schnitzer-lab/CORNN-public 37th Conference on Neural Information Processing Systems (NeurIPS 2023). APPENDIX A neuroscientists regularly collect and analyze large-scale neural recordings using electrophysiological or optical imaging recording methods. Both recording modalities have undergone rapid progress in recent years, yielding a steady increase in the numbers of cells that can be recorded and manipulated and trans- way of and of indi- valuable, To on

neural network-based modeling, is the training of dRNNs [24–29]. These models, part of a broader suite of computational strategies [11, 21, 30–32], align the dynamics of the neural units within the RNN to the activities of recorded neurons. This data-constrained methodology for RNN training has the potential to modernize estimations of functional connectivity [33]—still a mainstay in systems neuroscience [6]—and provides a viable path to distill underlying computations by extracting the state-space representations of the learned dRNNs [12]. Additionally, dRNNs may open the door to computational control for causal experimentation [9, 10], serving as a valuable computational adjunct to widely used techniques like optogenetics [8, 34]. Despite recent advances in the dRNN framework, how to perform fast and scalable reconstructions of neural activity traces from large-scale empirical recordings has remained unclear. Notably, the slowness of existing dRNN optimization algorithms may often necessitate the use of high-performance clusters and several days of computation in large-scale experiments. These limitations might have been a barrier to the widespread adoption of dRNN approaches. Moreover, to fully harness the potential of dRNNs, it is essential to extend their range of applicability from offline analyses to on-the-fly applications within individual recording sessions. An approach to training dRNNs in a nearly immediate manner would be a key advancement that would enable novel experimental approaches, such as theory-driven real-time interventions targeting individual cells with specific computational or functional roles (Fig. 1). However, realizing these benefits requires having an optimization routine that is fast, robust, and scalable to large networks. In this work, we present CORNN, a convex solver designed for fast and scalable network reproduction. We demonstrate its accuracy and efficiency compared to alternative network training methods such as FORCE [35] and back-propagation through time (BPTT), which have been employed in previous seminal works utilizing dRNNs [24–29]. Our main contributions include the development of CORNN (Fig. 2), the introduction of an initialization strategy to accelerate convergence (Figs. 2, S4, and S6), and the demonstration of CORNN’s scalable speed (Figs. 3, 5, S2, S3, S5, S8, and S9), which enables rapid training on the activity patterns of thousands of neurons (Figs. 4, 5, S7, S8, S9, S10, and S11). CORNN’s performance can be further enhanced by 1-2 orders of magnitude through the use of a standard graphical processing units (GPU) on a desktop computer (Figs. S5 and S8). Unlike BPTT and FORCE, CORNN does not require fine-tuning (Figs. S1, S6), making it a user-friendly technology for biologists. Lastly, we highlight CORNN’s robustness against non-idealities such as mismatches in the assumed dynamical system equations (Figs. 4, 5, S10, and S11), subsampling of neural populations (Fig. 5), differences in dynamical time-scales (Fig. S10), and existence of non-Gaussian (Fig. S3) or correlated noise (Fig. S11). By enabling user-friendly, subminute training on standard computers, CORNN represents a first necessary step in transforming data-constrained recurrent neural networks from a theoretical concept into an experimental/computational technology. While this work focuses on introducing and validating the fast solver on simulated benchmarks, future work with CORNN should focus on addressing the challenges that arise when analyzing large-scale neural recordings from experimental neuroscience. APPENDIX A Neural recordings Real-time interventions to neural populations Real-time extraction Behavioral experiment

RNN Desired behavior Simulate single-cell targeting Figure 1: Using data-constrained recurrent neural networks for the interpretation and manipulation of brain dynamics within a potential experimental pipeline. This approach centers around online modeling of network dynamics, which can enhance hypothesis testing at the single-cell level and support advancements in brain-machine interface research. The training process is motivated by three objectives: (i) predicting the patterns of neural populations, (ii) revealing inherent attractor structures, and (iii) formulating optimal control strategies for subsequent interventions. 2 Approach 2.1 Experimental setup for real-time interventions using dRNNs Previous work suggested the use of dRNNs for real-time feedback between experimental and com- putational research [8]. To date, however, the use of dRNNs has remained offline. Past studies extracted computational principles from neural recordings by fitting dRNNs after data collection [24–29]. Moreover, the use of dRNNs has generally required the expertise of computational neu- roscientists, likely due to the complications associated with the optimization and training of neural network models. However, with the advent of large-scale recordings and data-driven approaches, real-time brain-machine interface research is likely a forthcoming application [36–38]. In such use cases, neurobiologists might wish to train dRNNs via a user-friendly approach. Hence, a fast and straightforward dRNN training procedure would enable a new breed of interventional experiments to dissect the brain’s microcircuitry. To appreciate the advantages of fast data-driven reconstruction of neural dynamics, consider a hypothetical scenario in which dRNNs, facilitated by CORNN, enable real-time interventions at the single-cell level . In this neural activities from mice behavioral tasks are and pre-

from each mouse are reproduced by training dRNNs, which are then reverse-engineered to reveal underlying attractor structures. Rapid dRNN training allows for a tight feedback loop between incoming measured data and optimal experimental design, in order to refine the inferred model. Techniques similar to adversarial attacks [42] might be used to devise optimal cell targeting strategies from these dRNNs [8], allowing one to test hypotheses about the computational roles of individual neurons or to identify optimal neurons for use within brain-machine interfaces. Once a perturbation strategy is determined, it can be tested on subsequent trials of the experiment, using the same animal whose recorded neural dynamics led to the trained dRNN. Thus, a fast dRNN training algorithm could allow for better fitting of a dynamical system model to the experimental data and provide a natural testbed to probe hypotheses about the computational structure of the biological neural circuitry. Motivated by the goals described our focuses on dRNNs and as fast as possible. However,

constraints. Firstly, to ensure real-time

with the experimental apparatus, we require that the training process takes place on computers, not clusters. the computational complexity of real-time processing in large scale recordings,

of neurons, we expect that at least one of the GPUs is reserved for APPENDIX A Fixed point initialization & Activities Hessian inverse computation Generator Optimization loop Convex loss Convex RNN

Figure 2: CORNN: Convex and scalable solver for dRNNs via ADMM [51]. The CORNN algorithm optimizes the parameters of a recurrent neural network so that activity in hidden units align with activities measured from a ground-truth system, we refer to as a generator. The choice of objective function results in a convex loss landscape. The algorithm starts by finding a fixed point for initialization, where the Hessians of all subproblems are aligned to the correlation matrix of neural activities, which can be pre- predicting the neural checking for convergence. The the functional connections

from brain- [43] see [44] for a development in smaller scale experiments), or spikes from electrophysiological recordings [7], leaving the central processing unit (CPU), or perhaps a second GPU, for training dRNNs. Finally, maintaining biological interpretability via 1-1 matching of observed and modeled neurons is crucial, as arbitrariness brought by hidden units can explain away existing functional connectivity and/or contributions of observed neurons to underlying attractor structures. These concerns make many

learning architectures and optimization routines, typically trained on large clusters with billions of parameters, unsuitable for our purposes. This includes recent work on the dynamical system reconstruction paradigm [31, 45–50], which chiefly aims to maximize reconstruction accuracy of dynamical systems without regard to biological and 2.2 Inference

Recurrent neural networks are universal approximators for dynamical systems, whose internal computations can be reverse engineered and interpreted [12, 52]. Thus, we choose leaky firing-rate RNNs as the inference models of CORNN, which follow the dynamical system equations [22]: dr_i ^N∑ ∑ ^N τ = −r_i + n ^{rec in} (1) dt

where τ is the neural decay time, z_i is the total input to the recurrent unit i, r_i is the firing rate of the unit i, and tanh(.) is the pre-defined non-linearity, W_i ^r j^ec is the weight matrix for the recurrent connections, W_i ⁱ jⁿ is the weight matrix from the input to the recurrent units, u_j is the input vector, and ǫⁱ i^nput is the input noise to the recurrent units. Unlike previous literature, given our goal of reconstruction, we also consider the existence of a conversion noise ǫ^c i^on that accounts for errors and mismatches in the non-linearity. Additionally, the inference model has the constraint W_i ^r i^ec = 0 preventing any self-excitation of neurons. Suitable estimators for reconstructing such a model can be obtained (recursively) by discretizing the differential equation: rˆ_t+1,i = (1− α)rˆ_t,i + αf(zˆ_t,i), (2) where we define the time scale ^∆t zˆ_i,t explicitly linearly, and implicitly

We define the short-hand notations x = [r;u] and θˆ = [W ^rec;W ⁱⁿ]^T as concatenated matrices. The goal of the estimation problem is to find θˆ minimizing a loss function. APPENDIX A A reasonable choice for the loss function of this regression problem is the traditional L₂ loss: n∑∑ ^T L₂(θˆ) = (rˆ_t,i − r_t,i)². (3) i=1 t=1 However, due to the time recurrent definition of the general estimator rˆ_t,i, this minimization problem would require backpropagation through time (BPTT). It is not apriori clear if BPTT is the right optimization method for network reproduction, for which one has access to high information content regarding the internal in Instead of using BPTT, we

is teacher-forced to r_t,i when a single time-step one, which we call “single paradigm.” Under this condition, the L₂ loss function transforms to a simpler

n∑∑ ^T L₂(θˆ) = α2 (dˆ_t,i − d_t,i)², (4) i=1 t=1 where we define d_t,i = r −(1−α)r and correspondingly dˆ_t,i = tanh(zˆ_t,i) becomes an estimator consisting of a linear term

is not convex and leads to a 2). 2.3 Convexification of the loss function Rather than minimizing the naive L 1+dˆ 2 loss function, we observe that ₂ follows a linear + sigmoid form. As a result, we replace the L₂ loss function with a cross-entropy loss function that is well known to be convex under this estimator. In fact, in the limit of d_t,i → ±1, the problem reduces to logistic regression. Instead of a naive replacement, we use a weighted loss of the form: )

1 + dˆ L _t,i 1 + d_t,i _weighted(θˆ) = , . (5) i=1 t=1 2

Here, we use CE to denote the cross-entropy loss function and c_t,i are non-negative weighting factors chosen as c_t,i = [1 − d2 t_,i]⁻¹ for the CORNN solver. The specific choice of c_t,i, which resembles (but is not) preconditioning [53], is motivated from a theoretical view in Supplementary Section S1.3. The reader can verify that the Hessian of this loss, derived in Eq. (S4), is positive semi-definite. For now, we write down the regularized loss function for CORNN as: 1 λ L_CORNN(θˆ) = L (θˆ) + ||θˆ||² , (6) T _weighted 2 _F where ||θˆ||_F is the Frobenius norm of θˆ and λ is the regularization parameter. 2.4 Fundamental principle of CORNN: Subproblem Hessians align to the correlation matrix of neural activity traces Having convexified the los ∑s function, we next observe that the loss function itself is perfectly separable such that L_CORNN(θˆ) = _i L_i(θˆ_:,i), where θˆ_:,i denotes the ith column of the parameter matrix θˆ (See Eq. S1). This means that the original problem can be divided into n_rec sub-problems and solved independently. However, inspired by the fact that solving least-squares for a vector or a matrix target has the same complexity, we can devise a faster method. Specifically, the choice of the specific weighting factors, c_t,i = [1− d2 t_,i]⁻¹, leads to a shared (approximate) Hessian for all subproblems (See Eq. (S4)): 1 H ≈ _X ^T _{X + λI,} T ⁽⁷⁾ where X stands for x_t,i in the matrix form, I is the identity matrix, and λ is the regularization parameter. We note that the approximation is exact in the limit dˆ_t,i = d_t,i. In other words, we align the Hessian of each subproblem to the correlation matrix of neural activities. Unlike for APPENDIX A least-squares minimization, this is a local, not global, Hessian; but as long as we can initialize the network sufficiently close to the optimal solution, we expect that the descent direction with the approximate Hessian converges quickly. In this sense, the approximation in Eq. (7) shows similarities with quasi-Newton methods [53] but is distinct in the pre-computed nature of the Hessian. (⁸⁾ of (9)

we θ^k+1 = A⁺X^TXθ^k +A⁺X^TE^k, (10) where we highlight all quantities that can be pre-computed with blue. In essence, the update rule, which we call “Hessian aligned prediction error (HAPE) update,” consists of computing the prediction error followed by matrix multiplications; thus is well suited to be accelerated on a GPU if available, but still can be efficiently computed on CPU cores. The alignment of subproblem Hessians relies on the assumption that we can find a “close-enough” starting position, called the fixed-point, for the initial parameters θˆ⁽⁰⁾. As we show in Supplementary Section S1.2, a suitable candidate is the approximate least-squares solution: θ_ls := A⁺X^TZ. (11) ∑ Then, the initial set of predictions becomes dˆ_t,i := tanh( _j x_t,j(θ_ls)_j,i), which is subsequently refined through the HAPE iterations. This step is computationally negligible, as the blue colored matrix is already pre-computed for the iterations that follow. In this section, we focused on the leaky firing rate RNNs as described in Eq. (1). However, CORNN can be applied to the other widely used variant of RNNs, i.e., the leaky current RNN described in Eq. (S28) and regularly employed in neuroscience literature [12, 26, 52]. The reproduction of the leaky current RNNs can be made convex by observing the direct link between the firing rates, r, and the currents, z, via a linear plus non-linear relationship, i.e., r = tanh(z) and subsequently replacing d_t,i with r_t,i in Hessian and gradient calculations. To prevent introducing unnecessary complexity, we focus on leaky firing rate RNNs in this work. 3 Results 3.1 CORNN as a versatile base model While Hessian alignment, and subsequent HAPE updates, are the core principles driving the fast solver of CORNN, a versatile base model should be able to incorporate various regularization schemes and sparsity constraints. For example, the addition of constraints W_i ^r i^ec = 0 or a potential low-rank regularization would prevent the Hessian alignment and subsequently the use of fast HAPE updates. However, these constraints could be biologically relevant and perhaps vital to eliminate overfitting in the reproduced network. In fact, for the inference model given in Eq. (1), r_t,i ≈ r_t+1,i for small α; making the unconstrained problem ill-posed since self excitation would explain a good amount of the activity of each neuron. Thus, to incorporate the equality constraint and many others, we developed a solver using the alternating direction method of multipliers (ADMM) [51] (summarized in Fig. 2). In the ADMM solver, the primary problem is divided into two subproblems that are solved indi- vidually during subsequent iterations, yet linked through a consensus variable that ensures agree- ment between the solutions of both subproblems upon convergence. In this framework, the first subproblem solves the unconstrained optimization problem, whereas the second one enforces the constraints/regularization (See Supplementary Section S1 for derivations). This division allows the APPENDIX A Figure 3: CORNN reduces training times by several orders of magnitude. The plots il- CORNN GD on L2 loss GD on logistic loss lustrate the relationship between reconstruction FORCE CPU GPU accuracy (Pearson’s correlation coefficient be- 1 min 1 hr 1 day 1 tween ground truth and inferred weights) and n W training time, measured in seconds on a log oⁱ _{t f} cu ^o ^r _t ⁾ _r scale. The FORCE approach (here, on firing s ⁽ rates) is the default method in neuroscience lit- n y o c c erature for dRNN training [24, 25, 28, 29 e ^a _r ]. Pa- u R cc rameters: α = 0.1, n_rec = 5000, T = 30000. a 0 No input. ǫ^conv ∼ Poisson(10⁻³), ǫ^input ∼ 1⁰ ₁₀ ¹⁰ N (0, 10⁻⁴). Lines: median. Error bars: s.d. Training time (s) over 7 networks. use of fast HAPE updates in the computationally extensive first subproblem, where Hessians can be aligned. While we only considered the equality constraint in this work, the ADMM framework can be used to add additional L1 and/or nuclear norm regularizations to the learned connectivity matrix [51, 54], and utilize biological priors for further inductive biases. Finally, we note that given the highly irregular nature of the prediction error when d_t,i ≈ ±1, we supplied the fast solver with a simple automated outlier detection step explained in Supplementary Section S1.5. 3.2 Backpropagation through time on the L₂ loss is suboptimal In a traditional setting, in which neural activations are hidden and need to be learned through training to provide a correct output, backpropagation through time (BPTT) is necessary to learn long-term dependencies in the data. However, in the network reproduction where neural activations are no longer hidden, it is not clear whether BPTT would be beneficial for the learning, since long-term dependencies are inherently present in the state space of the dynamical system, whose equations are given in Eq. (1). Specifically, given that Eq. (1) has only a single time derivative, the time evolution is performed in a Markovian manner, i.e., given the current state, the next state can be computed without need for additional information (up to a corruption by random noise). To (See the with the cross- 2 that L₂ minute even orders of 3.3 that

we as extracted neural activities for a fixed duration of length T , and trained the dRNNs using CORNN, Newton’s solver, gradient descent, and FORCE on varying number of iterations to benchmark their speed and accuracy (See Supplementary Sections S2 and S3.1 for details). All algorithms were able to train, some achieving near perfect accuracies in reproducing the internal connectivity of the generator network (See Figs. 3, S2, and S3). However, FORCE, the current default method in neuroscience literature for dRNN training [24, 25, 28, 29], was 4 orders of magnitude slower (Fig. 3) and needed fine-tuned hyperparameters to barely outperform the fixed-point initialization (Fig. S4). Whether minimizing the weighted or logistic loss led to more accurate reconstruction depended on the underlying noise distribution (we tested two noise distributions: Poisson(mean) and N (mean, variance)) of the data as shown in Figs. 3, S2, and S3. However, CORNN outper- formed all other algorithms in speed in both cases with several orders of magnitude and converged APPENDIX A in sub-second times even with the CPU implementation, whereas others took up to several days for training (Fig. 3). We also observed increased efficiency for CORNN on a GPU (Fig. S5) and that fixed-point initialization could stabilize other algorithms (Figs. S4 and S6). 3.4 CORNN runtimes scale linearly with data size and polynomially with the network size To this point, we provided validity evidence for the speed and accuracy of CORNN, or single-step prediction error paradigm in general, via the randomly connected chaotic RNN generators. Next, we focused on RNNs with more structured connectivity matrices. Specifically, we trained a model RNN using BPTT to perform a 3-bit flip flop task following [12] as shown in Fig. S7 (See Supplementary Section S3.2 for details). Fig. S7A shows the output of an example test trial, not used for the generator RNN training or the CORNN reproduction, with 200 time points followed by 100 time points of inter trial interval. We observed that CORNN reproduced network performed the task the same way the original network did; potentially making similar mistakes. Thus, as a first step, we confirmed that CORNN can learn to output the same outputs as the original network. Next, we looked at underlying neural activities of an example (original) generator and reproduced (CORNN learned) network in Fig. S7B. The reproduction matched not only the output of the network, but also the internal firing dynamics of individual units. Yet, as shown in Fig. S7C, the synaptic connections were not perfectly learned even in such a simple scenario with matching models, which required more trials and potentially interventional data [15]. Thus, we next tested whether increasing the number of trials could lead to better reproduction of the original network and whether CORNN can scale well to fit the increased data and network size. Fig. S8 shows that accurately reproducing the large networks required significantly more trials. Moreover, the data size requirement for the estimation process increased with the number of units in the generator network, which emphasizes the importance of an optimization algorithm that can scale well to multiple day recordings [55]. Fortunately, CORNN had a near-linear scaling of training times vs. the increasing number of trials and hence can handle large datasets (See Fig. S8). In the other direction, we observed that the training times scaled polynomially (Fig. S9), following n^β _rec CORNN learned 10% subsampled, A B Generator network dynamics C time-scale and model miss-matched Generator Model (ℳ ) network dynamics Input window Distractor Output window Input window Distractor Output window D ^I D ference Model (_ℳ ) n ^I In ^o n r ^o u _r e ue N N tu t p_t u u ⁾ _. ^p _t ⁾ d subpopulation ^/ _t ^u u _. Observe O _. O ^u _. u ^a ₍ ^/ _tu ^a ₍ pⁿ _I p Time ( ⁿ _I D ms) Time (ms) E ℳ ≠ ℳ s N n τ ≠ τ N o Rⁱ _ti ^c despite against for the of

the learning of the timed-response task. C. Despite observing only 500 out of 5000 neurons and with an average 10% mismatch between α_G and α_I , the CORNN-learned network reproduced both the neural dynamics and the output. D. When reproducing the original network without subsampling or time- scale mismatches but with generator-inference mismatches, the learned weights exhibited imperfect, but non-zero, correlation with the ground truth. E. Similar to D, we reproduced the network shown in C with subsampling and time-scale mismatches. Parameters: n_G = 5000, ǫ^RNN ∼ N (0, 10⁻²) for 100 training trials, one example network, α_G = ∆t τ = 0.1, and ∆t = 1ms. B, D: n_O = 5000, α_I = 0.1. C, E: n_O = 500, α_I ∼ N (0.1, 10⁻⁴). APPENDIX A with β ∈ [1.3, 2.2], vs. the increasing number of neurons. Bringing both observations together, we concluded that there is an O(n^β _recT ) scaling of the empirical training times. 3.5 CORNN is robust to a diverse set of non-ideal conditions To this point, we had considered ideal scenarios in which the generator network shared the same dynamical system equations as the inference network, all neurons were observed, and the time-scales of the generator and inference RNN units matched perfectly. Next, we considered cases in which these assumptions are violated (See Fig. 4A). Specifically, we trained a generator RNN with different governing equations (See Supplementary Section S3.3) to perform a timed-response task (See Fig. 4B), which put the network in a limit-cycle like dynamical attractor state initiated by an input. We tested the existence of the underlying attractor through a novel distractor input, which was not present during the training of the task (See Fig. 4B, C). If CORNN was able to learn the underlying state-space geometry, the reproduced network activities should be able to return to their attractive trajectory under novel perturbations. When we subsampled the generator networks by 10% and introduced jitter in the time scales, we observed that the inference model reproduced the neural dynamics in the observed population in a novel perturbation trial (See Fig. 4C), even though the sub-connectivity matrix was not well reproduced (See Fig.4D, E). On the one hand, this observation reinforces the findings of previous literature that the connectivity matrices learned by dRNNs should be interpreted as functional connections, accounting for the behavior of networks at the level of neural dynamics, and not synaptic connections [15, 25]. On the other hand, the ability of the mismatched network to reproduce neural activities in a novel perturbation trial provided evidence that the effects of the underlying attractor structures can be robustly reproduced with CORNN despite severe experimental non-idealities. We quantified these robustness aspects in Figs. 5A,B (subsampling), S10 (time-scale mismatch), and S11 (correlated noise). Moreover, we observed the n^β r_ec scaling for the training times with increased number of network size (Fig. 5C), inline with the results of Fig. S9. 4 Discussion To place in perspective the experimental paradigm enabled by CORNN introduced in this work, we now return to the experimental scenario in Fig. 1. Imagine an experiment in which mice perform a predefined task several times, e.g., for an half hour, with imaging of 3000-4000 neurons. This experimental scenario yields roughly ten million parameters to be trained in the dRNN. The current study showcased CORNN’s efficacy in training networks of thousands of neurons, requiring just O(10) iterations, each comparable in complexity to gradient computation and taking seconds. Consequently, training such a network from the initial imaging session would take less than a minute. Once trained, the network can enable real-time planning of experimental interventions, testing A _B y_i ^C f ^tvⁱ _t ) n c _i ⁿ ⁱ _t W o f .9 a 0.9 o 0 ₎ R _l m ⁽ 10 c ^{o a} _r 0.8 e u^r _t ⁾ _r ⁽ ⁽ 0. y ue _i ^m sn y 8 c o c ^a _r ^{n t} 0.7 g 1 ce ^a _ru 0.7 u d c ^e _{t i} ⁿ c c c _i ⁿ R c a u 0.6 ^a _r a _ts ^r _t T 0.6 e s 0.1 1^% _{10 %} ^{100 %} T noc 1 % 10 % 100 % 100 1000 5000 Subsampling percentage ^e _r Subsampling percentage out of 5000 neurons out of 5000 neurons Number of neurons Figure 5: CORNN runtimes scale polynomially with increasing number of neurons, whereas reconstruction accuracies remain robust to subsampling. We use the timed-response task as a testbed for quantifying the robustness of CORNN to unobserved influences due to subsampling in neural populations. A Reconstruction accuracy of input weights are plotted as a function of subsampling ratio (fraction of observed neurons), B reconstruction accuracy of neural activities, measured as the R² between the ground truth and predicted activations, and C the linear scaling of the training times with asymptotic slopes ∈ [1.2, 1.3] vs increasing number of neurons on a log-log plot, hinting at polynomial scaling. The different colors in the plot correspond to different number of training trials. Parameters: α = 0.1, ǫ^RNN ∼ N (0, 10⁻²) for the training trials. See Supplementary Section S3.3 for further details. Data points: mean, error bars: s.e.m. over 10 networks. APPENDIX A multiple scenarios in parallel to identify several optimal targeting strategies, a template of potential interventions, for neuron combinations driving desired animal behavior. Moreover, with CORNN’s rapid convergence and typical inter-trial intervals of a few seconds in behavioral experiments, the network and the interventional strategy can be refined between trials with incoming data streams, facilitating real-time learning. To understand the timescales for real-time inference with a dRNN for within-trial intervention experiments using calcium imaging, we can look at typical values. Acquiring each image frame takes about 30 ms [2]. While new frames continue to be acquired, motion correction and neural activity extraction can happen in ∼5 ms per frame [56]. In addition to these ongoing imaging steps, the simulation and decision-making with the CORNN-fitted dRNN needs to be performed. The dRNN starts with the current observed activity and can quickly simulate multiple future responses under different input scenarios, i.e., potential interventions. Since the simulation involves only iterative matrix multiplications and point-wise nonlinearities, it finishes very quickly—in just a few ms. For instance, running a 1000neuron RNN forward for 10 timesteps (∼ a second), for 100 different initial conditions, takes < 6 ms on a GPU, and < 25 ms on a CPU, in our hands. Once a desired intervention is identified from the simulations, the phase mask on a spatial light modulator can be updated in ∼10 ms to optically stimulate the chosen neurons [56]. Putting all the steps together, in under 100 ms one could capture brain activity, simulate future responses, decide on an intervention, and update the optical stimulation parameters. This is fast enough for real-time closed-loop applications. By streamlining network training, CORNN provides dRNNs for these types of experiments. An important final point in considering the use of CORNN in experimental settings comes from the fact that neural networks are, in general, non-identifiable [57]. That is, for any given settings of the parameters, there are other settings which give the same input-output function. This means that CORNN does not aim to infer the true underlying synaptic connectivity matrix from a neural activity dataset. Instead, the main use of CORNN is to infer an RNN model which recapitulates the dynamical trajectories in a neural population. CORNN may also capture the underlying attractor structure of a system. However, we caution that any claim having to do with attractor structures must be experimentally validated with perturbation experiments that directly test for attractors. In our work, we simulated such experimental validation (Figs. 4 and 5) and found that in the setting tested, CORNN was able to predict the dynamical effects of perturbations on the neural population. 5 Conclusion In this work, we introduced a fast and scalable convex solver (CORNN) for rapid training of data- constrained recurrent neural networks. We showed that CORNN is as accurate yet faster than existing approaches by orders of magnitude and can easily scale to the large datasets of today, training in seconds to a few minutes depending on the data, network size and the availability of a GPU. CORNN, as a base model, lays the groundwork for future developments and can be improved by integrating experimentally relevant assumptions and regularizations within the ADMM framework. When applied to simulated data from structured networks, CORNN picked up the underlying attractor structure despite non-idealities (Figs. 4, 5, S10, and S11). Inspired by this observation, further devel- opment of data-constrained RNNs can support systems neuroscience research aiming to understand inter-area interactions and can complement or augment studies of canonical or pairwise correlations toward characterizations of functional connectivity [33]. Moreover, reverse engineering the learned network may provide a deeper understanding of how neuronal populations contribute to emergent computation, cognition, and memory [12, 30]. Finally, with the advent of interventional methods such as single-cell targeted optogenetics [10], dRNNs can be trained to test causal connections computationally and supply theoretical predictions at a scale previously unachievable for single-cell targeted experimentation [8]. This work constitutes a first step towards the application of CORNN to experimental data. However, several steps remain to apply dRNNs in real-time to interventional experiments. Some example steps may include the transformation of calcium traces or spike trains into traces of firing rates normalized within [−1, 1], applying the CORNN solver developed in this work into first offline and then online experimental scenarios, estimation of neuronal time-scales from the experimental data instead of tuning them as hyperparameters, and perhaps implementing a low-rank regularization approach that opens the door to interpreting the observed dynamics in terms of latent variables [26, 52]. APPENDIX A Acknowledgements We would like to thank Liam Storan, Udith Haputhanthri, Parth Nobel, Dr. Itamar Landau, Dr. Yoshihisa Yamamoto, Dr. Surya Ganguli, Dr. Jay Mclelland, Dr. John Duchi, and Dr. Stephen Boyd for valuable feedback and insightful discussions. FD receives funding from Stanford University’s Mind, Brain, Computation and Technology program, which is supported by the Stanford Wu Tsai Neuroscience Institute. MJS gratefully acknowledges funding from the Simons Collaboration on the Global Brain and the Vannevar Bush Faculty Fellowship Program of the U.S. Department of Defense. FD expresses gratitude for the valuable mentorship he received at PHI Lab during his internship at NTT Research. References [1] Christopher Langdon, Mikhail Genkin, and Tatiana A Engel. A unifying perspective on neural manifolds and circuits for cognition. Nature Reviews Neuroscience, pages 1–15, 2023. [2] Tony Hyun Kim and Mark J Schnitzer. Fluorescence imaging of large-scale neural ensemble dynamics. Cell, 185(1):9–41, 2022. [3] Anne E Urai, Brent Doiron, Andrew M Leifer, and Anne K Churchland. Large-scale neural recordings call for new insights to link brain and behavior. Nature neuroscience, 25(1):11–19, 2022. [4] Nicholas A Steinmetz, Christof Koch, Kenneth D Harris, and Matteo Carandini. Challenges and opportunities for large-scale electrophysiology with neuropixels probes. Current opinion in neurobiology, 50:92–100, 2018. [5] Nicholas James Sofroniew, Daniel Flickinger, Jonathan King, and Karel Svoboda. A large field of view two-photon mesoscope with subcellular resolution for in vivo imaging. elife, 5:e14472, 2016. [6] Sadegh Ebrahimi, Je´roˆme Lecoq, Oleg Rumyantsev, Tugce Tasci, Yanping Zhang, Cristina Irimia, Jane Li, Surya Ganguli, and Mark J Schnitzer. Emergent reliability in sensory cortical coding and inter-area communication. Nature, 605(7911):713–721, 2022. [7] Nicholas A Steinmetz, Cagatay Aydin, Anna Lebedeva, Michael Okun, Marius Pachitariu, Marius Bauza, Maxime Beau, Jai Bhagat, Claudia Bo¨hm, Martijn Broux, et al. Neuropixels 2.0: A miniaturized high-density probe for stable, long-term brain recordings. Science, 372(6539): eabf4588, 2021. [8] Matthew G Perich and Kanaka Rajan. Rethinking brain-wide interactions through multi-region ‘network of networks’ models. Current opinion in neurobiology, 65:146–151, 2020. [9] Shenqin Yao, Peng Yuan, Ben Ouellette, Thomas Zhou, Marty Mortrud, Pooja Balaram, Soumya Chatterjee, Yun Wang, Tanya L Daigle, Bosiljka Tasic, et al. Recv recombinase system for in vivo targeted optogenomic modifications of single cells or cell populations. Nature methods, 17 (4):422–429, 2020. [10] Kayvon Daie, Karel Svoboda, and Shaul Druckmann. Targeted photostimulation uncovers circuit motifs supporting short-term memory. Nature Neuroscience, 24(2):259–265, 2021. [11] David Sussillo, Mark M Churchland, Matthew T Kaufman, and Krishna V Shenoy. A neural net- work that finds a naturalistic solution for the production of muscle activity. Nature neuroscience, 18(7):1025–1033, 2015. [12] David Sussillo and Omri Barak. Opening the black box: low-dimensional dynamics in high- dimensional recurrent neural networks. Neural computation, 25(3):626–649, 2013. [13] John J Hopfield. Neural networks and physical systems with emergent collective computational abilities. Proceedings of the national academy of sciences, 79(8):2554–2558, 1982. [14] Sung Soo Kim, Herve´ Rouault, Shaul Druckmann, and Vivek Jayaraman. Ring attractor dynamics in the drosophila central brain. Science, 356(6340):849–853, 2017. APPENDIX A [15] Abhranil Das and Ila R Fiete. Systematic errors in connectivity inferred from activity in strongly recurrent networks. Nature Neuroscience, 23(10):1286–1296, 2020. [16] Mikail Khona and Ila R Fiete. Attractor and integrator networks in the brain. Nature Reviews Neuroscience, pages 1–23, 2022. [17] Sukbin Lim and Mark S Goldman. Balanced cortical microcircuitry for maintaining information in working memory. Nature neuroscience, 16(9):1306–1314, 2013. [18] Rodrigo Laje and Dean V Buonomano. Robust timing and motor patterns by taming chaos in recurrent neural networks. Nature neuroscience, 16(7):925–933, 2013. [19] Xiao-Jing Wang. Probabilistic decision making by slow reverberation in cortical circuits. Neuron, 36(5):955–968, 2002. [20] Hidehiko K Inagaki, Susu Chen, Kayvon Daie, Arseny Finkelstein, Lorenzo Fontolan, Sandro Romani, and Karel Svoboda. Neural algorithms and circuits for motor planning. Annual Review of Neuroscience, 45:249–271, 2022. [21] Daniel LK Yamins, Ha Hong, Charles F Cadieu, Ethan A Solomon, Darren Seibert, and James J DiCarlo. Performance-optimized hierarchical models predict neural responses in higher visual cortex. Proceedings of the national academy of sciences, 111(23):8619–8624, 2014. [22] Nicolas Y Masse, Guangyu R Yang, H Francis Song, Xiao-Jing Wang, and David J Freedman. Circuit mechanisms for the maintenance and manipulation of information in working memory. Nature neuroscience, 22(7):1159–1167, 2019. [23] Guangyu Robert Yang, Madhura R Joglekar, H Francis Song, William T Newsome, and Xiao- Jing Wang. Task representations in neural networks trained to perform many cognitive tasks. Nature neuroscience, 22(2):297–306, 2019. [24] Kanaka Rajan, Christopher D Harvey, and David W Tank. Recurrent network models of sequence generation and memory. Neuron, 90(1):128–142, 2016. [25] Matthew G Perich, Charlotte Arlt, Sofia Soares, Megan E Young, Clayton P Mosher, Juri Minxha, Eugene Carter, Ueli Rutishauser, Peter H Rudebeck, Christopher D Harvey, et al. Inferring brain-wide interactions using data-constrained recurrent neural network models. bioRxiv, pages 2020–12, 2021. [26] Adrian Valente, Jonathan W Pillow, and Srdjan Ostojic. Extracting computational mechanisms from neural data using low-rank rnns. Advances in Neural Information Processing Systems, 35: 24072–24086, 2022. [27] Lea Duncker and Maneesh Sahani. Dynamics on the manifold: Identifying computational dynamical activity from neural population recordings. Current opinion in neurobiology, 70: 163–170, 2021. [28] Zach Cohen, Brian DePasquale, Mikio C Aoi, and Jonathan W Pillow. Recurrent dynamics of prefrontal cortex during context-dependent decision-making. bioRxiv, pages 2020–11, 2020. [29] Arseny Finkelstein, Lorenzo Fontolan, Michael N Economo, Nuo Li, Sandro Romani, and Karel Svoboda. Attractor dynamics gate cortical information flow during decision-making. Nature Neuroscience, 24(6):843–850, 2021. [30] Valerio Mante, David Sussillo, Krishna V Shenoy, and William T Newsome. Context-dependent computation by recurrent dynamics in prefrontal cortex. nature, 503(7474):78–84, 2013. [31] Chethan Pandarinath, Daniel J O’Shea, Jasmine Collins, Rafal Jozefowicz, Sergey D Stavisky, Jonathan C Kao, Eric M Trautmann, Matthew T Kaufman, Stephen I Ryu, Leigh R Hochberg, et al. Inferring single-trial neural population dynamics using sequential auto-encoders. Nature methods, 15(10):805–815, 2018. [32] Ruxandra Barbulescu, Gonc¸alo Mestre, Arlindo L Oliveira, and Luı´s Miguel Silveira. Learning the dynamics of realistic models of c. elegans nervous system with recurrent neural networks. Scientific Reports, 13(1):467, 2023. APPENDIX A [33] Ed Bullmore and Olaf Sporns. Complex brain networks: graph theoretical analysis of structural and functional systems. Nature reviews neuroscience, 10(3):186, 2009. [34] Karl Deisseroth. Optogenetics. Nature methods, 8(1):26–29, 2011. [35] David Sussillo and Larry F Abbott. Generating coherent patterns of activity from chaotic neural networks. Neuron, 63(4):544–557, 2009. [36] Kelly B Clancy, Aaron C Koralek, Rui M Costa, Daniel E Feldman, and Jose M Carmena. Volitional modulation of optically recorded calcium signals during neuroprosthetic learning. Nature neuroscience, 17(6):807–809, 2014. [37] Riichiro Hira, Fuki Ohkubo, Yoshito Masamizu, Masamichi Ohkura, Junichi Nakai, Takashi Okada, and Masanori Matsuzaki. Reward-timing-dependent bidirectional modulation of cortical microcircuits during optical single-neuron operant conditioning. Nature communications, 5(1): 1–12, 2014. [38] Aaron C Koralek, Xin Jin, John D Long II, Rui M Costa, and Jose M Carmena. Corticostriatal plasticity is necessary for learning intentional neuroprosthetic skills. Nature, 483(7389):331– 335, 2012. [39] Andrea Giovannucci, Johannes Friedrich, Matt Kaufman, Anne Churchland, Dmitri Chklovskii, Liam Paninski, and Eftychios A Pnevmatikakis. Onacid: Online analysis of calcium imaging data in real time. In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vish- wanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, vol- ume 30. Curran Associates, Inc., 2017. URL https://proceedings.neurips.cc/paper/ 2017/file/4edaa105d5f53590338791951e38c3ad-Paper.pdf. [40] Andrea Giovannucci, Johannes Friedrich, Pat Gunn, Jeremie Kalfon, Brandon L Brown, Sue Ann Koay, Jiannis Taxidis, Farzaneh Najafi, Jeffrey L Gauthier, Pengcheng Zhou, Baljit S Khakh, David W Tank, Dmitri B Chklovskii, and Eftychios A Pnevmatikakis. Caiman: An open source tool for scalable calcium imaging data analysis. eLife, 8:e38173, 2019. [41] Hakan Inan, Murat A Erdogdu, and Mark Schnitzer. Robust estimation of neural signals in calcium imaging. Advances in neural information processing systems, 30, 2017. [42] Aleksander Madry, Aleksandar Makelov, Ludwig Schmidt, Dimitris Tsipras, and Adrian Vladu. Towards deep learning models resistant to adversarial attacks. arXiv preprint arXiv:1706.06083, 2017. [43] Hakan Inan, Claudia Schmuckermair, Tugce Tasci, Biafra O Ahanonu, Oscar Hernandez, Je´roˆme Lecoq, Fatih Dinc¸, Mark J Wagner, Murat A Erdogdu, and Mark J Schnitzer. Fast and statistically robust cell extraction from large-scale neural calcium imaging datasets. bioRxiv, 2021. [44] Zhe Chen, Garrett J Blair, Changliang Guo, Jim Zhou, Juan-Luis Romero-Sosa, Alicia Izquierdo, Peyman Golshani, Jason Cong, Daniel Aharoni, and Hugh T Blair. A hardware system for real-time decoding of in vivo calcium imaging data. Elife, 12:e78344, 2023. [45] Steven L Brunton, Joshua L Proctor, and J Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences, 113(15):3932–3937, 2016. [46] Kadierdan Kaheman, J Nathan Kutz, and Steven L Brunton. Sindy-pi: a robust algorithm for parallel implicit sparse identification of nonlinear dynamics. Proceedings of the Royal Society A, 476(2242):20200279, 2020. [47] Dominik Schmidt, Georgia Koppe, Zahra Monfared, Max Beutelspacher, and Daniel Durstewitz. Identifying nonlinear dynamical systems with multiple time scales and long-range dependencies. arXiv preprint arXiv:1910.03471, 2019. [48] Georgia Koppe, Hazem Toutounji, Peter Kirsch, Stefanie Lis, and Daniel Durstewitz. Identifying nonlinear dynamical systems via generative recurrent neural networks with applications to fmri. PLoS computational biology, 15(8):e1007263, 2019. APPENDIX A [49] Ibrahim Ayed, Emmanuel de Be´zenac, Arthur Pajot, Julien Brajard, and Patrick Gallinari. Learning dynamical systems from partial observations. arXiv preprint arXiv:1902.11136, 2019. [50] Daniel Durstewitz, Georgia Koppe, and Max Ingo Thurm. Reconstructing computational system dynamics from neural data with recurrent neural networks. Nature Reviews Neuroscience, pages 1–18, 2023. [51] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, Jonathan Eckstein, et al. Distributed opti- mization and statistical via the direction method of Foundations and [52] Manuel Os- tojic.

[53] Jorge Nocedal and Stephen J Wright. Numerical optimization. Springer, 1999.

[54] Jian Yang, Lei Luo, Jianjun Qian, Ying Tai, Fanlong Zhang, and Yong Xu. Nuclear norm based matrix regression with applications to face recognition with occlusion and illumination changes. IEEE transactions on pattern analysis and machine intelligence, 39(1):156–171, 2016. [55] Walter G Gonzalez, Hanwen Zhang, Anna Harutyunyan, and Carlos Lois. Persistence of neuronal representations through time and damage in the hippocampus. Science, 365(6455): 821–825, 2019. [56] Zihui Zhang, loop all-

2018. [57] A [58] Peter [59] Brian force: 2018.

[60] Haim networks. Physical review letters, 61(3):259, 1988.

APPENDIX A S1 Details on CORNN derivations and implementation S1.1 Separation into subproblems and derivation of the gradient and Hessian The loss function given in Eq. (6) can be written in details as: ^ ( ^ n^T ) ∑ ∑ 1 + dˆ 1 + d ⁿ∑ θˆ) = c_t,iCE _t,i , _t, λ L( _^ ¹ _i + θˆ_j ² ^ , (S1) 1 T t=1 2 _i i= 2 j=1 2 ︸︷︷︸ L (βˆ ) where n_tot = n_rec +n_in. We note that for a given i, the functions L_i depend only on βˆ⁽ⁱ⁾ := θˆ_:,i, i.e., the ith column of the θˆ matrix. In other words, the lack of cross-talk between the columns of θˆ in the loss function means each output can be regressed independently with respect to the corresponding part of the loss function, i.e., L_i(θˆ_i). Thus, for the rest of this section, without loss of generality, we simply considered one instance of the regression, i.e., fix i, derived the learning rule, and later generalized to the full problem. The loss function of interest for the subproblem is: T ( ) ⁿ ˆ 2 L(βˆ) =

t=1 j=1 The gradient can be computed as: T [ ( )] 1 ∑ 1 + dˆ (∇ _t 1 + d_t βˆL(βˆ))_k = c ∇ CE , + λβˆ , T _t βˆ _k t=1 2 2 k T [ ] 1 ∑ c_t 1− d 1− dˆ _t 1 + d = ( _t _t)(1 + dˆ ) − x + λβˆ , (S3) T 2 _{t t,k k} t=1 1− dˆ_t 1 + dˆ_t The

At this point, it is worth emphasizing that while x_t,i are shared across all subproblems, d_t are not. In fact, for any neuron i, a different set of d_t,i would contribute to the Hessian of the corresponding subproblem. Thus, in its current form, using the exact Hessian would require repeating matrix inversions for each of the N_rec neurons. We performed benchmarking with the exact Newton method, See Supplementary Section S2, though the repeated matrix inversion was no longer feasible for large networks (Fig. 3). To mitigate this, we took a fixed-point solver approach in CORNN, described below. First, we assumed that one can initialize the full problem with an initial point, call it d¯_t = tanh(βˆT f_px_t), such that d¯_t ≈ d_t. We discussed in Supplementary Section S1.2 on how to estimate βˆ_fp. Then, we replaced the Hessian with a local approximation∇^{2 2} βˆL(βˆ) ≈ ∇ βˆL(βˆ_fp) within the vicinity of d¯_t ≈ d_t, leading to the f_p. (S5) From this The local depends only and not

APPENDIX A the predicted or true outputs (dˆ_t or d_t). Moreover, it needs to be computed once, and at the beginning of the optimization only. Using the approximate Hessian and gradient, the parameter update rule simply follows: ∆_fp = −H−1 (βˆ). (S6) Even all (θˆ) to solve i.e., At this the a

semi-definiteness. Thus, even if the error between d¯_t and d_t is large for an instance, e.g., fixed-point is incorrectly chosen, it can only hurt the speed of the optimization. Bringing all together, fixed-point updates can be written as: θ^(k+1) := θ^(k) − γ_kH−1 f_p ∇_θL(θ^(k)), (S7) where γ_k is the user-controlled learning parameter chosen either as 1 (to bypass the line search if fixed-point was close enough) or through back-tracking line search. For this work, we simply picked γ_k = 1 without any line search as we did not observe any convergence issues. We will simplify Eq. (S7) to its final form when discussing the ADMM solver in Supplementary Section S1.4. S1.2 Fixed point initialization To derive the fixed point update in Eq. (11), we performed a Taylor approximation of the gradient and carried out a single update near the vicinity of the putative fixed-point: βˆ^(t+1) : = βˆ^(t) −H_f ⁻ p¹ ∇_βˆL(βˆ^(t)), ( )₋₁ [( ) ] (^t) 1 1 1 ≈ βˆ − X^TX + λI X^TX + λI βˆ^(t) − X^TZ , T T T ( ) (S8) t^{) (t)} 1 ₋₁ (^T 1 = βˆ − βˆ + X X + λI = where a single-step

problem on the currents Z. Thus, βˆ_ls is expected to be near to both the solution β^∗ and the regime described by the fixed-point Hessian H_fp; thus a suitable candidate to be the fixed-point itself. As an added bonus, this choice did has already been pre- S1.3 Interpretation of the

To interpret CORNN’s loss function, we consider the low error limit of the cross entropy CE(pˆ, p) around pˆ ≈ p: CE(pˆ, p) = −p

d ^∣ = −H(p)− [p log(pˆ) + (1− p) log(1− pˆ)]∣_∣ (pˆ− p) pˆ=p

APPENDIX A where −H(p) is a constant independent of the optimization variables, and thus can be ignored. Then, the low-error limit of the CORNN loss is: 1∑ ^T ( ) 1 + dˆ_t 1 + d_t 1∑ ^T (d_t − dˆ_t)2 c_tCE , ≈ 2, T 2 2 T (1− d² t=1 _t )

)[zˆ − f⁻¹(d) ² ≈ _t]) , (S9) T

t=1 2(1 1∑ ^T [zˆ − f⁻¹(d_t)]² = . T t=1 2 Here, in the second row, we

+ (dˆ^′ _t)|_{dˆ =d}(zˆ−f⁻¹(d_t)). In this low error limit, minimizing the scalable weighted loss approximately accounts to minimizing the L₂ loss on the currents. The main difference is, since the original weighted loss is on the firing rates and not currents, the weighted loss can account for the conversion noise, or non-linearity mismatches, that would induce bias in the explicit conversion z = f⁻¹(d).

So far, we discussed the fixed-point solver for the unconstrained CORNN loss with L2 regularization. In this section, we scenario, it is unexpected that all

one might be interested in introducing a = 0. To S1.4.1 In this the ADMM

1∑ ^T 1 + dˆ_t 1 + d ⁿ∑ (βˆ 2 minimize L _t

to ∀k ∈ B(k), βˆ_k = 0. We transformed

1∑ ^T 1 + dˆ 1 ⁿ∑ 2 c _t + d_t (βˆ) minimize _j _tCE , + λ , T t=1 2 2 j=1 2 (S11) subject to ∀k ∈ B(k), χˆ_k = 0, χˆ = βˆ. Defining the indicator function 1(x) as 1

The ADMM loss function became: 1∑ ^T ( ) 1 + dˆ_t 1 + d_t ⁿ∑ (βˆ)2 ∑ minimize L (βˆ, χˆ) _j _ADMM = cCE , + λ + 1(χˆ = 0) , T _t 2 2 _k _︸ t=1 ︷_︷ j=1 ² ︸

k∈B(k) ︸︷︷︸ (χˆ) subject to χˆ =

(S13) APPENDIX A ADMM Lagrangian has two separable loss functions of independent variables, connected through a linear equality constraint, which can be augmented to obtain: ρ(βˆ, χˆ, κ) = L_ADMM(βˆ, χˆ) + κ^T ρ L (βˆ − χ) + ||βˆ − χˆ||², (S14) 2 ₂ where κ is the dual variable. The ADMM steps (after re-defining κ → (1/ρ)κ) became (See Eqs. (3.5-3.7) in [51]): ( ) βˆ^k+1 := argmin L_CORNN(βˆ) + (ρ/2)||βˆ − χˆ^k + κ^k||² 2 , (S15a) βˆ (

) χˆ^k+1 := argmin − χˆ + κ^k||² 2 , (S15b) χˆ κ^k+1 := κ^k + βˆ^k+1 − χˆ^k+1. (S15c) S1.4.2 Solving first ADMM iteration via fixed-point updates Similar to the unconstrained case, we started by finding the gradient and the fixed-point Hessian: [ ∇ L (βˆ) + (ρ/2)||βˆ − X^TE ^{k k} βˆ _CORNN + (λ+ ρ)βˆ + ρ(κ − χˆ ), (S16a) [ ∇² βˆ L_CORNN(βˆ) + (ρ/2)||βˆ −

^TX + (λ+ ρ)I, (S16b) where we recall the prediction error matrix: d E _t,i − dˆ_t, _t,i = _i . (S17) 1− d² t_,i Then, a single fixed-point update step became: [ ]₋₁ [ ] 1 ∆_fp = − X^T 1 X + (λ+ ρ)I − X^TE + (λ+ ρ)βˆ + ρ(κ^k − χˆ^k) , T T [ _] ^[ ] T ₋₁ (S18) = X X + T (λ+ ρ)I X^TE − T (λ+ ρ)βˆ + Tρ(χˆ^k − κ^k) , = A⁺X^TE − (λ˜ + ρ˜)A⁺βˆ + ρ˜A⁺(χˆ^k − κ^k), with λ˜ = Tλ and ρ˜ = Tρ. Here, we defined the short-hand notation for the inverse matrix: [ ]₋₁ A⁺ = X^TX + (λ˜+ ρ˜)I . (S19) Consequently, the update rule for γ_k = 1 became: βˆ^k+1 = βˆ^k +A⁺X^TE^k − (λ˜+ ρ˜)A⁺βˆ^k + ρ˜A⁺(χˆ^k − κ^k), + (S20) = A X^TXβˆ^k +A⁺X^TE^k + ρ˜A⁺(χˆ^k − κ^k), where we highlighted all quantities that can be pre-computed with blue. For ρ = 0, this reproduces the update rule for the unconstrained problem in Eq. (10). As before, we promoted the vector βˆ to a matrix θˆ, hence vectorized across the full problem. S1.4.3 Solving second ADMM iteration via projection For this step, the goal is to minimize ∑ ⁿ minimize (ρ/2)(βˆ_j ^k+1 − χˆ_j + κk j )², j=1 (S21) subject to ∀k ∈ B(k), χˆ_k = 0. Fortunately, this problem was perfectly separable across the vector entries of χˆ_j . For cases where χ_k = 0 was constrained, this was the only feasible point and thus the answer to that particular problem. For when j ∈ / B(j), we obtained the minimum loss value by simply picking χˆ_j = βˆk+1 k j + κ_j . In short notation, the update rule for this subproblem became: ( ) ( ) χˆ^k+1 := argmin L_constraint(χˆ) + (ρ/2)||βˆ^k+1 − χˆ + κ^k||² 2 = Π_B βˆ^k+1 + κ^k , (S22) χˆ where we defined Π_B as the projection operator to the constraint satisfying subspace (i.e., ∀k ∈ B(k), χˆ_k = 0). Similar

APPENDIX A S1.4.4 Variable initialization and choice of ρ A priori, it was not clear how ρ should be chosen. However, looking at the Hessian inverse A⁺, we observed that ρ ≥ λ was an appropriate choice for the step size ρ to have a reasonable impact on the optimization problem. Empirically, we observed that ρ ≫ λ lead to faster convergence. We initialized the primal and dual variables with χ⁰ = βˆ⁰ = βˆ_fp, i.e. the pre-computed least-squares solution, βˆ_fp = A⁺X^TZ , without additional computational complexity. Moreover, given that ADMM ensured convergence χˆ = θˆ, the dual variable κ was bound to converge to non-zero values only at the diagonal such that the 2nd ADMM step would converge. To computational cost, we defined the initial conditions for all three

βˆ⁰ = A⁺X^TZ, (S23a) χˆ⁰ = A⁺X^TZ, (S23b) κ⁰ = 0. (S23c) S1.5 Automated outlier detection and the pseudo-code for the full solver Upon close inspection of the that the extremely large error values stemming from d_t ≈ 1 would noise and carry little-to-no learning signal. Thus, we implemented a simple errors larger than some value such that |E_t,i| > v =⇒ E_t,i = 0. In the CPU the gradient to account for the missing data points. Since these data points mostly < 1% of the full dataset, we did not perform the scaling for the GPU unnecessary overhead. It is worth emphasizing that this takes away the convexity of the problem since zeroing out part of the prediction

convex global loss. However, the outliers in

would automatically detect outliers. For the purpose of this work, we observed that the version of outlier detection we employed here performed well on empirical benchmarks. For all solvers, we also defined a working precision δ (= 10⁻⁶) and whenever |d_t| ≥ |1 − δ|, we projected d_t back to ±(1− δ) to prevent numerical overflow errors. For experiments in Figs. S1 and S6 with BPTT, we projected the firing rates |r_t| ≥ |1 − δ| to the boundary r_t = ±(1 − δ) for the same reasons. Empirically, this stabilized all solvers and was performed to ensure a fair comparison. Bringing all together, we state the full CORNN solver in Algorithm 1. We note that for the figures in this work, we simply used a maximum number of iterations (30-100 depending on experiment) instead of monitoring the convergence of CORNN; as CORNN runtimes remained sub-minute and a few more iterations than ideal were perfectly tolerable. S2 Implementation details for other solvers S2.1 Back-propagation through time We performed experiments, the (2), but we replace xˆ_t,i on θˆ implicitly, so

The probability of teacher forcing is the probability with which the replacement is performed at every single time point during the forward propagation of the network. Specifically, we use the following forward-propagation equations: { (1− α)r_t,i + αf(z_t,i), u < p, rˆ_t+1,i = (S24) (1− α)rˆ_t,i + αf(zˆ_t,i), u ≥ p, APPENDIX A Algorithm 1: Convex and scalable solver for CORNN via ADMM framework CORNN.fit (X,D, λ, nIter = 30, ρ_in = 100); % X contains the aggregated input, e.g. firing rates and input units. It is of size [T ×N_tot]; % D contains the targets. It is of size [T ×N_rec]; % To start, rescale the regularization parameter λ := Tλ ρ := ρ_inλ %ρ is picked to be comparable to entries of X^TX and rescaled λ; % Pre-compute the inverses once during initialization; A⁺ := (X^TX + (ρ+ λ)I)⁻¹; X⁺ := A⁺X^T ; X⁻ := A⁺X^TX = X⁺X; Z := f⁻¹(D) %Compute the currents from the targets ; θ_fp := X⁺Z%Initialize the fixed-point to the ls solution; % Initialize the primary and dual variables; θ := θ_fp %Initialize to the least-squares solution; χ := θ_fp %Initialize to the least-squares solution; u := 0 ∗ θ_fp %Initialize the dual variable; for i = 1 : nIter; Predict the targets Dˆ = f(Xθ^k); Compute the prediction error E^k = (D − Dˆ)⊘ (1−D ⊙D) and perform outlier detection % ⊘/⊙ stand for element-wise division/product ; U^{pdate first primal variable θk+1} _{:= X} ⁻ _θ ^k ₊ ( _X ⁺ _E ^k ₊) _ρA ⁺ _(χ ^k _{− κ} ^k _); Update second primal variable χ^k+1 := Π_B θ^k+1 + u^k ; Update the dual variable u^k+1 := κ^k + θ^k+1 − χ^k+1; Probe convergence by checking ||χ^k+1 − θ^k+1||, continue if not converged; end return χ where u is a uniform random number generated between 0 and 1. It is worth emphasizing that the teacher forcing in this context is distinct from [59], where the teacher signal is acquired not from the real-data but from a second network, though the overall idea is similar. S2.2 Current and firing rate based FORCE We followed the approach by [25] for training the recurrent weights to reproduce target neural activities via FORCE. For both network types, we defined the prediction errors as: e[t] = (rˆ[t]− r[t])/α, or e[t] = xˆ[t])− x[t]. (S25) Here, rˆ[t] is the predicted firing rate at time point t, vice versa for xˆ[t]. We no longer performed any teacher forcing, thus r[t− 1] and x[t− 1] needed to predict rˆ[t] came from the previous prediction step of the network, and were not the ground truth values. The recurrent weights and the corresponding inverse covariance matrix P [t] were updated by the following rule: P [t− 1]rˆ[t− 1]rˆ^T [t− 1]P [t− 1] P [t] = P [t− 1]− ^T , (S26a) 1 + rˆ [t− 1]P [t− 1]rˆ[t− 1] W_rec[t] = W_rec[t− 1]− e[t](P [t]rˆ[t− 1])^T . (S26b) Here, e[t] was an N_rec dimensional vector, rather than a scalar as was the case for the original work [35]. We initialized P [t] = λ⁻¹I , where λ is the regularization parameter. Short-comings of FORCE for dRNN training Despite its wide use in the literature for dRNN training [24, 25, 27–29], there are several short- comings of this of RNNs using FORCE that render it incompatible with the realities of experiments:

APPENDIX A 1. The training procedure requires a particular range of λ to converge (See Fig. S4). Even when FORCE does converge, this is for a very limited range of λ and epochs of training. 2. FORCE is initially designed for theoretical investigations and requires the training dataset to be continuous. Even subtle jumps in the data results in stability issues with FORCE training, as the premise of FORCE approach is to have small error from the start to the end [35]. Point 1 necessitates low-level hyperparameter tuning. The only work around to point 2 is to re- initialize the inverse covariance matrix P whenever there is a jump in data points. This practically means that if we want to perform multiple training steps on the same data, each epoch simply corresponds to picking a better initialization for the next epoch, whereas all the learned correlation structure, which is stored in the inverse covariance matrix P , needs to be forgotten. The latter point showcases that FORCE is fundamentally ill-suited to many aspects of experimental data, for example, when multiple imaging sessions need to be combined. S2.3 Newton’s method and gradient descent with single-step prediction error We implemented Newton’s method by using the exact Hessian given in Eq. (S5) for two different scenarios: i) = 1 = [1 − d² _t ]⁻¹

detection not only on the gradient, but also the Hessian level; since the division by [1 − d² _t ] could have destabilized the training. Given that each subproblem had a unique Hessian for the Newton descent, we solved the separate subproblems in parallel using multiprocessing tool in numpy. To implement the gradient descent, we used a single-layer Pytorch model with a linear + tanh structure. We implemented both a CPU and GPU version, and used the ADAM and SGD optimizers. Gradient descent was implemented for two types of loss functions: i) loss function, and ii) cross-entropy with c_t = 1. S3 Benchmarking details for generator RNNs S3.1 Randomly initialized chaotic RNNs When benchmarking dRNNs on randomly initialized chaotic RNNs, w (e picke )d the recurrent weights of the generator randomly from a normal distribution W_i ^r j^ec ∼ N 0, ^g n with zero mean and standard deviation of ^{√ g} n . For g > 1 and sufficiently large n_rec, these networks can sustain their activities without any external stimulus and exhibit chaotic behavior [60]. For these chaotic generators, we picked g = 3 throughout this work to ensure the networks were deep into the sustainable spontaneous activity regime. We chose the chaotic regime specifically, since previous work established that teacher-forcing when training an RNN to output a function can have potentially harmful consequences to the stability of the network training [35]. We aimed to test whether this generalized to the reproduction task we considered in our work. We note that in figures with randomly initialized chaotic RNNS, where we compared CORNN to other optimizers, we omitted the zero self-excitation constraint, i.e., W_i ^r i^ec = 0, to remain conservative in our speed comparisons. This is because 2nd order methods (FORCE and Newton descent) cannot incorporate the zero self-excitation constraint without additional optimization steps since projected 2nd order steps are not necessarily descent directions, whereas projected gradient descent is notoriously slower compared to the unconstrained version. In contrast, CORNN was approximately equally fast with or without the ADMM solver; since both ADMM and HAPE updates converged in more or less the same number of iterations. S3.2 Training of RNNs for the 3-bit flip flop task To benchmark CORNN on diverse conditions, on top of the randomly connected RNNs discussed above, we established ground truth data by training RNNs on a 3-bit flip flop task. The architecture of the generator RNN was as follows: r(t) = (1− α)r(t− 1) + α tanh [W_inu(t) +W_recr(t− 1)], (S27a) o(t) = W_outr(t), (S27b) APPENDIX A where r(t) is the hidden state of the RNN at time t, α is the time-constant associated with the dynamics, W_in is the weight matrix governing input into the RNN, u is the input into the RNN, W_rec are the recurrent weights, o is the output of the RNN and W_out are the weights which take the hidden state to the output. Since CORNN does not learn W_out, to compute the outputs for the CORNN learned networks in Fig. S8, we used the same projection matrix W_out as the original network. See the accompanying code for further details on network training and benchmarking experiments. S3.3 Training of RNNs for the timed-response task To test robustness of CORNN to mismatches between generator and inference network dynamics, we established ground truth data by training leaky current RNNs on a timed-response task with the following architecture: x(t) = (1− α)x(t− 1) + αW_inu(t) + αW_recr(t− 1) + αǫ_RNN, (S28a) r(t) = tanh(x(t)), (S28b) o(t) = W_outr(t). (S28c) Here, unlike the CORNN and 3-bit flip flop equations (See Eq. (S27)), the time-derivative was at the level of currents (x), not firing rates (r). Upon receiving an input during t ∈ [0, 100]ms, we trained 10 networks to output a Gaussian pulse centered at t = 500ms with standard deviation 20ms; otherwise zero everywhere. During test trials in Figs. 5, S10, and S11, we introduced a novel distractor input ten times stronger than the cue at t ∈ [400, 410]ms. To allow fair comparison between networks trained on different numbers of neurons, we computed the R² of the neural activities from the first 100 neurons across three 1s trials and performed a uniform average. For the illustrative Fig. 4, to allow visualization of finer details in time-activities after distraction, we instead plotted a scenario with an earlier and weaker (half strength than the cue) distractor at t ∈ [200, 210]ms. To obtain the correlated noise in Fig. S11, we first sampled random noise from a normal distribution, then convolved the i.i.d. noise with a 2D Gaussian kernel across time points (s.d. of 3ms) and neurons (s.d. of 5 neurons), and re-scaled by the targeted standard deviation. See accompanying code for further details. S4 Benchmarking CORNN’s speed and scalability To obtain the CORNN’s solver in Algorithm 1, we took five distinct steps: 1. Single-step prediction error paradigm, where we teacher-forced rˆ_t,i, the RNN unit activities at time t, to the observed neural activities, r_t,i, to predict r_t+1,i, with a single-step back- propagation. 2. Convexification of the loss function, where we replaced the L₂ loss function with a cross- entropy loss. 3. Hessian alignment to neural activity correlations, where we weighted individual samples to align subproblem Hessians to a pre-computable quantity, i.e., the correlation matrix of observed neural activities. 4. Initialization to a fixed-point, where we used the approximate least square solution to initialize the networks. 5. Using ADMM approach to enforce constraints without increasing algorithmic runtimes. In this work, we performed several benchmarking experiments to delineate the contribution of each step we take towards the final convex solver to the scalability and speed of CORNN. We proved, with targeted experiments, the need for each of these steps, as we discuss below. S4.1 Single-step prediction error paradigm As discussed in Section 2.2, the first step we took was to restrict the family of interest for the inference models to the models that can perform single-step predictions well. This simplifying assumption mapped the recurrent estimation problem to a single layer feed-forward network with a sigmoid/tanh non-linearity. Here, the inputs of the network are u_t,i and r_t,i , whereas the outputs are r_t+1,i. Since APPENDIX A this step eliminated the need for the back-propagation through time, both the training times and memory requirements decreased significantly (Compare Figs. S1 and S2). However, a priori, it was not clear whether the restrictive assumption of single-step prediction error minimization would lead to decreased performance. To test this, we performed the experiments in Fig. S1 and observed that teacher forcing increases the efficiency of the training, and was not deleterious. A potential explanation may lie in the observation that, unlike the traditional scenario where a sparse error signal needs to be propagated back in time to perform credit assignment, network reproduction training had access to rich variety of error signals at each time step and for each neuron. S4.2 Convexification of the loss function As discussed in Section 2.3, the convexification of the loss function plays an important role in speeding up the learning process. However, by itself alone, convexification is not the reason of the speed improvement. For example, in Fig. 3, gradient descent on L₂ loss converges faster than the gradient descent on the logistic loss. On the other hand, for smaller networks with 200 neurons (Figs. S2, S3, S5), convex loss function merits the use of Newton descent and leads to faster convergence than gradient descent. Yet, as shown in Figs. S2 and S5, CORNN is at least an order of magnitude faster even for these small networks. Hence, though convexification of the loss function is crucial, it cannot explain the speed or the scalability of CORNN by itself. S4.3 Hessian alignment The main contribution of our work is to pick c_t,i, which approximates the problem as a local least- squares as discussed in Supplementary Section S1.3, and aligns the Hessians of each subproblem to a pre-computable quantity, i.e., the correlation matrix of neural activities, as discussed in Section S1. This alignment allows the use of HAPE updates that have complexity of gradient descent but converge in O(10) iterations thanks to the approximate Hessian. Figs. S2 and S5 show that Newton’s descent with exact Hessian computations on the weighted CORNN loss is an order of magnitude slower than CORNN. Hence, the Hessian alignment was a cruical step to speed up the learning, but most importantly allow scaling to large networks due to its gradient-level complexity. S4.4 Fixed-point initialization The fixed-point initialization relies on an approximate solution when the final error is low (See Supplementary Sections S1.2 and S1.3). This initialization is a good first guess, since, if there was no conversion noise at all, one could have simply performed least-squares on the currents estimated via z = arctanh(r). This is akin to regressing to the inverse of the output before the non-linearity in a single layer network. Though such an approach can provide a good first guess, it is not robust to conversion noise. This was shown in Figs. 3 and S3, where non-Gaussian noise components lead to low accuracy for the initialization. Finally, the least-squares solution does not satisfy any additional constraints on the connectivity matrix, thus is not suitable as a valid solution itself. Yet, since the least-squares solution is computationally negligible, we used it as an initialization for CORNN. S4.5 Alternating direction method of multipliers We chose to use ADMM approach to enforce the no self-excitation constraints and present a ver- satile base model that can incorporate additional design choices such as L1 and/or nuclear norm regularization. In our comparison experiments, we omitted any constraints since other solvers were considerably slower, or outright incompatible. However, in our speed benchmarking results in Figs. 5, S8, and S9, we used the full CORNN solver with the ADMM described in Algorithm 1. The choice of ADMM for enforcing constraints is rooted in our observation that both the ADMM steps and the original HAPE updates take O(10) iterations to converge. Moreover, other parts of the ADMM updates are negligible compared to the primal problem, where the HAPE updates take place. In other words, the fact that HAPE updates and ADMM converge in similar steps creates a mutual synergy that increases accuracy and promotes versatility in our work without hurting convergence times. APPENDIX A S5 Reproducibility We provided the code to reproduce each figure in our Github. All experiments, except for those in and were run on a standard with Geforce RTX 3090 Ti GPU and Intel two due See the

APPENDIX A A B Runtime per experiment is 47 ± 1 seconds n Runtime per experiment is 47 ± 1 seconds no W i_t W_f 0.8 oⁱ _{t f} o 0.8 cu ^o ^r _t ⁾ _r cu ⁾ s ⁽ Cross entropy loss ^r _{t r} s ⁽ Cross entropy loss n y o c L2 l n yc L2 loss c ^a e _r oss u 0.2 oc ^a _r 0.2 Rc e u c a t Rc r o_i ⁿ r c t r n _h ^er _h ^e _h ^er _h ^er _h ^e a _i ^{n e}r _h ^{e e}r ^er _h ^er d_{p i} ^o _a ^c _a ^c _a ^c _a ^c _a ^c _dp ^o _i ^on _a ^c _{h a} ^c _a ^c _{h a} ^c _{h a} ^c Fⁱ _x ^e i_t ^ia _li^z _a ^t t N^{o e} _i _f ⁿg t o_r ^c ₂ ⁵%^e _i _f ⁿg t^e ^o _r ^c _i ₅ ⁰% _r ^cng t^e _i ₇ ⁵% _r ^cng _t ₀ ^e ^% _r ^c _i ⁿg _F ⁱ _x ^e i_t ^ia _l ^iz _a ^t t N^{o e} _i _f ⁿg ^te _i ⁿg t^e ^o _r ^c ₂ ⁵% _f ^o _r ^c _ig t^e _ig t^e _ig 5⁰% _r ^cn ⁵% _r ^cn ₀ ^% _r ^cn iⁿ _f ^o _f ^o ₁ ⁰ _f ^o in _f ^o _{7 f} ^o ₁ ⁰ _f ^o of loss

since cross-entropy with teacher forcing emerged as the superior method compared to BPTT. This superiority may be attributed to the rich information content available during network reproduction training. Parameters: A, B: α = 0.9, n_rec = 100, T = 100, n_epoch = 3000. ∀t = 1, ... , T, u(t) = 0.1 with W ⁱⁿ ∼ N (0, g²) and g = 3. A: ǫ^conv ∼ N (0, 10⁻⁸), ǫ^input ∼ N (0, 10⁻⁴). B: ǫ^conv ∼ N (0, 10⁻¹⁰), ǫ^input ∼ N (0, 10⁻²). Box plots are obtained from 10 random generator networks shared across optimization algorithms. Orange lines denote median values; boxes span the 25th to 75th percentiles; outliers are not shown. CPU implementations A B CORNN noⁱ _t W_f 1 Newton descent - weighted no c ^o cross entropy ⁱ _t W_f o 0.97 u^r _t ⁾ _r Newton descent - cross entropy cu ⁾ s ⁽ _r n y Gradient descent (ADAM) - ^r _ts ⁽y o c a 0.9 cross e no c ce _r ntropy u c ^a _r Rcc Original FORCE applied e u Rc 0.93 a to c 0.1 1 10 currents (z) a 0.1 1 10 Training time (s) Modified FORCE applied to firing rates (r) Training time (s) Figure S2: CORNN achieves a better combination of speed and accuracy compared to al- ternatives, dramatically accelerating convergence. The plots illustrate the relationship between accuracy (Pearson’s correlation coefficient between ground truth and inferred weights) and training time, measured in seconds on a log scale, with low (A) and high (B) input to conversion noise ratio. The blue line represents the performance of our newly developed method, which combines fixed-point initialization and CORNN, and achieves high accuracy while significantly accelerating convergence A: all is

noⁱ W t_f Newton descent - weighted c ^o cross entropy noise injections. Same u^r _t ⁾ _r Newton descent - cross entropy as in Figure S2, but with s ⁽ n y Gradient descent (ADAM) - oc parameters: α = 0.1, c ^a e_r cross entropy u n_rec = 200, T = 1500, Rcc 0 Original FORCE applied a to currents (z) No input. ǫ^conv ∼ 0.01 1 100 Poisson(10⁻¹), ǫ^input ∼ Training time (s) Modified FORCE applied to firing rates (r) Poisson(10⁻²). APPENDIX A show

not. Moreover, cases, cannot the fixed-point initialization we developed in this work. A with fixed point initialization and B with random initialization. Parameters are the same as in Fig. S2A. In all plots, data points indicate mean values, and error bars denote SEM over 10 runs of the simulation.

Figure S5: CORNN can be drammatically accelerated with on a standard GPU. The same as in Fig. S2, but with GPU acceleration for CORNN and gradient descent algorithms. APPENDIX A A Runtime per experiment is 56 ± 2 seconds B Runtime per experiment is 56 ± 2 seconds noⁱ _t W_f 0.8 noⁱ W t_f c ^o c ^o 0.6 u^r _t ⁾ _r u^r _t ⁾ _r but

∼ , ∼ . are ran on cover from the lower quartile to the upper quartile, with an orange line at the median. For C, in all plots, data points indicate median values, and error bars denote SEM over 100 runs of the simulation. Input and output for the fÄirst trial Trace of neural activations Weights Original CORNN learned

de n_ra l e a i ^k ⁿ _r s L o n gⁱ _r ^w _t ^o _r O e u n e N Recurrent Original dRNN

a 3-bit flip-flop task [12], where the network is given a short input pulse at a random time. If the input pulse is +1, the output of the network must transition to +1 or stay at +1. If the input pulse is -1, the network’s output must transition to -1 or stay at -1, ignoring input pulses that do not match. Through dRNN training, CORNN is able to accurately reproduce, left, the inputs and outputs of the task-trained RNN, middle, spatiotemporal pattern of the neural population dynamics, and, right, the inputs and (partially) recurrent weights of the synthetic RNN. Parameters: α = 0.9, n_rec = 1000, T = 200 trials (100 data points per trial), ǫ^conv ∼ N (0, 10⁻⁶), ǫ^input ∼ N (0, 10⁻⁴). APPENDIX A CPU implementation GPU implementation (64 bit) GPU implementation (32bit) 1 1 1 noⁱ _t W_f cu ^o ^r _t ⁾ _r s ⁽ n yc Number of oc neurons e ^a _ru 100 R cc a 0.2 1000 3000 0.2 0.2 fo tu 1 1

i^{m t} 1 g iⁿ 0.01 iⁿ 0.01 a_r T 0.01 1₀₀ ^One ₁₀₀₀ One 10 10 100 1000 ₁₀₀ One day month ^100One ₀ day month Number of trials Number of trials Number of trials Figure S8: CORNN runtimes scale linearly with increasing trial data, which can facilitate learning from larger networks. Top, reconstruction accuracy of network weights, measured as the correlation between recurrent weights, are plotted as a function of dataset size (number of trials), middle, reconstruction accuracy of network outputs, measured as the R² between the output units, and, bottom, the linear scaling of algorithmic wallclock runtimes with near-unity slope vs increasing number of trials on a log-log plot. The different colors in the plot correspond to different number of neurons in the generator network. Missing data points correspond to the limits of the GPU memory, which can be mitigated by computing prediction errors in temporal chunks in real-data applications, though not implemented for this case study. Parameters: α = 0.5, ǫ^conv ∼ N (0, 10⁻⁶), ǫ^input ∼ N (0, 10⁻⁴). In all plots, data points indicate mean values, and error bars denote SEM over 20 runs of the simulation.

Figure S9: CORNN scales polynomially to large network sizes. Using the networks with 3000 neurons trained for the 3-bit flip flop task in Figure S7, we subsampled neural populations and collected CORNN runtimes for CPU and GPU implementations. Parameters: α = 0.5, ǫ^conv ∼ N (0, 10⁻⁶), ǫ^input ∼ N (0, 10⁻⁴). A with 200 trials (20,000 data points) and B with a single trial (100 data points). APPENDIX A A y_i B f ^tvⁱ _t 1 C 0. 0.4 o c 8 ) a R_l 0 n (^a _r oⁱ _t W_f no W c ^{o i} _{t f} y u c e -1 u ⁾ ^a _r ^{n r} _{t r} c ( u ^o ^r _t ⁾ _r ⁽ 0.2 u d s y s c -2 n c n yc c ^e _t a c oc ^a _r oc ^a _r _t u Į^a^1^^^^^>V^G^@) e uc 0 e u s ^r _t -3 Rc .2 Rc e s c T noc 0.01 0.1 1 a 0.01 0.1 ¹ a 0.01 0.1 1 e _r 8QFHUWDLQW\^^V^G^^^LQ^WLPH^VFDOHV^^Į^ 8QFHUWDLQW\^^V^G^^^LQ^WLPH^VFDOHV^^Į^ 8QFHUWDLQW\^^V^G^^^LQ^WLPH^VFDOHV^^Į^ Figure S10: CORNN is robust to variations in the time-scales compared to the ground truth. Same as in Figure 5, but with 500 observed neurons and varying levels of uncertainties in time-scales.

1 y_i ^t f v o ⁱ _t 0 ) ca Unsuccessful reproduction level R _l a_r Type of

noise e 30 T no Independent noise Figure S11: Correlated noise presents a chal- c 50 Correlated noise e _r 100 lenge for CORNN, which can be mitigated by increasing the data size. Same as in Figure 5, 10 10 10 but with the fully observed network and varying Standard deviation of the injected noise levels of injected random and correlated noise (for details, see Section S3.3). APPENDIX A CODE FOR EXPERIMENT 1 #%% optimize_vaiables.py w1 = solve_newton_descent(r[:-1,:],r[1:,:],u_in = #!/usr/bin/env python3 None,alph = 0.1, # -*- coding: utf-8 -*- l2 = 1e-5,verbose = 1,threshold =0.2, """ initialize_fp = 0,num_iters = 10, Created on Wed Dec 2117:52:092022 solver_type = 'weighted') @author: dinc w_rec_nt = w1[:opts['n_rec'],:].T.ßatten() """ #%% import numpy as np from scipy.stats import pearsonr w4 = Þt_FORCE_gpu(r,None,alph = 0.1, from rnn_class import RNN lam = 100,g_in = 2,verbose = 2, import matplotlib.pyplot as plt initialize_fp = 1,num_iters = 100, from numpy.linalg import inv gnd = gnd,solver_type = 'Þring_rates') from utils import solve_corrn_gpu from utils import solve_corrn w_rec_fc = w4[:opts['n_rec'],:].T.ßatten() from utils import solve_newton_descent from utils import solve_gradient_descent #%% from utils import solve_pytorch_gpu plt.scatter(gnd,w_rec_pt) from utils import solve_pytorch plt.scatter(gnd,w_rec_co) from utils import Þt_FORCE from utils import Þt_FORCE_gpu print('NT:', pearsonr(w_rec_nt,gnd)[0]) opts = {}; print('GD:', pearsonr(w_rec_pt,gnd)[0]) opts['g'] = 3; print('CoRNN:', pearsonr(w_rec_co,gnd)[0]) opts['n_rec'] = 200 opts['n_in'] = 0 opts['sigma_input'] = 1e-2 opts['sigma_conversion'] = 1e-4 rnn_class.py T_data = 3000; #!/usr/bin/env python3 m1 = RNN(opts) # -*- coding: utf-8 -*- r = m1.get_time_evolution(T = T_data) """ gnd = m1.rnn['w_rec'].ßatten() Created on Fri Jun 1014:48:402022 #%% @author: dinc w = solve_corrn(r[:-1,:],r[1:,:],u_in = None,l2 = """ 1*1e-5, import numpy as np verbose = 2,gnd = gnd, threshold = 1, import multiprocessing num_iters = 30,check_convergence = 0, from joblib import Parallel, delayed solver_type = 'weighted',alph = 0.1) from scipy.stats import pearsonr w_rec_co = w[:opts['n_rec'],:].T.ßatten() from numpy.linalg import solve from scipy.linalg import solve as solve2 #%% w3 = solve_pytorch_gpu(r[:-1,:],r[1:,:],u_in = import time as time_now None,alph = 0.1, l2 = 1e-8,verbose = 2,gnd = gnd, class RNN(): initialize_fp = 0,num_iters = 5000, # DeÞne some global variables that are mostly learning_rate = 0.01,solver_type = shared across all mice 'logistic') w_rec_pt = w3[:opts['n_rec'],:].T.ßatten() APPENDIX A def __init__(self,opts=None): self.opts = self.get_defaults(opts); self.initialize_network() self.rnn = nn def get_defaults(self,opts = None): if opts is None: def get_time_evolution(self,r_in = None,u = opts = {}; None,T = None): if T is None: if 'n_rec' not in opts: T = 1000 opts['n_rec'] = 30 if r_in is None: r_in = np.random.uniform(- 1,1,self.opts['n_rec'])

if u is None: if self.opts['n_in']>0: if 'g' not in opts: u = np.zeros([T,self.opts['n_in']]) opts['g'] = 3 if 'sigma_input' not in opts: if self.opts['n_in']>0: opts['sigma_input'] = 0 w_in = self.rnn['w_in'] if 'sigma_conversion' not in opts: w_rec = self.rnn['w_rec'] opts['sigma_conversion'] = 0 alpha = self.opts['alpha'] sigma_input = self.opts['sigma_input'] if 'alpha' not in opts: sigma_conversion = opts['alpha'] = 0.1 self.opts['sigma_conversion'] if 'input_noise_type' not in opts: n_rec = w_rec.shape[0] opts['input_noise_type'] = 'Gaussian' r = np.zeros([T+1,n_rec]) if 'conversion_noise_type' not in opts: opts['conversion_noise_type'] = 'Gaussian' r[0,:] = r_in for i in range(T): return opts r_temp=r[i,:]; def initialize_network(self): if self.opts['input_noise_type'] == 'Gaussian': n_rec = self.opts['n_rec'] noise_input = n_in = self.opts['n_in'] np.random.normal(0,sigma_input,n_rec) g = self.opts['g'] elif self.opts['input_noise_type'] == 'Laplace': # Initialize the input weights noise_input = nn={} np.random.laplace(0,sigma_input,n_rec) if self.opts['n_in']>0: elif self.opts['input_noise_type'] == 'Poisson': noise_input = nn['w_in']=np.random.normal(0,g/np.sqrt(n_in),n_ np.random.poisson(sigma_input,n_rec) in*n_rec).reshape(n_rec,n_in) # Initialize the recurrent weights if self.opts['n_in']>0: z = np.dot(w_rec,r_temp) + temp=np.random.normal(0,g/np.sqrt(n_rec),n_rec np.dot(w_in,u[i,:]) \ *n_rec).reshape(n_rec,n_rec) + noise_input; np.Þll_diagonal(temp,0) nn['w_rec']=temp APPENDIX A else: z = np.dot(w_rec,r_temp) + noise_input; r_out = alpha*np.tanh(z) + (1-alpha)*r; return r_out if self.opts['conversion_noise_type'] == 'Gaussian': run_experiment_1a.py noise_conversion = #!/usr/bin/env python3 np.random.normal(0,sigma_conversion,n_rec) # -*- coding: utf-8 -*- elif self.opts['conversion_noise_type'] == """ 'Laplace': Created on Wed Dec 2117:52:092022 noise_conversion = np.random.laplace(0,sigma_conversion,n_rec) @author: dinc elif self.opts['conversion_noise_type'] == """ 'Poisson': noise_conversion = import numpy as np np.random.poisson(sigma_conversion,n_rec) from scipy.stats import pearsonr from scipy.stats import spearmanr temp_val = alpha*np.tanh(z) + (1- from rnn_class import RNN alpha)*r_temp + alpha*noise_conversion; import matplotlib.pyplot as plt from numpy.linalg import inv # Newly added error could throw Þring rates from utils import solve_corrn out of bounds, correct for it! from utils import solve_corrn_gpu #temp_val[temp_val >= 1] = 1; from utils import solve_newton_descent #temp_val[temp_val<=-1] = -1; from utils import solve_gradient_descent from utils import solve_pytorch r[i+1,:] = temp_val; from utils import solve_pytorch_gpu from utils import Þt_FORCE import time as time_now return r time_limit = 120; def run_forward_propagation(self,r,u = None): opts = {}; """ opts['g'] = 3; When inputting r, make sure you input it as r[:- opts['n_rec'] = 200 # 300 for short 1,:], as we want to compute the Þring rates opts['n_in'] = 1 at time t+1, but r[T+1,:] is not deÞned, as r is opts['sigma_input'] = 1e-2 deÞned between 0 and T, both inclusive! opts['sigma_conversion'] = 1e-4 """ opts['alpha'] = 0.1 T = r.shape[0] opts['input_noise_type'] = 'Gaussian' if u is not None and u.shape[0] != T: opts['conversion_noise_type'] = 'Gaussian' raise Exception('There is a missmatch opts['verbose'] = False; between u and r dimensions!!!') opts['lambda_reg'] = 1e-5 opts['num_cores'] = 4 if self.opts['n_in']>0: opts['parallel'] = 1 w_in = self.rnn['w_in'] T_data = 3000; # 3000 for short w_rec = self.rnn['w_rec'] alpha = self.opts['alpha'] num_algs = 9; num_exps = 100; if self.opts['n_in']>0: iter_list = np.linspace(0,20,21).astype(int) z = np.dot(r,w_rec.T) + np.dot(u,w_in.T) slopes = else: np.zeros([num_algs,num_exps,iter_list.shape[0]]) z = np.dot(r,w_rec.T) + np.nan APPENDIX A rmse = initialize_fp = 1,num_iters = np.zeros([num_algs,num_exps,iter_list.shape[0]]) iters, + np.nan solver_type = 'weighted') rmedse = times[3,idx_exp,iters] = np.zeros([num_algs,num_exps,iter_list.shape[0]]) time_now.perf_counter() - start_time; + np.nan w_rec = w[:opts['n_rec'],:].T.ßatten() cor_p = slopes[3,idx_exp,iters] = (gnd @ w_rec) / np.zeros([num_algs,num_exps,iter_list.shape[0]]) (gnd@gnd) + np.nan rmse[3,idx_exp,iters] = np.sqrt(np.mean( cor_s = (w_rec-gnd)**2 )) np.zeros([num_algs,num_exps,iter_list.shape[0]]) rmedse[3,idx_exp,iters] = + np.nan np.sqrt(np.median( (w_rec-gnd)**2 )) times = cor_p[3,idx_exp,iters] = np.zeros([num_algs,num_exps,iter_list.shape[0]]) pearsonr(w_rec,gnd)[0] + np.nan cor_s[3,idx_exp,iters] = spearmanr

for idx_exp in range(num_exps): %(cor_p counter = np.zeros(num_algs) if time_now.perf_counter() - start_time > temp = time_now.localtime() time_limit: current_time = time_now.strftime("%H:%M:%S", counter[3] =1 temp) except: print('%s: Running Experiment %d' w = 0 %(current_time,idx_exp)) else: m1 = RNN(opts) print('\t \t \t NT, weighted skipped due to r = m1.get_time_evolution(T = T_data) timeout error' ) gnd = m1.rnn['w_rec'].ßatten() #%% for iters in iter_list: temp = time_now.localtime() current_time = start_time = time_now.perf_counter() time_now.strftime("%H:%M:%S", temp) if counter[5] == 0: print('\t %s: Running iteration %d' try: %(current_time,iters)) w = solve_newton_descent(r[:- 1,:],r[1:,:],u_in = None,alph = 0.1, if iters == 0: l2 = 1e-5,verbose = 0,threshold try: =0.2, w = solve_newton_descent(r[:- initialize_fp = 1,num_iters = 1,:],r[1:,:],u_in = None,alph = 0.1, iters, l2 = 1e-5,verbose = 0,threshold solver_type = 'logistic') =0.2, times[5,idx_exp,iters] = initialize_fp = 1,num_iters = iters, time_now.perf_counter() - start_time; solver_type = 'weighted') w_rec = w[:opts['n_rec'],:].T.ßatten() except: slopes[5,idx_exp,iters] = (gnd @ w_rec) / w = 0 (gnd@gnd) rmse[5,idx_exp,iters] = np.sqrt(np.mean( start_time = time_now.perf_counter() (w_rec-gnd)**2 )) if counter[3] == 0: rmedse[5,idx_exp,iters] = try: np.sqrt(np.median( (w_rec-gnd)**2 )) w = solve_newton_descent(r[:- cor_p[5,idx_exp,iters] = 1,:],r[1:,:],u_in = None,alph = 0.1, pearsonr(w_rec,gnd)[0] l2 = 1e-5,verbose = 0,threshold cor_s[5,idx_exp,iters] = =0.2, spearmanr(w_rec,gnd)[0] APPENDIX A print('\t \t \t NT,logistic. Cor %.3f' else: %(cor_p[5,idx_exp,iters])) print('\t \t \t CoRNN, CPU skipped due to if time_now.perf_counter() - start_time > timeout error') time_limit: counter[5] =1 except: #%% print('\t \t \t NT,logistic terminated due to if iters == 0: timeout error') w = solve_corrn_gpu(r[:-1,:],r[1:,:],u_in = counter[5] = 1 None,l2 = 1*1e-5, else: verbose = 0,gnd = gnd, threshold = print('\t \t \t NT,logistic skipped due to 0.2, timeout error') num_iters = 5*iters, solver_type = 'weighted') start_time = time_now.perf_counter() if counter[8] == 0: try: w = solve_corrn_gpu(r[:-1,:],r[1:,:],u_in = None,l2 = 1*1e-5, #%% verbose = 0,gnd = gnd, threshold = start_time = time_now.perf_counter() 0.2, if counter[0] == 0: num_iters = 5*iters, try: solver_type = 'weighted') times[8,idx_exp,iters] = w = solve_corrn(r[:-1,:],r[1:,:],u_in = time_now.perf_counter() - start_time; None,l2 = 1*1e-5, w_rec = w[:opts['n_rec'],:].T.ßatten() verbose = 0,gnd = gnd, threshold = slopes[8,idx_exp,iters] = (gnd @ w_rec) / 0.2, (gnd@gnd) num_iters = rmse[8,idx_exp,iters] = np.sqrt(np.mean( 5*iters,check_convergence = 0, (w_rec-gnd)**2 )) solver_type = 'weighted') rmedse[8,idx_exp,iters] = times[0,idx_exp,iters] = np.sqrt(np.median( (w_rec-gnd)**2 )) time_now.perf_counter() - start_time; cor_p[8,idx_exp,iters] = w_rec = w[:opts['n_rec'],:].T.ßatten() pearsonr(w_rec,gnd)[0] slopes[0,idx_exp,iters] = (gnd @ w_rec) / cor_s[8,idx_exp,iters] = (gnd@gnd) spearmanr(w_rec,gnd)[0] rmse[0,idx_exp,iters] = np.sqrt(np.mean( print('\t \t \t CoRNN, GPU. Cor %.3f' (w_rec-gnd)**2 )) % [8,idx_exp,iters])) rmedse[0,idx_exp,iters] =

if time_now.perf_counter() - start_time > np.sqrt(np.median( (w_rec-gnd)**2 )) time_limit: cor_p[0,idx_exp,iters] = counter[8] =1 pearsonr(w_rec,gnd)[0] except: cor_s[0,idx_exp,iters] = print('\t \t \t CoRNN, GPU terminated due spearmanr(w_rec,gnd)[0] to timeout error') print('\t \t \t CoRNN, CPU. Cor %.3f' counter[8] = 1 %(cor_p[0,idx_exp,iters])) else: if time_now.perf_counter() - start_time > print('\t \t \t CoRNN, GPU skipped due to time_limit: timeout error') counter[0] =1 except: print('\t \t \t CoRNN, CPU terminated due #%% to timeout error') start_time = time_now.perf_counter() counter[0] = 1 if counter[1] == 0: APPENDIX A try: slopes[2,idx_exp,iters] = (gnd @ w_rec) / (gnd@gnd) w = solve_pytorch_gpu(r[:-1,:],r[1:,:],u_in = rmse[2,idx_exp,iters] = np.sqrt(np.mean( None,alph = 0.1, (w_rec-gnd)**2 )) l2 = 1e-8,verbose = 0,gnd = rmedse[2,idx_exp,iters] = gnd, #1e-7 for short np.sqrt(np.median( (w_rec-gnd)**2 )) initialize_fp =1,num_iters = cor_p[2,idx_exp,iters] = iters*1000, pearsonr(w_rec,gnd)[0] learning_rate = .01,solver_type cor_s[2,idx_exp,iters] = = 'logistic') spearmanr(w_rec,gnd)[0] times[1,idx_exp,iters] = print('\t \t \t Pytorch,logistic, CPU. Cor time_now.perf_counter() - start_time; %.3f' %(cor_p[2,idx_exp,iters])) w_rec = w[:opts['n_rec'],:].T.ßatten() if time_now.perf_counter() - start_time > slopes[1,idx_exp,iters] = (gnd @ w_rec) / time_limit: (gnd@gnd) counter[2] =1 rmse[1,idx_exp,iters] = np.sqrt(np.mean( except: (w_rec-gnd)**2 )) print('\t \t \t Pytorch,logistic, CPU rmedse[1,idx_exp,iters] = terminated due to timeout error') np.sqrt(np.median( (w_rec-gnd)**2 )) counter[2] = 1 cor_p[1,idx_exp,iters] = else: pearsonr(w_rec,gnd)[0] print('\t \t \t Pytorch,logistic, CPU skipped cor_s[1,idx_exp,iters] = due to timeout error') spearmanr(w_rec,gnd)[0] print('\t \t \t Pytorch,logistic,GPU. Cor %.3f' %(cor_p[1,idx_exp,iters])) #%% if time_now.perf_counter() - start_time > start_time = time_now.perf_counter() time_limit: if counter[4] == 0: counter[1] =1 try: except: print('\t \t \t Pytorch,logistic,GPU w = solve_pytorch_gpu(r[:-1,:],r[1:,:],u_in = terminated due to timeout error') None,alph = 0.1, counter[1] = 1 l2 = 1e-8,verbose = 0,gnd = else: gnd, #1e-7 for short print('\t \t \t Pytorch,logistic,GPU skipped initialize_fp =1,num_iters = due to timeout error') iters*1000, learning_rate = .01,solver_type = 'l2') #%% times[4,idx_exp,iters] = start_time = time_now.perf_counter() time_now.perf_counter() - start_time; if counter[2] == 0: w_rec = w[:opts['n_rec'],:].T.ßatten() try: slopes[4,idx_exp,iters] = (gnd @ w_rec) / (gnd@gnd) w = solve_pytorch(r[:-1,:],r[1:,:],u_in = rmse[4,idx_exp,iters] = np.sqrt(np.mean( None,alph = 0.1, (w_rec-gnd)**2 )) l2 = 1e-8,verbose = 0,gnd = rmedse[4,idx_exp,iters] = gnd, np.sqrt(np.median( (w_rec-gnd)**2 )) initialize_fp = 1,num_iters = cor_p[4,idx_exp,iters] = iters * 1000, pearsonr(w_rec,gnd)[0] learning_rate = .01,solver_type cor_s[4,idx_exp,iters] = = 'logistic') spearmanr(w_rec,gnd)[0] times[2,idx_exp,iters] = print('\t \t \t Pytorch,l2,GPU. Cor %.3f' time_now.perf_counter() - start_time; %(cor_p[4,idx_exp,iters])) w_rec = w[:opts['n_rec'],:].T.ßatten() APPENDIX A if time_now.perf_counter() - start_time > w = Þt_FORCE(r,None,alph = 0.1, time_limit: lam = 200,g_in = 2,verbose = 0, counter[4] =1 initialize_fp = 1,num_iters = 5*iters, except: gnd = gnd,solver_type = print('\t \t \t Pytorch,l2,GPU terminated 'Þring_rates') due to timeout error') times[7,idx_exp,iters] = counter[4] = 1 time_now.perf_counter() - start_time; else: w_rec = w[:opts['n_rec'],:].T.ßatten() print('\t \t \t Pytorch,l2,GPU skipped due to slopes[7,idx_exp,iters] = (gnd @ w_rec) / timeout error') (gnd@gnd) rmse[7,idx_exp,iters] = np.sqrt(np.mean( (w_rec-gnd)**2 )) #%% rmedse[7,idx_exp,iters] = start_time = time_now.perf_counter() np.sqrt(np.median( (w_rec-gnd)**2 )) if counter[6] == 0: cor_p[7,idx_exp,iters] = try: pearsonr(w_rec,gnd)[0] cor_s[7,idx_exp,iters] = w = Þt_FORCE(r,None,alph = 0.1, spearmanr(w_rec,gnd)[0] lam = 200,g_in = 2,verbose = 0, print('\t \t \t Force, Þring rates. Cor %.3f' initialize_fp = 1,num_iters = 5*iters, %(cor_p[7,idx_exp,iters])) gnd = gnd,solver_type = 'currents') if time_now.perf_counter() - start_time > times[6,idx_exp,iters] = time_limit: time_now.perf_counter() - start_time; counter[7] =1 w_rec = w[:opts['n_rec'],:].T.ßatten() except: slopes[6,idx_exp,iters] = (gnd @ w_rec) / print('\t \t \t Force, Þring rates terminated (gnd@gnd) due to timeout error') rmse[6,idx_exp,iters] = np.sqrt(np.mean( counter[7] = 1 (w_rec-gnd)**2 )) else: rmedse[6,idx_exp,iters] = print('\t \t \t Force, Þring rates skipped due to np.sqrt(np.median( (w_rec-gnd)**2 )) timeout error') cor_p[6,idx_exp,iters] = pearsonr(w_rec,gnd)[0] if iters == 0: cor_s[6,idx_exp,iters] = counter = np.zeros(num_algs) spearmanr(w_rec,gnd)[0] print('\t \t \t Force,currents. Cor %.3f' %(cor_p[6,idx_exp,iters])) np.savez('experiment1a_results_fp.npz',slopes = if time_now.perf_counter() - start_time > slopes, time_limit: rmse = rmse, rmedse = rmedse,cor_p = counter[6] =1 cor_p,cor_s = cor_s, except: times = times) print('\t \t \t Force,currents terminated due to timeout error') counter[6] = 1 run_experiment_1b.py else: #!/usr/bin/env python3 print('\t \t \t Force,currents skipped due to # -*- coding: utf-8 -*- timeout error') """ Created on Wed Dec 2117:52:092022 #%% @author: dinc start_time = time_now.perf_counter() """ if counter[7] == 0: try: import numpy as np from scipy.stats import pearsonr APPENDIX A from scipy.stats import spearmanr from rnn_class import RNN for idx_exp in range(num_exps): import matplotlib.pyplot as plt counter = np.zeros(num_algs) from numpy.linalg import inv temp = time_now.localtime() from utils import solve_corrn current_time = time_now.strftime("%H:%M:%S", from utils import solve_corrn_gpu temp) from utils import solve_newton_descent print('%s: Running Experiment %d' from utils import solve_gradient_descent %(current_time,idx_exp)) from utils import solve_pytorch m1 = RNN(opts) from utils import solve_pytorch_gpu r = m1.get_time_evolution(T = T_data) from utils import Þt_FORCE gnd = m1.rnn['w_rec'].ßatten() import time as time_now #% for iters in iter_list: time_limit = 120; temp = time_now.localtime() current_time = opts = {}; time_now.strftime("%H:%M:%S", temp) opts['g'] = 3; print('\t %s: Running iteration %d' opts['n_rec'] = 200 # 300 for short %(current_time,iters)) opts['n_in'] = 1 opts['sigma_input'] = 1e-1 if iters == 0: opts['sigma_conversion'] = 1e-5 try: opts['alpha'] = 0.1 w = solve_newton_descent(r[:- opts['input_noise_type'] = 'Gaussian' 1,:],r[1:,:],u_in = None,alph = 0.1, opts['conversion_noise_type'] = 'Gaussian' l2 = 1e-5,verbose = 0,threshold opts['verbose'] = False; =0.2, opts['lambda_reg'] = 1e-5 initialize_fp = 1,num_iters = iters, opts['num_cores'] = 4 solver_type = 'weighted') opts['parallel'] = 1 except: T_data = 3000; # 3000 for short w = 0 start_time = time_now.perf_counter() num_algs = 9; if counter[3] == 0: num_exps = 100; try: iter_list = np.linspace(0,20,21).astype(int) w = solve_newton_descent(r[:- slopes = 1,:],r[1:,:],u_in = None,alph = 0.1, np.zeros([num_algs,num_exps,iter_list.shape[0]]) l2 = 1e-5,verbose = 0,threshold + np.nan =0.2, rmse = initialize_fp = 1,num_iters = np.zeros([num_algs,num_exps,iter_list.shape[0]]) iters, + np.nan solver_type = 'weighted') rmedse = times[3,idx_exp,iters] = np.zeros([num_algs,num_exps,iter_list.shape[0]]) time_now.perf_counter() - start_time; + np.nan w_rec = w[:opts['n_rec'],:].T.ßatten() cor_p = slopes[3,idx_exp,iters] = (gnd @ w_rec) / np.zeros([num_algs,num_exps,iter_list.shape[0]]) (gnd@gnd) + np.nan rmse[3,idx_exp,iters] = np.sqrt(np.mean( cor_s = (w_rec-gnd)**2 )) np.zeros([num_algs,num_exps,iter_list.shape[0]]) rmedse[3,idx_exp,iters] = + np.nan np.sqrt(np.median( (w_rec-gnd)**2 )) times = cor_p[3,idx_exp,iters] = np.zeros([num_algs,num_exps,iter_list.shape[0]]) pearsonr(w_rec,gnd)[0] + np.nan spearmanr

APPENDIX A print('\t \t \t NT,weighted. Cor %.3f' %(cor_p[3,idx_exp,iters])) if time_now.perf_counter() - start_time > time_limit: #% counter[3] =1 start_time = time_now.perf_counter() except: if counter[0] == 0: w = 0 try: else: print('\t \t \t NT, weighted skipped due to w = solve_corrn(r[:-1,:],r[1:,:],u_in = timeout error' ) None,l2 = 1*1e-5, verbose = 0,gnd = gnd, threshold = 0.2, num_iters = 5*iters,check_convergence = 0, start_time = time_now.perf_counter() solver_type = 'weighted') if counter[5] == 0: times[0,idx_exp,iters] = try: time_now.perf_counter() - start_time; w = solve_newton_descent(r[:- w_rec = w[:opts['n_rec'],:].T.ßatten() 1,:],r[1:,:],u_in = None,alph = 0.1, slopes[0,idx_exp,iters] = (gnd @ w_rec) / l2 = 1e-5,verbose = 0,threshold (gnd@gnd) =0.2, rmse[0,idx_exp,iters] = np.sqrt(np.mean( initialize_fp = 1,num_iters = (w_rec-gnd)**2 )) iters, rmedse[0,idx_exp,iters] = solver_type = 'logistic') np.sqrt(np.median( (w_rec-gnd)**2 )) times[5,idx_exp,iters] = cor_p[0,idx_exp,iters] = time_now.perf_counter() - start_time; pearsonr(w_rec,gnd)[0] w_rec = w[:opts['n_rec'],:].T.ßatten() cor_s[0,idx_exp,iters] = slopes[5,idx_exp,iters] = (gnd @ w_rec) / spearmanr(w_rec,gnd)[0] (gnd@gnd) print('\t \t \t CoRNN, CPU. Cor %.3f' rmse[5,idx_exp,iters] = np.sqrt(np.mean( %(cor_p[0,idx_exp,iters])) (w_rec-gnd)**2 )) if time_now.perf_counter() - start_time > rmedse[5,idx_exp,iters] = time_limit: np.sqrt(np.median( (w_rec-gnd)**2 )) counter[0] =1 cor_p[5,idx_exp,iters] = except: pearsonr(w_rec,gnd)[0] print('\t \t \t CoRNN, CPU terminated due cor_s[5,idx_exp,iters] = to timeout error') spearmanr(w_rec,gnd)[0] counter[0] = 1 print('\t \t \t NT,logistic. Cor %.3f' else: %(cor_p[5,idx_exp,iters])) print('\t \t \t CoRNN, CPU skipped due to if time_now.perf_counter() - start_time > timeout error') time_limit: counter[5] =1 except: #% print('\t \t \t NT,logistic terminated due to if iters == 0: timeout error') w = solve_corrn_gpu(r[:-1,:],r[1:,:],u_in = counter[5] = 1 None,l2 = 1*1e-5, else: verbose = 0,gnd = gnd, threshold = print('\t \t \t NT,logistic skipped due to 0.2, timeout error') num_iters = 5*iters, solver_type = 'weighted') start_time = time_now.perf_counter() if counter[8] == 0: try: APPENDIX A rmse[1,idx_exp,iters] = np.sqrt(np.mean( w = solve_corrn_gpu(r[:-1,:],r[1:,:],u_in = (w_rec-gnd)**2 )) None,l2 = 1*1e-5, rmedse[1,idx_exp,iters] = verbose = 0,gnd = gnd, threshold = np.sqrt(np.median( (w_rec-gnd)**2 )) 0.2, cor_p[1,idx_exp,iters] = num_iters = 5*iters, pearsonr(w_rec,gnd)[0] solver_type = 'weighted') cor_s[1,idx_exp,iters] = times[8,idx_exp,iters] = spearmanr(w_rec,gnd)[0] time_now.perf_counter() - start_time; print('\t \t \t Pytorch,logistic,GPU. Cor w_rec = w[:opts['n_rec'],:].T.ßatten() %.3f' %(cor_p[1,idx_exp,iters])) slopes[8,idx_exp,iters] = (gnd @ w_rec) / if time_now.perf_counter() - start_time > (gnd@gnd) time_limit: rmse[8,idx_exp,iters] = np.sqrt(np.mean( counter[1] =1 (w_rec-gnd)**2 )) except: rmedse[8,idx_exp,iters] = print('\t \t \t Pytorch,logistic,GPU np.sqrt(np.median( (w_rec-gnd)**2 )) terminated due to timeout error') cor_p[8,idx_exp,iters] = counter[1] = 1 pearsonr(w_rec,gnd)[0] else: cor_s[8,idx_exp,iters] = print('\t \t \t Pytorch,logistic,GPU skipped spearmanr(w_rec,gnd)[0] due to timeout error') print('\t \t \t CoRNN, GPU. Cor %.3f' %(cor_p[8,idx_exp,iters])) if time_now.perf_counter() - start_time > #% time_limit: start_time = time_now.perf_counter() counter[8] =1 if counter[2] == 0: except: try: print('\t \t \t CoRNN, GPU terminated due to timeout error') w = solve_pytorch(r[:-1,:],r[1:,:],u_in = counter[8] = 1 None,alph = 0.1, else: l2 = 1e-8,verbose = 0,gnd = print('\t \t \t CoRNN, GPU skipped due to gnd, timeout error') initialize_fp = 1,num_iters = iters * 1000, learning_rate = .01,solver_type #% = 'logistic') start_time = time_now.perf_counter() times[2,idx_exp,iters] = if counter[1] == 0: time_now.perf_counter() - start_time; try: w_rec = w[:opts['n_rec'],:].T.ßatten() slopes[2,idx_exp,iters] = (gnd @ w_rec) / w = solve_pytorch_gpu(r[:-1,:],r[1:,:],u_in = (gnd@gnd) None,alph = 0.1, rmse[2,idx_exp,iters] = np.sqrt(np.mean( l2 = 1e-8,verbose = 0,gnd = (w_rec-gnd)**2 )) gnd, #1e-7 for short rmedse[2,idx_exp,iters] = initialize_fp =1,num_iters = np.sqrt(np.median( (w_rec-gnd)**2 )) iters*1000, cor_p[2,idx_exp,iters] = learning_rate = .01,solver_type pearsonr(w_rec,gnd)[0] = 'logistic') cor_s[2,idx_exp,iters] = times[1,idx_exp,iters] = spearmanr(w_rec,gnd)[0] time_now.perf_counter() - start_time; print('\t \t \t Pytorch,logistic, CPU. Cor w_rec = w[:opts['n_rec'],:].T.ßatten() %.3f' %(cor_p[2,idx_exp,iters])) slopes[1,idx_exp,iters] = (gnd @ w_rec) / if time_now.perf_counter() - start_time > (gnd@gnd) time_limit: counter[2] =1 APPENDIX A except: print('\t \t \t Pytorch,logistic, CPU w = Þt_FORCE(r,None,alph = 0.1, terminated due to timeout error') lam = 200,g_in = 2,verbose = 0, counter[2] = 1 initialize_fp = 1,num_iters = 5*iters, else: gnd = gnd,solver_type = 'currents') print('\t \t \t Pytorch,logistic, CPU skipped times[6,idx_exp,iters] = due to timeout error') time_now.perf_counter() - start_time; w_rec = w[:opts['n_rec'],:].T.ßatten() slopes[6,idx_exp,iters] = (gnd @ w_rec) / #% (gnd@gnd) start_time = time_now.perf_counter() rmse[6,idx_exp,iters] = np.sqrt(np.mean( if counter[4] == 0: (w_rec-gnd)**2 )) try: rmedse[6,idx_exp,iters] = np.sqrt(np.median( (w_rec-gnd)**2 )) w = solve_pytorch_gpu(r[:-1,:],r[1:,:],u_in = cor_p[6,idx_exp,iters] = None,alph = 0.1, pearsonr(w_rec,gnd)[0] l2 = 1e-8,verbose = 0,gnd = cor_s[6,idx_exp,iters] = gnd, #1e-7 for short spearmanr(w_rec,gnd)[0] initialize_fp =1,num_iters = print('\t \t \t Force,currents. Cor %.3f' iters*1000, %(cor_p[6,idx_exp,iters])) learning_rate = .01,solver_type if time_now.perf_counter() - start_time > = 'l2') time_limit: times[4,idx_exp,iters] = counter[6] =1 time_now.perf_counter() - start_time; except: w_rec = w[:opts['n_rec'],:].T.ßatten() print('\t \t \t Force,currents terminated due slopes[4,idx_exp,iters] = (gnd @ w_rec) / to timeout error') (gnd@gnd) counter[6] = 1 rmse[4,idx_exp,iters] = np.sqrt(np.mean( else: (w_rec-gnd)**2 )) print('\t \t \t Force,currents skipped due to rmedse[4,idx_exp,iters] = timeout error') np.sqrt(np.median( (w_rec-gnd)**2 )) cor_p[4,idx_exp,iters] = pearsonr(w_rec,gnd)[0] #% cor_s[4,idx_exp,iters] = start_time = time_now.perf_counter() spearmanr(w_rec,gnd)[0] if counter[7] == 0: print('\t \t \t Pytorch,l2,GPU. Cor %.3f' try: %(cor_p[4,idx_exp,iters])) if time_now.perf_counter() - start_time > w = Þt_FORCE(r,None,alph = 0.1, time_limit: lam = 200,g_in = 2,verbose = 0, counter[4] =1 initialize_fp = 1,num_iters = 5*iters, except: gnd = gnd,solver_type = print('\t \t \t Pytorch,l2,GPU terminated 'Þring_rates') due to timeout error') times[7,idx_exp,iters] = counter[4] = 1 time_now.perf_counter() - start_time; else: w_rec = w[:opts['n_rec'],:].T.ßatten() print('\t \t \t Pytorch,l2,GPU skipped due to slopes[7,idx_exp,iters] = (gnd @ w_rec) / timeout error') (gnd@gnd) rmse[7,idx_exp,iters] = np.sqrt(np.mean( (w_rec-gnd)**2 )) #% rmedse[7,idx_exp,iters] = start_time = time_now.perf_counter() np.sqrt(np.median( (w_rec-gnd)**2 )) if counter[6] == 0: cor_p[7,idx_exp,iters] = try: pearsonr(w_rec,gnd)[0] APPENDIX A cor_s[7,idx_exp,iters] = opts = {}; spearmanr(w_rec,gnd)[0] opts['g'] = 3; print('\t \t \t Force, Þring rates. Cor %.3f' opts['n_rec'] = 500 # 300 for short %(cor_p[7,idx_exp,iters])) opts['n_in'] = 1 if time_now.perf_counter() - start_time > opts['sigma_input'] = 1e-1 time_limit: opts['sigma_conversion'] = 1e-5 counter[7] =1 opts['alpha'] = 0.1 except: opts['input_noise_type'] = 'Gaussian' print('\t \t \t Force, Þring rates terminated opts['conversion_noise_type'] = 'Gaussian' due to timeout error') opts['verbose'] = False; counter[7] = 1 opts['lambda_reg'] = 1e-5 else: opts['num_cores'] = 4 print('\t \t \t Force, Þring rates skipped due to opts['parallel'] = 1 timeout error') T_data = 10000; # 3000 for short if iters == 0: counter = np.zeros(num_algs) num_algs = 9; num_exps = 100; iter_list = np.linspace(0,10,11).astype(int) np.savez('experiment1b_results_fp.npz',slopes = slopes = slopes, np.zeros([num_algs,num_exps,iter_list.shape[0]]) rmse = rmse, rmedse = rmedse,cor_p = + np.nan cor_p,cor_s = cor_s, rmse = times = times) np.zeros([num_algs,num_exps,iter_list.shape[0]]) + np.nan rmedse = run_experiment_1c.py np.zeros([num_algs,num_exps,iter_list.shape[0]]) + np.nan #!/usr/bin/env python3 cor_p = # -*- coding: utf-8 -*- np.zeros([num_algs,num_exps,iter_list.shape[0]]) """ + np.nan Created on Wed Dec 2117:52:092022 cor_s = np.zeros([num_algs,num_exps,iter_list.shape[0]]) @author: dinc + np.nan """ times = np.zeros([num_algs,num_exps,iter_list.shape[0]]) import numpy as np + np.nan from scipy.stats import pearsonr from scipy.stats import spearmanr from rnn_class import RNN for idx_exp in range(num_exps): import matplotlib.pyplot as plt counter = np.zeros(num_algs) from numpy.linalg import inv temp = time_now.localtime() from utils import solve_corrn current_time = time_now.strftime("%H:%M:%S", from utils import solve_corrn_gpu temp) from utils import solve_newton_descent print('%s: Running Experiment %d' from utils import solve_gradient_descent %(current_time,idx_exp)) from utils import solve_pytorch m1 = RNN(opts) from utils import solve_pytorch_gpu r = m1.get_time_evolution(T = T_data) from utils import Þt_FORCE gnd = import time as time_now #%

for iters time_limit = 200; temp = time_now.localtime() APPENDIX A current_time = time_now.strftime("%H:%M:%S", temp) print('\t %s: Running iteration %d' start_time = time_now.perf_counter() %(current_time,iters)) if counter[5] == 0: try: if iters == 0: w = solve_newton_descent(r[:- try: 1,:],r[1:,:],u_in = None,alph = 0.1, w = solve_newton_descent(r[:- l2 = 1e-5,verbose = 0,threshold 1,:],r[1:,:],u_in = None,alph = 0.1, =1, l2 = 1e-5,verbose = 0,threshold initialize_fp = 0,num_iters = =0.2, iters, initialize_fp = 0,num_iters = iters, solver_type = 'logistic') solver_type = 'weighted') times[5,idx_exp,iters] = except: time_now.perf_counter() - start_time; w = 0 w_rec = w[:opts['n_rec'],:].T.ßatten() slopes[5,idx_exp,iters] = (gnd @ w_rec) / start_time = time_now.perf_counter() (gnd@gnd) if counter[3] == 0: rmse[5,idx_exp,iters] = np.sqrt(np.mean( try: (w_rec-gnd)**2 )) w = solve_newton_descent(r[:- rmedse[5,idx_exp,iters] = 1,:],r[1:,:],u_in = None,alph = 0.1, np.sqrt(np.median( (w_rec-gnd)**2 )) l2 = 1e-5,verbose = 0,threshold cor_p[5,idx_exp,iters] = =1, pearsonr(w_rec,gnd)[0] initialize_fp = 0,num_iters = cor_s[5,idx_exp,iters] = iters, spearmanr(w_rec,gnd)[0] solver_type = 'weighted') _print('\t \t \t NT,logistic. Cor %.3f. %.1f times[3,idx_exp,iters] = mins' %(cor_p[3,idx_exp,iters],\ time_now.perf_counter() - start_time; ( time_now.perf_counter() - w_rec = w[:opts['n_rec'],:].T.ßatten() start_time)/60 )) slopes[3,idx_exp,iters] = (gnd @ w_rec) / if time_now.perf_counter() - start_time > (gnd@gnd) time_limit: rmse[3,idx_exp,iters] = np.sqrt(np.mean( counter[5] =1 (w_rec-gnd)**2 )) except: rmedse[3,idx_exp,iters] = print('\t \t \t NT,logistic terminated due to np.sqrt(np.median( (w_rec-gnd)**2 )) timeout error') cor_p[3,idx_exp,iters] = counter[5] = 1 pearsonr(w_rec,gnd)[0] else: cor_s[3,idx_exp,iters] = print('\t \t \t NT,logistic skipped due to spearmanr(w_rec,gnd)[0] timeout error') print('\t \t \t NT,weighted. Cor %.3f. %.1f mins' %(cor_p[3,idx_exp,iters],\ ( time_now.perf_counter() - start_time)/60 )) if time_now.perf_counter() - start_time > time_limit: counter[3] =1 except: #% w = 0 start_time = time_now.perf_counter() else: if counter[0] == 0: print('\t \t \t NT, weighted skipped due to try: timeout error' ) w = solve_corrn(r[:-1,:],r[1:,:],u_in = None,l2 = 1*1e-5, APPENDIX A verbose = 0,gnd = gnd, threshold = w_rec = w[:opts['n_rec'],:].T.ßatten() 1, initialize_fp = 0, slopes[8,idx_exp,iters] = (gnd @ w_rec) / num_iters = (gnd@gnd) 300*iters,check_convergence = 0, rmse[8,idx_exp,iters] = np.sqrt(np.mean( solver_type = 'weighted') (w_rec-gnd)**2 )) times[0,idx_exp,iters] = rmedse[8,idx_exp,iters] = time_now.perf_counter() - start_time; np.sqrt(np.median( (w_rec-gnd)**2 )) w_rec = w[:opts['n_rec'],:].T.ßatten() cor_p[8,idx_exp,iters] = slopes[0,idx_exp,iters] = (gnd @ w_rec) / pearsonr(w_rec,gnd)[0] (gnd@gnd) cor_s[8,idx_exp,iters] = rmse[0,idx_exp,iters] = np.sqrt(np.mean( spearmanr(w_rec,gnd)[0] (w_rec-gnd)**2 )) print('\t \t \t CoRNN, GPU. Cor %.3f' rmedse[0,idx_exp,iters] = %(cor_p[8,idx_exp,iters])) np.sqrt(np.median( (w_rec-gnd)**2 )) if time_now.perf_counter() - start_time > cor_p[0,idx_exp,iters] = time_limit: pearsonr(w_rec,gnd)[0] counter[8] =1 cor_s[0,idx_exp,iters] = except: spearmanr(w_rec,gnd)[0] print('\t \t \t CoRNN, GPU terminated due print('\t \t \t CoRNN, CPU. Cor %.3f' to timeout error') %(cor_p[0,idx_exp,iters])) counter[8] = 1 if time_now.perf_counter() - start_time > else: time_limit: print('\t \t \t CoRNN, GPU skipped due to counter[0] =1 timeout error') except: print('\t \t \t CoRNN, CPU terminated due to timeout error') #% counter[0] = 1 start_time = time_now.perf_counter() else: if counter[6] == 0: print('\t \t \t CoRNN, CPU skipped due to try: timeout error') w = solve_corrn(r[:-1,:],r[1:,:],u_in = None,l2 = 1*1e-5, #% verbose = 0,gnd = gnd, threshold = if iters == 0: 1, initialize_fp = 1, w = solve_corrn_gpu(r[:-1,:],r[1:,:],u_in = num_iters = None,l2 = 1*1e-5, 5*iters,check_convergence = 0, verbose = 0,gnd = gnd, threshold = 1, solver_type = 'weighted') num_iters = 500*iters, ßoat_type = times[6,idx_exp,iters] = '64bit', time_now.perf_counter() - start_time; solver_type = 'weighted') w_rec = w[:opts['n_rec'],:].T.ßatten() start_time = time_now.perf_counter() slopes[6,idx_exp,iters] = (gnd @ w_rec) / if counter[8] == 0: (gnd@gnd) try: rmse[6,idx_exp,iters] = np.sqrt(np.mean( (w_rec-gnd)**2 )) w = solve_corrn_gpu(r[:-1,:],r[1:,:],u_in = rmedse[6,idx_exp,iters] = None,l2 = 1*1e-5, np.sqrt(np.median( (w_rec-gnd)**2 )) verbose = 0,gnd = gnd, threshold = cor_p[6,idx_exp,iters] = 1, initialize_fp = 0, pearsonr(w_rec,gnd)[0] num_iters = 300*iters, ßoat_type = cor_s[6,idx_exp,iters] = '32bit', spearmanr(w_rec,gnd)[0] solver_type = 'weighted') print('\t \t \t CoRNN, CPU. Cor %.3f' times[8,idx_exp,iters] = %(cor_p[6,idx_exp,iters])) time_now.perf_counter() - start_time; APPENDIX A if time_now.perf_counter() - start_time > else: time_limit: print('\t \t \t CoRNN, GPU skipped due to counter[6] =1 timeout error') except: print('\t \t \t CoRNN, CPU terminated due to timeout error') counter[6] = 1 #% else: start_time = time_now.perf_counter() print('\t \t \t CoRNN, CPU skipped due to if counter[1] == 0: timeout error') try: w = solve_pytorch_gpu(r[:-1,:],r[1:,:],u_in = #% None,alph = 0.1, if iters == 0: l2 = 1e-8,verbose = 0,gnd = w = solve_corrn_gpu(r[:-1,:],r[1:,:],u_in = gnd, #1e-7 for short None,l2 = 1*1e-5, initialize_fp =0,num_iters = verbose = 0,gnd = gnd, threshold = 1, iters*500, num_iters = 500*iters, ßoat_type = learning_rate = .01,solver_type '64bit', = 'logistic') solver_type = 'weighted') times[1,idx_exp,iters] = start_time = time_now.perf_counter() time_now.perf_counter() - start_time; if counter[7] == 0: w_rec = w[:opts['n_rec'],:].T.ßatten() try: slopes[1,idx_exp,iters] = (gnd @ w_rec) / (gnd@gnd) w = solve_corrn_gpu(r[:-1,:],r[1:,:],u_in = rmse[1,idx_exp,iters] = np.sqrt(np.mean( None,l2 = 1*1e-5, (w_rec-gnd)**2 )) verbose = 0,gnd = gnd, threshold = rmedse[1,idx_exp,iters] = 1, initialize_fp = 1, np.sqrt(np.median( (w_rec-gnd)**2 )) num_iters = 5*iters, ßoat_type = cor_p[1,idx_exp,iters] = '32bit', pearsonr(w_rec,gnd)[0] solver_type = 'weighted') cor_s[1,idx_exp,iters] = times[7,idx_exp,iters] = spearmanr(w_rec,gnd)[0] time_now.perf_counter() - start_time; print('\t \t \t Pytorch,logistic,GPU. Cor w_rec = w[:opts['n_rec'],:].T.ßatten() %.3f' %(cor_p[1,idx_exp,iters])) = @ w_rec) / if time_now.perf_counter() - start_time > time_limit:

( counter[1] =1 (w_rec-gnd)**2 )) except: rmedse[7,idx_exp,iters] = print('\t \t \t Pytorch,logistic,GPU np.sqrt(np.median( (w_rec-gnd)**2 )) terminated due to timeout error') cor_p[7,idx_exp,iters] = counter[1] = 1 pearsonr(w_rec,gnd)[0] else: cor_s[7,idx_exp,iters] = print('\t \t \t Pytorch,logistic,GPU skipped spearmanr(w_rec,gnd)[0] due to timeout error') print('\t \t \t CoRNN, GPU. Cor %.3f' %(cor_p[7,idx_exp,iters])) if time_now.perf_counter() - start_time > #% time_limit: start_time = time_now.perf_counter() counter[7] =1 if counter[2] == 0: except: try: print('\t \t \t CoRNN, GPU terminated due to timeout error') w = solve_pytorch(r[:-1,:],r[1:,:],u_in = counter[7] = 1 None,alph = 0.1, APPENDIX A l2 = 1e-8,verbose = 0,gnd = rmedse[4,idx_exp,iters] = gnd, np.sqrt(np.median( (w_rec-gnd)**2 )) initialize_fp = 0,num_iters = cor_p[4,idx_exp,iters] = iters *500, pearsonr(w_rec,gnd)[0] learning_rate = .01,solver_type cor_s[4,idx_exp,iters] = = 'logistic') spearmanr(w_rec,gnd)[0] times[2,idx_exp,iters] = print('\t \t \t Pytorch,l2,GPU. Cor %.3f' time_now.perf_counter() - start_time; %(cor_p[4,idx_exp,iters])) w_rec = w[:opts['n_rec'],:].T.ßatten() if time_now.perf_counter() - start_time > slopes[2,idx_exp,iters] = (gnd @ w_rec) / time_limit: (gnd@gnd) counter[4] =1 rmse[2,idx_exp,iters] = np.sqrt(np.mean( except: (w_rec-gnd)**2 )) print('\t \t \t Pytorch,l2,GPU terminated rmedse[2,idx_exp,iters] = due to timeout error') np.sqrt(np.median( (w_rec-gnd)**2 )) counter[4] = 1 cor_p[2,idx_exp,iters] = else: pearsonr(w_rec,gnd)[0] print('\t \t \t Pytorch,l2,GPU skipped due to cor_s[2,idx_exp,iters] = timeout error') spearmanr(w_rec,gnd)[0] print('\t \t \t Pytorch,logistic, CPU. Cor %.3f' %(cor_p[2,idx_exp,iters])) if time_now.perf_counter() - start_time > if iters == 0: time_limit: counter = np.zeros(num_algs) # Force takes counter[2] =1 hours to run iteration 1.... except: print('\t \t \t Pytorch,logistic, CPU terminated due to timeout error') np.savez('experiment1c_results_fp.npz',slopes = counter[2] = 1 slopes, else: rmse = rmse, rmedse = rmedse,cor_p = print('\t \t \t Pytorch,logistic, CPU skipped cor_p,cor_s = cor_s, due to timeout error') times = times) #% try_timeout.py start_time = time_now.perf_counter() #!/usr/bin/env python3 if counter[4] == 0: # -*- coding: utf-8 -*- try: """ Created on Wed Dec 2117:52:092022 w = solve_pytorch_gpu(r[:-1,:],r[1:,:],u_in = None,alph = 0.1, @author: dinc l2 = 1e-8,verbose = 0,gnd = """ gnd, #1e-7 for short initialize_fp =0,num_iters = import numpy as np iters*500, from scipy.stats import pearsonr learning_rate = .01,solver_type from rnn_class import RNN = 'l2') import matplotlib.pyplot as plt times[4,idx_exp,iters] = from numpy.linalg import inv time_now.perf_counter() - start_time; from utils import solve_corrn_gpu w_rec = w[:opts['n_rec'],:].T.ßatten() from utils import solve_corrn slopes[4,idx_exp,iters] = (gnd @ w_rec) / from utils import solve_newton_descent (gnd@gnd) from utils import solve_gradient_descent rmse[4,idx_exp,iters] = np.sqrt(np.mean( from utils import solve_pytorch_gpu (w_rec-gnd)**2 )) from utils import solve_pytorch APPENDIX A from utils import Þt_FORCE print('function terminated') from utils import Þt_FORCE_gpu import signal import time as time_now #%% signal.alarm(TIMEOUT) try: class TimeoutException(Exception): w1 = solve_newton_descent(r[:-1,:],r[1:,:],u_in = pass None,alph = 0.1, l2 = 1e-5,verbose = 1,threshold =0.2, def timeout_handler(signum, frame): initialize_fp = 0,num_iters = 10, raise TimeoutException solver_type = 'weighted') w_rec_nt = w1[:opts['n_rec'],:].T.ßatten() signal.signal(signal.SIGALRM, timeout_handler) except TimeoutException: print('function terminated') TIMEOUT = 10 # seconds signal.alarm(TIMEOUT) #%% signal.alarm(TIMEOUT) opts = {}; try: opts['g'] = 3; w4 = Þt_FORCE_gpu(r,None,alph = 0.1, opts['n_rec'] = 1000 lam = 100,g_in = 2,verbose = 2, opts['n_in'] = 0 initialize_fp = 1,num_iters = 100, opts['sigma_input'] = 1e-2 gnd = gnd,solver_type = 'Þring_rates') opts['sigma_conversion'] = 1e-4 T_data = 3000; w_rec_fc = w4[:opts['n_rec'],:].T.ßatten() m1 = RNN(opts) signal.alarm(0) r = m1.get_time_evolution(T = T_data) except TimeoutException: gnd = m1.rnn['w_rec'].ßatten() print('function terminated') #%% #%% signal.alarm(TIMEOUT) plt.scatter(gnd,w_rec_pt) try: plt.scatter(gnd,w_rec_co) w = solve_corrn_gpu(r[:-1,:],r[1:,:],u_in = None,l2 = 1*1e-5, print('NT:', pearsonr(w_rec_nt,gnd)[0]) verbose = 2,gnd = gnd, threshold = 0.2, print('GD:', pearsonr(w_rec_pt,gnd)[0]) num_iters = 10,check_convergence = 0, print('CoRNN:', pearsonr(w_rec_co,gnd)[0]) solver_type = 'weighted',alph = 0.1) w_rec_co = w[:opts['n_rec'],:].T.ßatten() except TimeoutException: print('function terminated') utils.py #!/usr/bin/env python3 #%% # -*- coding: utf-8 -*- signal.alarm(TIMEOUT) """ try: Created on Fri Jun 1014:48:402022 w3 = solve_pytorch_gpu(r[:-1,:],r[1:,:],u_in = None,alph = 0.1, @author: dinc l2 = 1e-8,verbose = 2,gnd = gnd, """ initialize_fp = 1,num_iters = 5000, import numpy as np learning_rate = 0.01,solver_type = from scipy.linalg import solve 'logistic') import re w_rec_pt = w3[:opts['n_rec'],:].T.ßatten() import time as time_now except TimeoutException: from numpy.linalg import inv APPENDIX A from scipy.stats import pearsonr Xm = Xp @ x; import multiprocessing z = np.arctanh(d); from joblib import Parallel, delayed theta_fp = Xp @ z; import torch if initialize_fp: import torch.nn as nn theta = theta_fp; import torch.optim as optim else: theta = def solve_corrn(r_in,r_out,u_in = None, alph = 0.1, np.random.normal(0,1/np.sqrt(n_rec),n_rec * l2 = 1e-4, threshold = 1, initialize_fp = 1, n_tot).reshape(n_tot,n_rec) verbose = 0,check_convergence = 0, num_iters = 30,gnd = None,solver_type = if verbose == 2: 'weighted'): temp = time_now.localtime() if verbose: current_time = temp = time_now.localtime() time_now.strftime("%H:%M:%S", temp) current_time = print('\t %s: Starting iterations' time_now.strftime("%H:%M:%S", temp) %(current_time)) print('%s: Running CoRNN solver.' %(current_time)) for idx in range(num_iters): start_time = time_now.perf_counter() theta_old = theta.copy() # Compute predictions and prediction errors T_data = r_in.shape[0] dhat = np.tanh(x @ theta); n_rec = r_in.shape[1] if solver_type == 'weighted': if u_in is not None: E_pred = (d-dhat) /(1-d**2); n_in = u_in.shape[1] scale = np.sum(abs(E_pred)>threshold,0); else: n_in = 0 if verbose: per_not_used = np.mean(scale) / T_data; # Get the inputs x and targets d if n_in >0: scale = T_data/(T_data-scale); x = np.c_[r_in,u_in]; else: E_pred[abs(E_pred)>threshold] = 0 x = r_in.copy() E_pred = E_pred * scale; d = (r_out - (1-alph)*r_in)/alph elif solver_type == 'standard': d[d<=-1+1e-6] = -1 + 1e-6; E_pred = (d-dhat); d[d>=1-1e-6] = 1-1e-6; per_not_used = 0; elif solver_type == 'robust': if verbose == 2: E_pred = (d-dhat) /(1-d**2); temp = time_now.localtime() per_not_used = 0; current_time = E_pred[abs(E_pred)>threshold] = threshold * time_now.strftime("%H:%M:%S", temp) np.sign(E_pred[abs(E_pred)>threshold]) print('\t %s: Performing initial inverse operations' %(current_time)) # Scale the regularization and compute the Þxed point l2 = l2*T_data; # Perform the Þrst primal variable update theta = Xm @ theta + Xp @ E_pred; conv = np.sqrt(n_rec) * np.sqrt(np.mean((theta [1])); - theta_old ) **2 ))

Ainv = inv(A); if ( (verbose == 2) & (gnd is not None) ) : Xp = Ainv @ x.T; prd = (theta[:n_rec,:].T).ßatten(); APPENDIX A p_cor = pearsonr(gnd,prd)[0]; n_in = u_in.shape[1] else: else: p_cor = np.nan; =

# Get the inputs x and targets d if ((verbose == 2) & (np.mod(idx,300) == 299)): if n_in >0: temp = time_now.localtime() x = np.c_[r_in,u_in]; current_time = else: time_now.strftime("%H:%M:%S", temp) x = r_in.copy() slope = (gnd @ prd) / (gnd@gnd) d = (r_out - (1-alph)*r_in)/alph rmse = np.sqrt(np.mean( (prd-gnd)**2 )) d[d<=-1+1e-6] = -1 + 1e-6; print('\t \t %s: Iteration %d Þnished. RMSE: d[d>=1-1e-6] = 1-1e-6; %.4f. Correlation %.3f. Slope %.3f. Not used %.2f. Convergence %.7f.' \ %(current_time,idx + # Scale the regularization and compute Þxed 1,rmse,p_cor,slope,100*per_not_used,conv)) point l2_scaled = l2*T_data; if ((check_convergence>0) & (idx > 10)): if initialize_fp: if conv < 1e-5: A = x.T @ x + (l2_scaled) * break np.diag(np.ones(x.shape[1])); if verbose: temp = time_now.localtime() Ainv = inv(A); current_time = Xp = Ainv @ x.T; time_now.strftime("%H:%M:%S", temp) z = np.arctanh(d); tt = time_now.perf_counter() - start_time theta_fp = Xp @ z; print('%s: CoRNN solver Þnished in %.2f mins' %(current_time,tt/60)) if initialize_fp: theta = theta_fp return theta else: theta = np.zeros([n_tot,n_rec ]) for idx in range(num_iters): def solve_gradient_descent(r_in,r_out,u_in,alph = 0.1, l2 = 1e-5,verbose = 0, grad = initialize_fp = 0,num_iters = 1e4, compute_gradient(x,d,theta,l2,threshold,solver_ty threshold = 1,learning_rate = 0.001, pe) gnd = None,momentum = del_gd = - grad * learning_rate ; 0,solver_type = 'weighted'): if idx == 0: del_gd_old = del_gd.copy() if verbose: temp = time_now.localtime() delta_gd = momentum * (del_gd_old) + (1- current_time = momentum) * del_gd time_now.strftime("%H:%M:%S", temp) print('%s: Running GD solver.' del_gd_old = delta_gd.copy() %(current_time)) start_time = time_now.perf_counter() theta = theta + delta_gd if ((verbose == 2) & (np.mod(idx,300) == 299) ): T_data = r_in.shape[0] prd = (theta[:n_rec,:].T).ßatten(); n_rec = r_in.shape[1] if gnd is not None: if u_in is not None: p_cor = pearsonr(gnd,prd)[0]; APPENDIX A else: d[d<=-1+1e-6] = -1 + 1e-6; p_cor = np.nan d[d>=1-1e-6] = 1-1e-6; temp = time_now.localtime() current_time = time_now.strftime("%H:%M:%S", temp) # Scale the regularization and compute Þxed rmse = np.sqrt(np.mean( (prd-gnd)**2 )) point slope = (gnd @ prd) / (gnd@gnd) l2_scaled = l2*T_data; print('\t %s: Iteration %d Þnished. RMSE: if initialize_fp: %.4f. Correlation %.3f. Slope %.3f.' \ A = x.T @ x + (l2_scaled) * %(current_time,idx + np.diag(np.ones(x.shape[1])); 1,rmse,p_cor,slope)) Ainv = inv(A); if verbose: Xp = Ainv @ x.T; temp = time_now.localtime() z = np.arctanh(d); current_time = theta_fp = Xp @ z; time_now.strftime("%H:%M:%S", temp) tt = time_now.perf_counter() - start_time num_cores = print('%s: gradient descent Þnished in %.2f np.array(multiprocessing.cpu_count()/2).astype(i mins' %(current_time,tt/60)) nt) if verbose: return theta temp = time_now.localtime() current_time = time_now.strftime("%H:%M:%S", temp) print('%s: NT-descent runs with %d cores.' def )

l2 = 1e-5,verbose = 0, def solve_descent(iters): initialize_fp = 0,num_iters = 10, if initialize_fp: threshold = 1,solver_type = theta_in = theta_fp[:,iters].copy() 'weighted'): else: if verbose: theta_in = None temp = time_now.localtime() beta = current_time = run_nt_algorithm(x,d[:,iters],num_iters,l2,theta_in time_now.strftime("%H:%M:%S", temp) ,threshold,solver_type) print('%s: Running NT solver.' %(current_time)) return beta start_time = time_now.perf_counter() (result) = T_data = r_in.shape[0] Parallel(n_jobs=num_cores)(delayed(solve_desce n_rec = r_in.shape[1] nt)(iters) for iters in range(n_rec)) if u_in is not None: n_in = u_in.shape[1] theta = np.array(result).T else: n_in = 0 if verbose: n_tot = n_rec + n_in; temp = time_now.localtime() theta = np.zeros([n_tot,n_rec]); current_time = time_now.strftime("%H:%M:%S", temp) # Get the inputs x and targets d tt = time_now.perf_counter() - start_time if n_in >0: print('%s: NT-descent Þnished in %.2f mins' x = np.c_[r_in,u_in]; %(current_time,tt/60)) else: x = r_in.copy() return theta d = (r_out - (1-alph)*r_in)/alph APPENDIX A def if solver_type == 'weighted': run_nt_algorithm(x,d_tar,n_iter,l2,theta_in,thresho E_pred = ( (dhat-d_tar) / (1-d_tar**2) ).T; ld,solver_type): scale = np.sum(abs(E_pred)>threshold,0); if np.max(scale) > 0.9 * T_data: if theta_in is not None: scale = beta = theta_in np.sum(abs(E_pred)>10*threshold,0); else: threshold = 10* threshold; beta = np.zeros([x.shape[1]]) scale = T_data/(T_data-scale); for idx in range(n_iter): E_pred[abs(E_pred)>threshold] = 0 E_pred = E_pred * scale; E_pred = E_pred @ x / T_data grad = else: compute_gradient(x,d_tar,beta,l2,threshold,solver E_pred = ( (dhat-d_tar) ).T; _type) E_pred = E_pred @ x / T_data; hes =

gradient = E_pred.T + reg * beta compute_hessian _type) del_nt = -solve(hes, grad,assume_a='pos') return gradient kappa = 1 gamma = 0.4; def chi = 0.6; compute_hessian(x,d_tar,beta,reg,threshold,solve count = 0 r_type): T_data = x.shape[0] while( compute_loss_function(x,d_tar,beta + kappa * del_nt,l2) >= \ dhat = np.tanh( x @ beta ) compute_loss_function(x,d_tar,beta,l2) + gamma * kappa* (grad @ del_nt) ): if solver_type == 'weighted': E_pred = ( (dhat-d_tar) / (1-d_tar**2) ); kappa = chi * kappa scale = np.sum(abs(E_pred)>threshold,0); count = count + 1 if np.max(scale) > 0.9 * T_data: if count == 4: scale = kappa = 1; np.sum(abs(E_pred)>10*threshold,0); break threshold = 10* threshold; scale = T_data/(T_data-scale); delta_gd = kappa * del_nt temp = np.sqrt(1-dhat**2) / np.sqrt(1- beta = beta + delta_gd d_tar**2) temp[abs(E_pred)>threshold] = 0 temp = temp * scale; else: return beta temp = np.sqrt(1-dhat**2); def r = x * temp[:,None]; compute_gradient(x,d_tar,beta,reg,threshold,solv er_type): T_data = x.shape[0] H = r.T @ r /r.shape[0] + reg * np.diag(np.ones(r.shape[1])) return H dhat = np.tanh( x @ beta ) APPENDIX A class Model(nn.Module): def compute_loss_function(x,d,beta,reg): def __init__(self, n_rec, n_in): super(Model, self).__init__() dhat = np.tanh( x @ beta ) self.n_rec = n_rec dhat[dhat>=1] = 1-1e-16 self.n_in = n_in dhat[dhat<=-1] = -1+1e-16 self.linear = nn.Linear(n_rec + n_in, loss = - np.mean( 0.5*(d+1)* n_rec,bias=False) np.log(0.5*(dhat+1)) + 0.5*(1-d)* np.log(0.5*(1- self.tanh = nn.Tanh() dhat)) ) + reg * np.sum(beta**2)/2 def forward(self, x_in): return loss #x_in = torch.tensor(x_in, dtype=torch.ßoat32) d_out = self.tanh(self.linear(x_in)) def solve_pytorch(r_in,r_out,u_in,alph = 0.1, return d_out l2 = 1e-5,verbose = 0, initialize_fp = 0,num_iters = 1e4, d = torch.tensor(d, dtype=torch.ßoat32) threshold = 2,learning_rate = 0.001, x = torch.tensor(x, dtype=torch.ßoat32) gnd = None,solver_type = 'logistic'): model = Model(n_rec, n_in) if verbose: # deÞne the multi target logistic regression loss temp = time_now.localtime() function current_time = if solver_type == 'logistic': time_now.strftime("%H:%M:%S", temp) criterion = nn.BCELoss() print('%s: Running Pytorch-GD solver.' else: %(current_time)) criterion = nn.MSELoss() start_time = time_now.perf_counter() optimizer = optim.Adam(model.parameters(), lr=learning_rate, weight_decay=l2) n_rec = r_in.shape[1] if initialize_fp: if u_in is not None: model.linear.weight.data = n_in = u_in.shape[1] torch.tensor(theta_fp.T, dtype=torch.ßoat32) else: for idx in range(num_iters): n_in = 0 d_out = model(x) # Get the inputs x and targets d # compute the loss if n_in >0: loss = criterion((1 + d_out)/2, (1 + d)/2) x = np.c_[r_in,u_in]; else: # backward x = r_in.copy() loss.backward() d = (r_out - (1-alph)*r_in)/alph d[d<=-1+1e-6] = -1 + 1e-6; d[d>=1-1e-6] = 1-1e-6; # update the weights if initialize_fp: optimizer.step() T_data = x.shape[0] # zero the gradients A = x.T @ x + l2 *T_data *1e4* optimizer.zero_grad() np.diag(np.ones(x.shape[1])); if ((verbose == 2) & (np.mod(idx,300) == 299) ): Ainv = inv(A); weights = model.linear.weight Xp = Ainv @ x.T; prd = weights[:, z = np.arctanh(d); :n_rec].detach().numpy().ßatten() theta_fp = Xp @ z; if gnd is not None: p_cor = pearsonr(gnd,prd)[0]; else: # DeÞne the network p_cor = np.nan APPENDIX A slope = (gnd @ prd) / (gnd@gnd) x = r_in.copy() d = (r_out - (1-alph)*r_in)/alph temp = time_now.localtime() d[d<=-1+1e-6] = -1 + 1e-6; current_time = d[d>=1-1e-6] = 1-1e-6; time_now.strftime("%H:%M:%S", temp) rmse = np.sqrt(np.mean( (prd-gnd)**2 )) # DeÞne the network print('\t %s: Iteration %d Þnished. RMSE: class Model(nn.Module): %.4f Correlation %.3f. Slope %.3f.' \ def __init__(self, n_rec, n_in): %(current_time,idx + super(Model, self).__init__() 1,rmse,p_cor,slope)) self.n_rec = n_rec self.n_in = n_in self.linear = nn.Linear(n_rec + n_in, if verbose: n_rec,bias=False) temp = time_now.localtime() self.tanh = nn.Tanh() current_time = time_now.strftime("%H:%M:%S", temp) def forward(self, x_in): tt = time_now.perf_counter() - start_time #x_in = torch.tensor(x_in, print('%s: Pytorch Þnished in %.2f mins' dtype=torch.ßoat32) %(current_time,tt/60)) d_out = self.tanh(self.linear(x_in)) return d_out return model.linear.weight.detach().numpy().T d = torch.tensor(d, dtype=torch.ßoat32,device = def solve_pytorch_gpu(r_in,r_out,u_in,alph = 0.1, gpu_device) l2 = 1e-5,verbose = 0, x = torch.tensor(x, dtype=torch.ßoat32,device = initialize_fp = 0,num_iters = 1e4, gpu_device) threshold = 2,learning_rate = 0.001, gnd = None,solver_type = 'logistic'): if initialize_fp: if verbose: T_data = x.shape[0] temp = time_now.localtime() n_tot = n_rec+ n_in current_time = reg_term = torch.tensor((l2) * T_data * 1e4 * time_now.strftime("%H:%M:%S", temp) np.diag(np.ones(n_tot)),\ print('%s: Running Pytorch-GD solver.' device=gpu_device,dtype = %(current_time)) torch.ßoat32) start_time = time_now.perf_counter() A = x.T @ x + reg_term; n_rec = r_in.shape[1] if u_in is not None: try: n_in = u_in.shape[1] Ainv = torch.linalg.inv(A); else: except: n_in = 0 Ainv = torch.linalg.inv(A.cpu()).to(gpu_device); try: Xp = Ainv @ x.T; gpu_device = torch.device("cuda:0"); z = torch.arctanh(d); a = torch.tensor(1,device = gpu_device); theta_fp = Xp @ z; except: theta_fp = theta_fp.cpu().detach().numpy() gpu_device = torch.device("mps"); model = Model(n_rec, n_in) # Get the inputs x and targets d model.to(gpu_device) if n_in >0: # deÞne the multi target logistic regression loss x = np.c_[r_in,u_in]; function else: if solver_type == 'logistic': APPENDIX A criterion = nn.BCELoss() else: return criterion = nn.MSELoss() model.linear.weight.cpu().detach().numpy().T optimizer = optim.Adam(model.parameters(), lr=learning_rate, weight_decay=l2) if initialize_fp: def Þt_FORCE(r,u=None,alph = 0.1, model.linear.weight.data = lam = 100,lam_in = 1,g_in = 3,verbose = 0, torch.tensor(theta_fp.T, initialize_fp = 0,num_iters = 1e3, dtype=torch.ßoat32).to(gpu_device) gnd = None,solver_type = 'currents'): for idx in range(num_iters): if verbose: d_out = model(x) temp = time_now.localtime() current_time = # compute the loss time_now.strftime("%H:%M:%S", temp) loss = criterion((1 + d_out)/2, (1 + d)/2) print('%s: Running FORCE solver on %s.' %(current_time,solver_type)) # backward loss.backward() def

= None): if x = np.dot(w_rec,r_before) + # update the weights np.dot(w_in,u_in); optimizer.step() else: x = np.dot(w_rec,r_before); # zero the gradients rout = (1-alph) * r_before + alph * np.tanh(x) optimizer.zero_grad() return rout,x if ((verbose == 2) & (np.mod(idx,300) == 299)): if initialize_fp: weights = model.linear.weight if u is not None: prd = weights[:, n_in = u.shape[1] :n_rec].cpu().detach().numpy().ßatten() else: if gnd is not None: n_in = 0 p_cor = pearsonr(gnd,prd)[0]; r_in = r[:-1,:] else: r_out = r[1:,:] p_cor = np.nan d = (r_out - (1-alph)*r_in)/alph slope = (gnd @ prd) / (gnd@gnd) d[d<=-1+1e-6] = -1 + 1e-6; d[d>=1-1e-6] = 1-1e-6; temp = time_now.localtime() if n_in >0: current_time = x = np.c_[r_in,u]; time_now.strftime("%H:%M:%S", temp) else: rmse = np.sqrt(np.mean( (prd-gnd)**2 )) x = r_in.copy() print('\t %s: Iteration %d Þnished. RMSE: T_data = x.shape[0] %.4f Correlation %.3f. Slope %.3f.' \ A = x.T @ x + lam/T_data * %(current_time,idx + np.diag(np.ones(x.shape[1])); 1,rmse,p_cor,slope)) Ainv = inv(A); if verbose: Xp = Ainv @ x.T; temp = z = np.arctanh(d);

theta_fp = Xp @ z; tt = time_now.perf_counter() - start_time print('%s: Pytorch Þnished in %.2f mins' T = r.shape[0]-1 %(current_time,tt/60)) n_rec= r.shape[1] APPENDIX A if u is not None and u.shape[0] != T: elif solver_type == 'Þring_rates': raise Exception('There is a missmatch e_min = (r_out - r[i+1,:])/alph; between u and r dimensions!!!') Pxr = P @ r_bef; rxPxr = r_bef @ (P @ r_bef) + 1 delP = - np.outer(Pxr, Pxr) / rxPxr P = P + delP; ,n_ delW = - np.outer(e_min, P @ r_bef) rec*n_rec)

np.Þll_diagonal(temp,0) w_rec = temp w_rec = w_rec + delW; if u is not None: if n_in > 0: n_in = u.shape[1] Pxu = P_in @ u[i] w_in uxPxu = u[i] @ (P_in @ u[i]) + 1 =np.random.normal(0,g_in/np.sqrt(n_in),n_in*n_r delP = - np.outer(Pxu, Pxu) / uxPxu ec).reshape(n_rec,n_in) P_in = P_in + delP; else: delW = - np.outer(e_min, P_in @ u[i]) n_in = 0 w_in = w_in + delW; r_bef,x_pred = predict_single(r_bef,u[i]) if initialize_fp: else: w_rec = theta_fp[:n_rec,:].T r_bef,x_pred = predict_single(r_bef) if n_in >0: w_in = if solver_type == 'currents': theta_fp[n_rec:,:].T.reshape(n_rec,n_in) e_pls = x_pred - x_now if (verbose == 2): elif solver_type == 'Þring_rates': prd = w_rec.ßatten() e_pls = (r_bef - r[i+1,:])/alph; if gnd is not None: p_cor = pearsonr(gnd,prd)[0]; if ((verbose == 2)): else: prd = w_rec.ßatten() p_cor = np.nan if gnd is not None: print('\t %s: FP initialized: Correlation %.3f.' \ p_cor = pearsonr(gnd,prd)[0]; %(current_time,p_cor)) else: p_cor = np.nan slope = (gnd @ prd) / (gnd@gnd) for k in range(num_iters): P = np.diag(np.ones(n_rec))/lam; temp = time_now.localtime() if n_in > 0: current_time = P_in = np.diag(np.ones(n_in))/lam_in; time_now.strftime("%H:%M:%S", temp) r_bef = r[0,:]; rmse = np.sqrt(np.mean( (prd-gnd)**2 )) for i in range(T): print('\t %s: Iteration %d Þnished. RMSE: temp = (r[i+1,:] - (1-alph)*r[i,:])/alph %.4f Correlation %.3f. Slope %.3f. Conv. %.3f.' \ temp[temp>=1-1e-6] = 1-1e-6 %(current_time,k + temp[temp<=-1+1e-6] = -1+1e-6 1,rmse,p_cor,slope,np.min(e_pls/e_min))) x_now = np.arctanh(temp) if if n_in > 0: theta r_out,x_out = predict_single(r_bef,u[i]) else:

else: theta = w_rec r_out,x_out = predict_single(r_bef) return theta.T if solver_type == 'currents': e_min = x_out - x_now; APPENDIX A def solve_corrn_gpu(r_in,r_out,u_in = None, alph = # Scale the regularization and compute the Þxed 0.1, ßoat_type = '32bit', point l2 = 1e-4, threshold = 0.2,initialize_fp = 1, l2 = l2*T_data; verbose = 0,mask = None, check_convergence = 0, num_iters = 30,gnd = None,solver_type = 'weighted'): if mask is None: if verbose: mask = temp = time_now.localtime() np.r_[np.diag(np.ones(n_rec)),np.zeros([n_in,n_re current_time = c]) ] time_now.strftime("%H:%M:%S", temp) mask = (mask == 1) solver.' % reg_term = torch.tensor((l2) *

() np.diag(np.ones(n_tot)),\ try: device=gpu_device, dtype=dtype) gpu_device = torch.device("cuda:0"); a = torch.tensor(1,device = gpu_device); A = x.T @ x + reg_term; except: gpu_device = torch.device("mps"); num_iters = int(num_iters); try: T_data = r_in.shape[0] Ainv = torch.linalg.inv(A); n_rec = r_in.shape[1] except: if u_in is not None: Ainv = torch.linalg.inv(A.cpu()).to(gpu_device); n_in = u_in.shape[1] Xp = Ainv @ x.T; else: Xm = Xp @ x; n_in = 0 z = torch.arctanh(d); theta_fp = Xp @ z; # Get the inputs x and targets d if initialize_fp: if n_in >0: theta = theta_fp; x = np.c_[r_in,u_in]; else: else: w_rec = x = r_in.copy() np.random.normal(0,1/np.sqrt(n_rec),n_rec*n_tot d = (r_out - (1-alph)*r_in)/alph ).reshape(n_tot,n_rec) d[d<=-1+1e-6] = -1 + 1e-6; theta = d[d>=1-1e-6] = 1-1e-6; torch.tensor(w_rec,device=gpu_device, if ßoat_type == '32bit': dtype=dtype); dtype = torch.ßoat32 if ( (verbose == 2) & (gnd is not None) ) : else: prd = (theta[:n_rec,:].T).ßatten().cpu().numpy(); dtype = torch.ßoat64 p_cor = pearsonr(gnd,prd)[0]; n_tot = x.shape[1] print('\t \t Fixed point correlation: %.3f' %p_cor) x = torch.tensor(x,device=gpu_device, else: dtype=dtype) p_cor = np.nan; d = torch.tensor(d,device=gpu_device, dtype=dtype) if verbose == 2: temp = time_now.localtime() if verbose == 2: current_time = temp =

time_now.strftime("%H:%M:%S", temp) current_time = print('\t %s: Starting iterations' time_now.strftime("%H:%M:%S", temp) %(current_time)) print('\t %s: Performing initial inverse operations' %(current_time)) for idx in range(num_iters): # Compute predictions and prediction errors APPENDIX A dhat = torch.tanh(x @ theta); print('%s: CoRNN solver Þnished in %.2f mins' if solver_type == 'weighted': %(current_time,tt/60)) E_pred = (d-dhat) /(1-d**2); per_not_used = return theta.cpu().numpy() torch.sum(abs(E_pred)>threshold) / (n_rec * T_data) def Þt_FORCE_gpu(r,u=None,alph = 0.1, E_pred[abs(E_pred)>threshold] = 0 lam = 100,lam_in = 1,g_in = 3,verbose = 0, elif solver_type == 'standard': initialize_fp = 0,num_iters = 1e3, E_pred = (d-dhat); gnd = None,solver_type = 'currents'): per_not_used = 0; if verbose: elif solver_type == 'robust': temp = time_now.localtime() E_pred = (d-dhat) /(1-d**2); current_time = per_not_used = 0; time_now.strftime("%H:%M:%S", temp)

E_pred[abs(E_pred)>threshold] = threshold * print('%s: Running FORCE solver on np.sign(E_pred[abs(E_pred)>threshold]) %(current_time,solver_type)) try: gpu_device = torch.device("cuda:0"); a = torch.tensor(1,device = gpu_device); except: gpu_device = torch.device("mps"); # Perform the Þrst primal variable update def predict_single(r_before,u_in = None): theta = Xm @ theta + Xp @ E_pred if u_in is not None: x = w_rec @ r_before +w_in @u_in; else: if ( (verbose == 2) & (gnd is not None) ) : x = w_rec @r_before; prd = rout = (1-alph) * r_before + alph * torch.tanh(x) (theta[:n_rec,:].T).ßatten().cpu().numpy(); p_cor = pearsonr(gnd,prd)[0]; return rout,x else: p_cor = np.nan; if initialize_fp: if u is not None: n_in = u.shape[1] else: if ((verbose == 2) & (np.mod(idx,300) == 299) ): n_in = 0 temp = time_now.localtime() r_in = r[:-1,:] current_time = r_out = r[1:,:] time_now.strftime("%H:%M:%S", temp) d = (r_out - (1-alph)*r_in)/alph slope = (gnd @ prd) / (gnd@gnd) d[d<=-1+1e-6] = -1 + 1e-6; rmse = np.sqrt(np.mean( (prd-gnd)**2 )) d[d>=1-1e-6] = 1-1e-6; print('\t \t %s: Iteration %d Þnished. RMSE: if n_in >0: %.4f. Correlation %.3f. Slope %.3f. Not used %.2f.' x = np.c_[r_in,u]; \ else: %(current_time,idx + x = r_in.copy() 1,rmse,p_cor,slope,100*per_not_used)) x = torch.tensor(x,device=gpu_device,dtype = torch.ßoat32) if verbose: d = torch.tensor(d,device=gpu_device,dtype = temp = time_now.localtime() torch.ßoat32) current_time = time_now.strftime("%H:%M:%S", temp)

tt = time_now.perf_counter() - start_time APPENDIX A reg_term = torch.tensor(lam/T_data * w_in = np.diag(np.ones(x.shape[1])),\ torch.tensor(theta_fp[n_rec:,:].clone().detach().T.r device=gpu_device,dtype = eshape(n_rec,n_in),device=gpu_device,dtype = torch.ßoat32) torch.ßoat32) if (verbose == 2): A = x.T @ x + reg_term; prd = w_rec.ßatten().cpu().numpy(); if gnd is not None: p_cor = pearsonr(gnd,prd)[0]; try: else: Ainv = torch.linalg.inv(A); p_cor = np.nan except: print('\t %s: FP initialized: Correlation %.3f.' \ Ainv = %(current_time,p_cor)) torch.linalg.inv(A.cpu()).to(gpu_device); Xp = Ainv @ x.T; z = torch.arctanh(d); for k in range(num_iters): theta_fp = Xp @ z; P = torch.tensor(np.diag(np.ones(n_rec))/lam,device= gpu_device,dtype = torch.ßoat32); r = torch.tensor(r,device=gpu_device,dtype = if n_in > 0: torch.ßoat32) P_in = alph = torch.tensor(np.diag(np.ones(n_in))/lam_in,device torch.tensor(alph,device=gpu_device,dtype = =gpu_device,dtype = torch.ßoat32); torch.ßoat32) r_bef = r[0,:]; T = r.shape[0]-1 for i in range(T): n_rec= r.shape[1] temp = (r[i+1,:] - (1-alph)*r[i,:])/alph if u is not None and u.shape[0] != T: temp[temp>=1-1e-6] = 1-1e-6 raise Exception('There is a missmatch temp[temp<=-1+1e-6] = -1+1e-6 between u and r dimensions!!!') x_now = torch.arctanh(temp) if n_in > 0: r_out,x_out = predict_single(r_bef,u[i]) else: temp=np.random.normal(0,g_in/np.sqrt(n_rec),n_ r_out,x_out = predict_single(r_bef) rec*n_rec).reshape(n_rec,n_rec) np.Þll_diagonal(temp,0) if solver_type == 'currents': w_rec = e_min = x_out - x_now; torch.tensor(temp,device=gpu_device,dtype = elif solver_type == 'Þring_rates': torch.ßoat32) e_min = (r_out - r[i+1,:])/alph; if u is not None: Pxr = P @ r_bef; n_in = u.shape[1] rxPxr = r_bef @ (P @ r_bef) + 1 w_in delP = - torch.outer(Pxr, Pxr) / rxPxr =np.random.normal(0,g_in/np.sqrt(n_in),n_in*n_r P = P + delP; ec).reshape(n_rec,n_in) delW = - torch.outer(e_min, P @ r_bef) w_in = torch.tensor(w_in,device=gpu_device,dtype = torch.ßoat32) w_rec = w_rec + delW; else: n_in = 0 if n_in > 0: Pxu = P_in @ u[i] if initialize_fp: uxPxu = u[i] @ (P_in @ u[i]) + 1 w_rec = theta_fp[:n_rec,:].clone().detach().T delP = - torch.outer(Pxu, Pxu) / uxPxu if n_in >0: P_in = P_in + delP; APPENDIX A delW = - torch.outer(e_min, P_in @ u[i]) f = np.load('experiment1a_results_fp.npz') w_in = w_in + delW; slopes = f['slopes'] r_bef,x_pred = predict_single(r_bef,u[i]) rmse = f['rmse'] else: rmedse = f['rmedse'] r_bef,x_pred = predict_single(r_bef) cor_p = f['cor_p'] cor_s = f['cor_s'] if solver_type == 'currents': times = f['times'] e_pls = x_pred - x_now elif solver_type == 'Þring_rates': # 0 is CoRNN, CPU e_pls = (r_bef - r[i+1,:])/alph; # 1 is Pytorch, logistic, GPU # 2 is Pytorch, logistic, CPU if ((verbose == 2)): # 3 is NT, weighted prd = w_rec.ßatten().cpu().numpy(); # 4 is Pytorch, l2, GPU if gnd is not None: # 5 is NT, logistic p_cor = pearsonr(gnd,prd)[0]; # 6 is Force, currents else: # 7 is Force, Þring rates p_cor = np.nan # 8 is CoRNN, GPU slope = (gnd @ prd) / (gnd@gnd) temp = time_now.localtime() current_time = labels_all = ["Fixed point, weighted, time_now.strftime("%H:%M:%S", temp) CPU","Pytorch, CE, GPU", \ (np.mean( (prd-gnd)**2 )) "Pytorch, CE, CPU", "Newton, weighted", %d Þnished. RMSE: "Pytorch, l2, GPU", \

Slope %.3f. Conv. %.3f.' \ "Newton, CE", "Force, currents", "Force, %(current_time,k + Þring rates", \ 1,rmse,p_cor,slope,torch.min(e_pls/e_min))) "Fixed point, weighted, GPU"] w_rec = w_rec.cpu().numpy() for pick in [8,0,3,5,1,2,4,6,7]: if n_in>0: temp = cor_p[pick,:,:].mean(0) w_in = w_in.cpu().numpy() temp_std = np.std(cor_p[pick,:,:],0) / theta = np.c_[w_rec,w_in] np.sqrt(cor_p.shape[1]) else: x = np.nanmean(times[pick,:,:],0) theta = w_rec x_std = np.std(times[pick,:,:],0) / np.sqrt(cor_p.shape[1]) return theta.T plt.errorbar(temp,x,x_std,temp_std,label = labels_all[pick]) Figure scripts - plot_subÞgure.py #!/usr/bin/env python3 plt.yscale('log') # -*- coding: utf-8 -*- """ plt.xlim([0.84,1.01]) Created on Thu Jan 1913:56:082023 plt.xticks([0.9,1],fontsize=14) plt.yticks([0.1,1,10,100],fontsize=14) @author: dinc #plt.xlabel('Correlation between ground truth and """ inferred weights',fontsize = 15) plt.ylabel('Time (s)',fontsize = 15) import numpy as np import matplotlib.pyplot as plt plt.subplot(3,1,1) plt.subplot(3,1,2)

Claims

APPENDIX A f = np.load('experiment1b_results_fp.npz') rmse = f['rmse'] slopes = f['slopes'] rmedse = f['rmedse'] rmse = f['rmse'] cor_p = f['cor_p'] rmedse = f['rmedse'] cor_s = f['cor_s'] cor_p = f['cor_p'] times = f['times'] cor_s = f['cor_s'] times = f['times'] # 0 is CoRNN, CPU # 1 is Pytorch, logistic, GPU # 0 is CoRNN, CPU # 2 is Pytorch, logistic, CPU # 1 is Pytorch, logistic, GPU # 3 is NT, weighted # 2 is Pytorch, logistic, CPU # 4 is Pytorch, l2, GPU # 3 is NT, weighted # 5 is NT, logistic # 4 is Pytorch, l2, GPU # 6 is Force, currents # 5 is NT, logistic # 7 is Force, Þring rates # 6 is Force, currents # 8 is CoRNN, GPU # 7 is Force, Þring rates # 8 is CoRNN, GPU labels_all = ["Fixed point, weighted, CPU","Pytorch, CE, GPU", \ labels_all = ["Fixed point, weighted, "Pytorch, CE, CPU", "Newton, weighted", CPU","Pytorch, CE, GPU", \ "Pytorch, l2, GPU", \ "Pytorch, CE, CPU", "Newton, weighted", "Newton, CE", "Force, currents", "Force, "Pytorch, l2, GPU", \ Þring rates", \ "Newton, CE", "Force, currents", "Force, "Fixed point, weighted, GPU"] Þring rates", \ "Fixed point, weighted, GPU"] for pick in [8,0,3,5,1,2,4,6,7]: for pick in [8,0,3,5,1,2,4,6,7]: temp = cor_p[pick,:,:].mean(0) temp_std = np.std(cor_p[pick,:,:],0) / temp = cor_p[pick,:,:].mean(0) np.sqrt(cor_p.shape[1]) temp_std = np.std(cor_p[pick,:,:],0) / x = np.nanmean(times[pick,:,:],0) np.sqrt(cor_p.shape[1]) x_std = np.std(times[pick,:,:],0) / x = np.nanmean(times[pick,:,:],0) np.sqrt(cor_p.shape[1]) x_std = np.std(times[pick,:,:],0) / plt.errorbar(temp,x,x_std,temp_std,label = np.sqrt(cor_p.shape[1]) labels_all[pick]) plt.errorbar(temp,x,x_std,temp_std,label = labels_all[pick]) plt.xlabel('Correlation between ground truth and plt.xlim([0.84,1.01]) inferred weights',fontsize = 15) plt.xticks([0.9,1]) plt.ylabel('Time (s)',fontsize = 15) plt.yticks([0.1,1,10,100],fontsize=14) plt.yscale('log') plt.ylabel('Time (s)',fontsize = 15) plt.yscale('log') plt.legend(loc='center left', bbox_to_anchor=(1, 0.5)) plt.subplot(3,1,3) f = np.load('experiment1c_results_fp.npz') slopes = f['slopes'] Attorney Docket No.: NTTR-2023003pct 14. The method of claim 13, wherein the constraints include preventing self-excitation of neurons within the neural network. 15. The method of any one of claims 2-4 and 11-14, wherein the convex solver is configured to utilize biological priors for inductive biases during training. 16. The method of any one of claims 2-4 and 11-15, wherein the convex solver is configured to handle input noise and conversion noise during the training of the neural network. 17. The method of any one of claims 2-4 and 11-16, wherein the convex solver is configured to perform outlier detection during the training process. 18. The method of any one of claims 2-4 and 11-17, wherein the convex solver is configured to scale to datasets collected from more than 10,000 neurons for offline processing. 19. A non-transitory computer-readable medium comprising instructions that, when executed by one or more processors, cause a computer system to perform the method of any one of claims 1- 18. 20. A computer system for training a neural network to interpret neural dynamics, the system comprising: one or more processors; and a non-transitory computer-readable medium including instructions that, when executed by the one or more processors, cause the computer system to perform the method of any one of claims 1-18. 21. The computer system of claim 22, wherein the neural network is a data constrained recurrent neural network. 88