Continuous multinomial logistic regression for neural decoding
Authors: Anuththara Rupasinghe
Jonathan W. Pillow
Venue: ICLR 2026 Poster (2026)
Keywords: Neural propulation coding
Conditional Density Estimation
Gaussian processes
Variational inference
Probabilistic models
Primary area: applications to neuroscience & cognitive science
TL;DR: We introduce Continuous Multinomial Logistic Regression (CMLR), a Gaussian-process-regularized exponential-family model for scalable and flexible conditional density estimation in continuous neural decoding tasks.

Abstract: Multinomial logistic regression (MLR) is a classic model for multi-class classification that has been widely used for neural decoding. However, MLR requires a finite set of discrete output classes, limiting its applicability to settings with continuous-valued outputs (e.g., time, orientation, velocity, or spatial position), which are common in neuroscience settings. To address this limitation, we propose Continuous Multinomial Logistic Regression (CMLR), a generalization of logistic regression to continuous output spaces. CMLR represents a novel exponential-family model for conditional density estimation (CDE), mapping neural population activity to a full probability density over external covariates. It captures the influence of each neuron’s activity on the decoded variable through a smooth, interpretable tuning function, regularized by a Gaussian process prior. The resulting nonparametric decoding model flexibly captures asymmetric and multimodal densities, and accommodates both linear and circular variables. To illustrate the performance of CMLR, we applied it to large-scale datasets from mouse and monkey visual cortex, mouse hippocampus, and monkey motor cortex, where it generally outperformed a wide variety of other decoding methods, including deep neural networks (DNNs), XGBoost, and FlexCode. It also outperformed a closely-related correlation-blind decoder, highlighting the importance of correlations for accurate neural decoding. The CMLR model provides a scalable, flexible, and interpretable method for decoding continuous variables from diverse brain regions.

Published as a conference paper at ICLR 2026
CONTINUOUS MULTINOMIAL LOGISTIC REGRESSION FOR NEURAL DECODING
Anuththara Rupasinghe & Jonathan W. Pillow Princeton Neuroscience Institute, Princeton University (a (ar0621, p 11ow) @princetoon.co edu
ABSTRACT
Multinomial logistic regression (MLR) is a classic model for multi-class classifica- tion that has been widely used for neural decoding. However, MLR requires a finite set of discrete output classes, limiting its applicability to settings with continuous- valued outputs (e.g., time, orientation, velocity, or spatial position), which are com- mon in neuroscience settings. To address this limitation, we propose Continuous Multinomial Logistic Regression (CMLR), a generalization of logistic regression continuous output spaces. CMLR represents a novel exponential-famme for conditional density estimation (CDE), mapping neural population activity to a full probability density neuron's 's activity on the decoded variable through a smooth, interpretable tuning function, regularized by Gaussian process prior. The resulting nonparametric decoding model flexibly captures asymmetric and multimodal densities, and ac- commodates both linear and circular variables. To illustrate the performance of CMLR, we applied it to large-scale datasets from mouse and monkey visual cortex, mouse hippocampus, and monkey motor cortex, where it generally outperformed a wide variety of other decoding methods, including deep neural networks (DNNs), XGBoost, and FlexCode. It also outperformed closely-related correlation-bling decoder, highlighting the importance of correlations for accurate neural decod- ing. The CMLR model provides a scalable, flexible, and interpretable method for decoding continuous variables from diverse brain regions.
model
external covariates. It captures the influence of each
INTRODUCTION
Neural decoding refers to the problem of estimating behavioral or sensory variables from neural activity, a central challenge in neuroscience (Georgopoulos et al. 1982] Brown al. 1998). Logistic regression is foundational model for binary classification that has been widely applied to neural decoding problems with two alternatives (Ryali et al. 2010; Glaser et al., 2020). For tasks with multiple alternatives, multinomial logistic regression (MLR) provides natural extension, defining class probabilities based on the linear projection of neural activity onto weight vectors, one for each output class (Huttunen al., 2013 Song et al.] 2014 Greenidge et al., 2024). However, many neural decoding tasks involve continuous variables such as time, orientation, head direction, position, or velocity (Paninski & Cunningham] 2018). In such settings, standard regression models are frequently inadequate because they yield only point predictions and cannot represent multimodal, circular, or asymmetric output distributions (Hyndman et al.] 1996). To obtain full predictive distributions, researchers often adapt MLR- classifiers by discretizing the continuous output into bins (Stringer al.) 202T] Greenidge et al., (024). However, this discretization reduces effective resolution, introduces quantization artifacts, and typically necessitates additional regularization to prevent overfitting (Altman & Royston
2006 Nojavan A. et al. 2017).
To overcome this limitation, we introduce the Continuous Multinomial Logistic Regression (CMLR) model, which generalizes MLR to continuous output spaces. Whereas MLR defines a discrete probability distribution over finite set of output classes using log-linear combination of weight vectors, the CMLR model defines a probability density function (pdf) over a continuous output variable using log-linear combination of weight functions. These weight functions conceptually analogous to neural tuning curves, as they characterize how each neuron's 's activity influences the predicted density over the variable of interest and often resemble the underlying tuning, although they

---

Published as a conference paper at ICLR 2026
input vector
weight functions
output density
w(y)
ply x) x exp(w(u) x)
x
21
wi (v)
12
exp
X
:
DD
ymin y
ymax
Umin
mmax
oulput variable
oulput
Figure I: Continuous Multinomial Logistic Regression (CMLR) model schematic. Given an input feature vector RDD, the model defines logistic density over y using weight functions w(y) {wa(v))amr The probability density at output value y given x is proportional to exp(w(y) T x). For each feature d, the function wa(y) has independent Gaussian process prior to induce smoothness.
originate from decoding model rather than encoding model. CMLR thus provides a framework for conditional density estimation (CDE), a problem setting that involves mapping a feature vector to probability density (Hyndman et al. (996). For neural decoding applications, this density can be summarized by its mean or mode to obtain point estimates of the external variable of interest.
To enforce smoothness, we place Gaussian process (GP) priors on the weight functions, approach closely related to Logistic Gaussian Process density estimation (Tokdar et al.] 2004 Rithimäki & Vehtari 2014). Given vector of neural activities, the resulting model can represent complex conditional densities over the output variable, including multimodal, asymmetric, and circular forms often observed in neural decoding. To make inference tractable, we develop a Fourier-domain stochastic variational inference algorithm (Hoffman et 20133,, leveraging GP kernel stationarity for diagonalization and spectral truncation (Hensman et al. (018). The resulting approach handles high- dimensional inputs, extends naturally to multivariate-output CDE via anisotropic kernels (Rasmussen Williams 2006), and handles large-scale neural datasets.
We evaluate CMLR on large-scale datasets from mouse and monkey primary visual cortex (V1) (Stringer et al. 2021; Graf et al., (01I). mouse hippocampus CAI (Hazon et al., 2022 Jercog & Schnitzer 2022). and monkey motor cortex (Glaser et al. 2018 2020). Across datasets, CMLR generally outperforms strong decoding baselines, including Extreme Gradient Boosting (XGBoost) (Chen & Guestrin 2016), deep neural networks (DNNs) (Orbcc (962). and the leading non- parametric CDE method Clex (Izbicki & Leel 2017), while producing well-calibrated posterior densities. These advantages are most pronounced in low-data regimes, where CMLR's GP-based functional priors and additive structure provide strong regularization and limit overfitting, and in tasks with structured or multimodal outputs such as circular orientation decoding, where accurate representation of periodicity and multimodality is essential. In these settings, access to the full conditional density allows CMLR to capture the underlying structure far more faithfully than point- estimate models such as XGBoost and DNNs, which cannot naturally model circular or multimodal output distributions. Furthermore, extending prior comparisons of correlation-awan and correlation- blind decoders (Nirenberg & Latham 2003 Graf et al. [011 Greenidge et al. (024), we show that CMLR (correlation- aware), by modeling shared variability, surpasses Naive Bayes (correlation- blind) in predictive accuracy across V1, CA1, and motor cortex. These results establish CMLR as a flexible, interpretable, and scalable framework for high-resolution neural decoding across brain regions.
CONTINUOUS MULTINOMIAL LOGISTIC REGRESSION (CMLR) MODEL
The CMLR model (schematized in Fig. defines a mapping from input vector x E IR (e.g.. vector of spike counts) to probability density over an output variable y E Ω defined over some compact domain Q. The model parameters consist of a set of D weight functions, one for each element of x, where wa(y)) describes the additive influence of input feature xa on the log-density over y. The resulting conditional density is given by:
exp
W (1)
x
P(Y = x =
(I)
Jn exp (w(y') dy""
2

---

Published as a conference paper at ICLR 2026
where w(y) = (y), wyyy) denotes the vector obtained by evaluating all weight functions at y. and the denominator provides the normalizing constant, ensuring the density integrates to 1 (Tokdar et al. 2004; Riihimäii & Vehtari 2014).
Note that this generalizes standard multinomial logistic regression (MLR) with K discrete classes, where the probability of class k is given by:
exp(w x)
p(Y = k xx)
(2)
Lj=1
exp(w]
where we denotes the weight vector for class CMLR can thus be viewed as the continuous limit of MLR as K + so, replacing a set of ordered discrete class weights with smooth weight functions defined over a continuous output space.
To enforce smoothness in the decoding weights w(y), we place an independent Gaussian Process (GP) prior on each weight function:
wa(y) ~ GP (0, Ka);
where Ka the covariance function governing the weight function: cov(wa(s"), wa(y"))
Kaly', ("") for any pair of output values y'.y" & Q2. Here, we employ the standard radial basis func-
tion (RBF) covariance (Rasmussen & Williams, 2006): Ka (y', ,y") Pa exp (- - )²
(- 222
where Pa and la are hyperparameters controlling the marginal variance (i.e., amplitude) and length scale (i.e., smoothness) of the weight function, respectively. For circular output spaces (e.g., orien- tations with Ω [0, 2T)). we use the periodic version of the RBF kernel defined as (Schölkopf
Smola 2001): y Ka (1/,"")) v/"") Pa Em=-oo exp (-(1' - yy " + 2mm) /(2()) around the unit circle.
ensuring continuity
3 EFFICIENT LEARNING VIA STOCHASTIC VARIATIONAL. INFERENCE (SVI)
Given dataset of N input-output pairs D = {(x;: Un)]n= aim to jointly fit the weight functions w(y) {wa(t)) and the hyperparameters = {ps, fa)a- governing the GP prior over each weight function. The joint log- proobabbayyl of the observations and the decoding weights is given by:
log po (D, w())
w (yn) Xn -2n)
Σ (log |aal+ + w.KKa Ka wa)
where Zn log (Lg exp (w(y) wwy))xx)dy )dy is the normalizing constant; Wa
[wa(v1), wa(v)).... "T(N)) denotes the weight function of input feature d evaluated at the observed outputs {ym, mmzi and Ka E NxN is the corresponding GP covariance matrix over these outputs.
This formulation intractable due to the need to marginalize over weight functions when computing each Zn, To jointly infer the latent functions w(y) and hyperparameters 6, thus introduce an efficient inference scheme based on Variational Inference (VI) (Beal) 2003) Jordan et al. 1999) Hoffman et al. 2013 Blei al., 2017). VI approximates the true intractable posterior distribution p(wix x) with a tractable distribution (w), by optimizing the Evidence Lower Bound (ELBO):
C(0, v) = RB [1og po ) (w) ]
DEL. (1p (wp(w)lpp
where Dxz (gllp) denotes the Kullback-Leible divergence between the variational distribution and the prior (Blei et al. 2017). This ELBO serves as the objective that CMLR maximizes during training.
3.1 RIEMANN INTEGRAL APPROXIMATION OF THE NORMALIZING CONSTANT
To compute the ELBO, approximate the intractable normalizing integral in Zn using Riemann integration. We partition the output range Ω = [ymin ymax] into T uniform bins of width A = (ymax //mn)/T, and let De denote the center of the 1th bin. The integral is then approximated as: T
exp (w(y) Xx dy R
exp
w (T.) Xn
(6)
a
3

---

Published as a conference paper at ICLR 2026
3.2 FOURIER-DOMALI REPRESENTATION OF WEIGHT FUNCTIONS
To improve computational efficiency, we parameterize the decoding weights in the frequency domain, exploiting the fact that the RBF covariance is diagonalized in Fourier basis. This allows efficient inference by eliminating matrix inversions and operating in a reduced-ding deibcedddimmessioall space (Hensman)
2018 Gondur et al. 2024 Keeley et al., 2020 Aoi Pillow, 2017). Let B E N and denote discrete orthonormal Fourier basis matrices evaluated at the sample outputs (u Jn=l and the Riemann grid points {T, J respectively. The decoding weights are then related
al.
B
R 77XM
N.
their frequency-domain representations via:
Wd (vn) (B wa), and Wa (V+) (Bwa),
where Wa = [Wl,ds W2,ds W.ddd are Fourier coefficients for feature d, drawn independently as: W ~ N (0, km,a), kim,d = Dd exp 1- for m = 1, , M,,
(8)
with Pd = V2Tpala, and fm denoting the mth Fourier frequency. By truncating to M < T, N basis functions, we drastically reduce the dimensionality of inference. In this formulation, inference reduces to estimating oow-dimensional Fourier coefficients {wa) (4=)2=1> d== which parameterize the decoding weight functions via a fixed orthonormal basis. This reparameterization turns the original infnite-diimenionall inference problem into tractable finiteddimensional one in the frequency domain.
3.3 STOCHASTIC VARIATIONAL OPTIMIZATION IN THE FREQUENCY DOMAIN
We assume the variational distribution over frequency-domen domain weights is fully factorized across features and frequencies, with each component following univariate Gaussian:
M
((w)
II ((wm.d) q(um.d) N (Jm,, Fm.dd) d=1 mmm
Under this assumption, the expected log-likes likelikoi term in the ELBO (Eq. .V.
5 becomes:
D
N
T
D
B (a) and )
Egg |oggo ( ((xm) 1 [w)
ΣΣ
pa), Fndd BBg
Σ
log
A
Σ
exp
d=1
where the second expectation approximated via Monte Carlo sampling (Robert et al. 1999.. The KL divergence term in the ELBO (Eq. 5), due to the Gaussian assumptions, admits
DEL (po(w))))(())
Σ
Σ (log (km,a/oma (mm.)) + Σ Gm,dd /.m + Σ (am,a)
/o//////mnm /kmm.
= 2
M.D We jointly optimize the variational parameters (1mm,d mm.a) m,d=1 and GP hyperparameters 0
{pa, la) f using stochastic variational inference with mini-batches of size N' < N. thus facilitating the implementation to scale with N (Hoffman et al.) 2013). We used the Adam optimizer (Kingma & Ba 2017) for training, with all scale parameters (including Ø) optimized in log-space to ensure positivity. Extensive simulations (Appendix I show that the proposed inference procedure reliably recovers both the decoding weights and hyperparameters across diverse settings.
Computational scaling: We evaluated the computational efficiency of CMLR across a range of settings for the number of neurons/features D, samples N, and Fourier components M, shown in Appendix A.3] Fig. S3 Training time scales linearly with the number of neurons D: for example, convergence requires roughly 102 seconds for D RS 200 and 10s seconds for D z 2000 on standard laptop (Intel i7 2.4 GHz CPU, 16 GB RAM). Thanks to stochastic variational inference, training time increases only modestly with the number of samples N. In addition, both runtime and decoding performance are largely insensitive to the number of Fourier components M, indicating that relatively small M values are sufficient for accurate estimation. Overall, Fig. [33 shows that CMLR computationally efficient and robust across broad range of dataset sizes and configurations.
3.4 EXTENSION TO MULTI-DIMENSIONAL OUTPUTS
While previous sections focused on scalar outputs, CMLR naturally extends to multi-dimennional settings, enabling multivariate conditional density estimation. For example, suppose the output
4

---

Published as a conference paper at ICLR 2026
variable is wiv-dimensional, y [y(*) (), y We model the GP prior covariance of the decoding weight functions using an anisotropic RBF kernel (Rasmussen & Williams 2006), which assigns separate length scales to each output dimension:
(-2((11))
y - ? ((9)) +
Ui 1 ynn
Ka (ymr Mnn) Pa exp
yn
(10)
This results in the following modified frequency-doman representation of the prior variance (Eq. 8):
2 kimdd Da exp (-2 ((12))
(2)
2
2
for m 1, , M,
(11)
where AJ = 2mpaed ( ((2 and pm f denote the Fourier frequencies in the two output dimensions. Our stochastic variational inference framework generalizes naturally this multi-dimensional setting. We jointly estimate the variational parameters {Um,d, d M,D along with the GP hyperparam-
{pss
D Fig. S4 presents simulation results demonstrating that the proposed
eters 9
d=1 inference procedure accurately recovers the true underlying parameters in the 2D stimulus setting.
4 OUTPUT PREDICTION USING THE TRAINED CMLR MODEL
Once the CMLR model is trained, the learned parameters can be used to predict outputs for test data at any desired resolution. For target resolution (ymnn Umax] into total of J classes, where = ceil ((Umax ymmi)) /6). Let ün denote the center of the jth bin. Based on this uniform grid, we construct the corresponding decoding Fourier basis matrix Bdec IR Jx and compute the decoding weights for features d 1,..,, D resolution 6 as:
uniformly partition the output range
wa(ū;) = (Bdec x pa);
Given test sample X,, we then compute the posterior over the output grid via the softmax function:
exp (B,) xn.
p (Yn = ÚJ X,, w()))
E/= exp (w(ū))) Tx,))
This yields the full posterior conditional distribution. We then decode the output using either the posterior mcan or the posterior mode, depending on the task objective. The posterior mean yields smooth, high-resolution predictions:
Umnan
Eū; j=1
Új PVv = رs Xn. w(y)).
while the posterior mode selects the most likely output:
Tmode
arg max
p (yn- = ūj Xn, w(y)
jE(I:T)
The mean decoder is particularly well-suited for regression-style tasks, whereas the mode decoder preferable for minimizing classification error. Note that we can set N to be arbitrarily small to achieve a desired level of resolution in output space.
5 APPLICATIONS TO NEURAL DECODING
We evaluated CMLR on four real world neural datasets spanning diverse brain regions and species: mouse VI (Stringer et 2021), monkey V1 (Graf et al. 2011), mouse hippocampal CA1 (Hazon et al. 2022, Jercog & Schnitzer 2022). and monkey motor cortex (Glaser et al.) 2018 2020). For each dataset, we performed 5-fold cross-validation. In each fold, the data were split into 80% training and 20% test sets. For DNN and XGBoost, we further reserved 20% of the training set (i.e., 16% of the total data) as a validation set for hyperparameter tuning via Bayesian optimization (Gardner et al. 2014). For CMLR and the other baselines, hyperparameters were fixed per dataset. The selected CMLR design parameters and their rationale are provided in Appendix A.3] Final performance was averaged over the test sets across all folds. We compared CMLR with the following four baselines.
5

---

Published as a conference paper at ICLR 2026
Naive Bayes (NB): We implemented continuous- variant of the Naive Bayes decoder that shares the same Gaussian process (GP) priors over decoding weight functions wa(y) as CMLR, using RBF kernels and Fourier-domain inference. However, unlike CMLR, NB assumes conditional independence of neural responses across features given the output, yielding the following likelihood:
p ((x))== (w(w))
II p p( ((nnalla(()))
d=1
We used Gaussian observation models for calcium imaging data and Poisson observation models for spike count data, following Greenidge et al. (2024). This formulation preserves the functional interpretability of weights as tuning curves and supports CDE, but ignores correlations across neurons.
Flexible nonparametric conditional density estimation (FlexCode): FlexCode is a state-of-the-art nonparametric CDE method that reformulates CDE as a series expansion problem, estimating the basis coefficients via regression (Izbicki & Lee 2017). We use the publicly available implementation from Tzbicki Lee (2017) with Random Forest regression to estimate coefficients.
Extreme Gradient Boosting (XGBoost): XGBoost is a widely used implementation of gradient- boosted decision trees (Chen & Guestrin 2016). We used the implementation from Glaser et al. (2020), which was specifically tuned for neural decoding and used multiple adjacent time bins as input. Similar to Glaser et al.] 2200), we optimized hyperparameters using Bayesian optimization (Gardner al. 2014) based on validation-set R², searching over tree depth, number of trees, and learning rate.
Deep neural network (DNN): DNNs consist of multiple layers of nonlinear transformations that map inputs outputs (Orbach 1962; Goodfellow al. 2016). We used the architecture and training procedure from Glaser et al. 2020). which was specifically designed for neural decoding. Following this work, we tuned the number of hidden units, dropout rate, and number of epochs using Bayesian optimization based on validation R2, and we used multiple adjacent time bins as input.
5.1 PRIMARY VISUAL CORTEX (VI): DECODING DRIFTING GRATING ORIENTATION
We first applied our method to three two-photon calcium imaging datasets recorded from mouse
primary visual cortex during drifting grating stimuli (Stringer et al. 202T). ). Each dataset included D=11311-20616 neurons (input features) and N 4282 4469 tnals, with stimulus orientations (outputs) uniformly sampled from [0, 2T).. We used the same features x as in SStrnnger et al. (2021).
Fig. 2 summarizes the results. Fig. A compares decoding weights from the CMLR and Naive Bayes models for selected neurons, normalized to [0, 1] in amplitude. CMLR weights appear smoother, consistent with the larger inferred length scales shown in Fig. 2 We next assessed mean absolute circular error across decoding grid resolutions, controlled by the number of decoding classes J (Fig. 2C). CMLR consistently outperformed FlexCode and Naive Bayes, with performance gains saturating beyond J = 5000. We also compared CMLR to XGBoost and DNN decoders (Fig. (DD-)). These models are omitted from Fig.2C because point-estimate regressors lack conditional density and do not naturally admit resolution-dependent analysis. Fig. P2 shows decoded versus true orientation, with CMLR predictions clustering tightly around the identity line, with most large errors occurring near 180", reflecting underlying bimodality. Fig. 2 shows that CMLR achieves the lowest decoding error (mean = SD: 3.1 = 9.3", median: 2.1 Inter-Quartily range (IQR): 2.7%), followed by FlexCode (3.2 + 5.5%. 2.2", 2.8%). Naive Bayes (4.9 + 10.8", 3.3", 4.4%), XGBoost (13.6 2 6.9 11.19), and DNN (18.3 + 23.6". 11.99, 16.5)). Moreover, the scalability analysis in Appendix A.6 (Fig. (55) shows that CMLR experiences only modest reductions in accuracy as D and N decrease, with an even larger performance gap over XGBoost and DNN in low-data settings.
Next, we applied our method to five electrophysiological datasets from monkey primary visual cortex, recorded under drifting grating stimuli (Graf et al., 201T). Fach dataset contained spiking activity from D= = 113-148 neurons across 72 discrete grating orientations (spaced at 5" intervals), with a total of N = 3600 trials. Although CMLR is designed for continuous outputs, it naturally accommodates discrete settings without modification. Using the same inference settings before, we found that CMLR and FlexCode consistently outperformed all baseline models (Fig. (S6).
As correlation-awan decoder, CMLR stands in contrast to the correlation-bline Naive Bayes model (Greenigee et al. 2024), reinforcing the importance of modeling noise correlations for accurate decoding in VI. It is also noteworthy that the discrete MLR model in Grrenidge et al. (2024) can

---

Published as a conference paper at ICLR 2026
A
Inferred weights
B
Length-scale
C
Error vs decoding resolution
comparisons
mean abs eircular error (9)
CMLR
Naive Bayes (NB)
CMLR
FlexCodo
X Naive Bayes (NB)
Naive Bayes
0.0
2m
10'
10
y
CMLR
(number of decoding classes)
D
True decoded
E
Error comparisons
CMLR
3.1 3.215.5
1.9 13.6 28.4
decoded orientation
100
FlexCode
abs circular error
NB
XGBoost
DNN
Ideal
true orientation 300
CMLR
FlexCode
NB
XGBoost
DNN
Figure 2: Application of the proposed method to mouse VI data (data from Stringer et al. (002) (A) Decoding weights inferred by the CMLR and Naive Bayes (NB) models for selected neurons, normalized to [o, 1]. (B) Scatter plots comparing inferred length scales by CMLR model versus NB model for individual neurons. (C) Mean absolute circular error function of the number of decoding classes J for CMLR, FlexCode, and NB. (D) Scatter plots of decoded versus true stimulus orientations across five models: CMLR, FlexCode, NB, XGBoost, and DNN. (E) Box plots comparing the absolute circular decoding error across all methods. Mean errors + standard deviation indicated as insets.
be seen as a special case of CMLR the low-J regime. As increases, CMLR enables principled evaluation of resolution-deeendent performance in the high-resolution limit, rather than relying on arbitrary discretization. This especially relevant for continuous or circular variables, such as orientation, where fine- gne-grained distinctions are behaviorally meaningful. The learned decoding functions (a( act like tuning curves, providing smooth and interpretable weights that link individual neurons to specific stimuli or behavioral outputs. This makes CMLR not only a powerful decoder but also a useful tool for probing population coding in V1.
5.2 HIPPOCAMPUS CA1: DECODING POSITION DURING SPATTAL NAVIGATION
We next applied CMLR to eight calcium imaging datasets recorded from pyramidal neurons in the mouse hippocampus CAI region, while the animals navigated 120-cm linear track to collect water rewards (data from Hazon et al. (2022); Jercog & Schnitzer (2022)). Each dataset contained D 151-497 neurons (input features). and we used the same preprocessed neural features
in Hazon et al. Jercog & Schnitzer (022). downsampled by factor of 10, yielding N = 3600 3600-5524 samples. Position outputs were normalized to [0. 1] for ease comparison.
Fig. summarizes the results. Panel A shows that CMLR consistently outperformed FlexCode and Naive Bayes mean absolute error across decoding grid resolutions. This reinforces the benefit of modeling population-level structure in hippocampal circuits. Panel B displays scatter plots of true versus decoded positions for all models, with CMLR predictions aligning most closely with the identity line. Panel C compares Euclidean decoding errors: CMLR achieved the lowest error (mean + SD: 0.15 = 0.31, median: 0.01, IQR: 0.09), followed by FlexCode (0.16 0.30, 0.01, 0.07). Naive Bayes (0.16 + 0.31, 0.02, 0.06), XGBoost (0.16 + 0.13, 0.11. 0.16), and DNN (0.18 + 0.16, 0.15. 0.19). Together, these findings demonstrate the flexibility and robustness of CMLR across continuous decoding tasks in hippocampal circuits.
7

---

Published as a conference paper at ICLR 2026
A Error vs decoding resolution
B
True vs decoded
C
Error comparisons
0.19
CMLR
decoded position
154.31 165.3 J64.31 164.13 48
FlexCode
abs docoding error
Naive Bayes (NB)
10P
true position
CMLR FlexCode NB
XGBBoost DNN
(number of decoding classes)
FlexCode DNN --. Ideal
XGBoost
Figure 3: Application of the proposed method to mouse hippocampal data (Hazon et al. 2022 Jercog
Schnitzer, 2022). (A) Mean absolute error for CMLR, FlexCode, and Naive Bayes (NB) function of the number of decoding classes J. (B) Scatter plots comparing true versus decoded spatial positions for CMLR, FlexCode, NB, XGBoost, and DNN. (C) Box plots comparing absolute decoding errors across models (insets: mean + standard deviation).
5.3
MOTOR CORTEX: DECODING 2D VELOCITY IN REACHING TASK
We then applied CMLR to a dataset recorded from the motor cortex of monkey performing target-reaching task (data from Glaser et al. (2018 2020 ) In this task, the monkey controlled cursor on screen using a manipulandum (Glaser et al. 2018), and our goal was to decode the cursor's two-dimensional (2D) velocity in the x and y directions. We used the same preprocessed neural features x in Glaser et al. (2020). The dataset comprised N 25,299 samples over 21-minute session, with spiking activity from D = 164 neurons.
Fig. 4 summarizes the results. Panel A shows inferred decoding weights for seven example neurons, revealing diverse spatial tuning profiles in the 2D velocity space, with all neurons shown in Fig. [S7 Additional analyses in Fig. (S8 show that the inferred weight functions reflect empirical spike-velocity structure and remain smooth and coherent even when trained on restricted portions of the velocity space. Panel B compares true and predicted velocities for all models; CMLR aligns closely with the identity line. Panel C presents Euclidean decoding errors and reports the coefficient of determination (R²): CMLR achieves strong performance (0.53), outperforming FlexCode (0.35) and Naive Bayes (- 0.43), and approaching XGBoost (0.55) and DNN (0.58). The higher scores of XGBoost and DNN are expected, given the large size of this dataset, which favors high-capacity nonlinear models. CMLR
A
Inferred decoding weights
neurce
112
147
veloeity
201
-20
-25
25
velocity
C
Error comparisons
B
True vs decoded
x velocity
y velocity
= 0.53
0.35
-0.43
0.55
0.58
CMLR
504
20-
FlexCode
Ruelidean error
decoded
NB
XGBooot
DNN
Ideal
tru
25
FlexCode
NB
XGBoost
DNN
Figure 4: Decoding 2D cursor velocity from motor cortex data (Glaser et al. 2018 2020). (A) Decoding weights inferred by CMLR for seven representative neurons (full set shown in Fig. S7 (B) Scatter plots comparing true versus decoded z (left) and y (right) cursor velocities for all methods: CMLR, FlexCode, Naive Bayes (NB), XGBoost, and deep neural networks (DNN). (C) Box plots of Euclidean decoding error across all methods (insets report the R² score).
8

---

Published as a conference paper at ICLR 2026
nevertheless remains competitive while also providing full conditional densities and interpretable tuning functions, capabilities that are not available in point-estimate regressors. Overall, these findings demonstrate that CMLR extends naturally to multivariate-output CDE and continuous decoding.
5.4 CALIBRATION OF PREDICTIVE UNCERTAINTY
Finally, we evaluated the calibration of the predictive posteriors from the two CDE methods, CMLR and FlexCode, using PIT histograms and quantile calibration curves (Fig. S9 see Appendix A.10 for details). Across mouse VI, macaque VI, and mouse hippocampal CAI, CMLR consistently produces well-calibrated posteriors: PIT histograms are close to uniform, quantile calibration curves track the identity line, and PIT values vary smoothly with decoding error. FlexCode, in contrast, shows systematic miscalibration, with peaked or multimodal PIT distributions, under-coverage in quantile calibration, and irregular PIT-error relationships. These results indicate that CMLR provides more reliable and interpretable uncertainty estimates, which is important in neural decoding applications where posterior uncertainty guides scientific interpretation and downstream decisions, such as identifying ambiguous stimuli or estimating decoding confidence.
5.5 RUN TIME COMPARISONS
Table reports average training times for all methods on a 2.4 GHz Intel i7 CPU with 16 GB RAM (all without parallelization). CMLR trains within minutes to a few hours, depending on dataset size, with run times comparable to FlexCode and faster than Naive Bayes. XGBoost and DNN are generally faster, but they do not provide full conditional densities or calibrated uncertainty. Overall, CMLR achieves competitive computational efficiency while offering richer probabilistic outputs.
Table Run time comparisons for all methods. Values are averaged over folds.
Dataset (# neurons)
N (# trials)
CMLR
FlexCode
NB
XGBoost
DNN
Mouse V1
11-21 k
4.3-4.5 k
3.3 h
2.3 h
22 h
4.4 h
33 m
Monkey V1
113-148
3.6 k
5.5 m
25
2.6 m
2 m
1.9 m
Mouse CA1
151-497
3.6-5.5 k
10 m
47
4 m
4 m
2 m
Monkey MC
164
k
10
11 m
30 m
23
m
CONCLUSION AND DISCUSSION
We introduce the Continuous Multinomial Logistic Regression (CMLR) model, novel exponential- family framework for scalable, nonparametric conditional density estimation (CDE). CMLR extends multinomial logistic regression (MLR) to continuous outputs by replacing discrete class weights with smooth, output-indexed functions drawn from Gaussian process (GP) priors, yielding normalized densities while preserving MLR's s additive, interpretable structure. For efficient training, we develop memory- efferreet stochastic variational inference algorithm in the Fourier domain, leveraging GP stationarity for kernel diagonalization and basis truncation. Applied to large-scale neural datasets from mouse and monkey V1, hippocampus CA1, and motor cortex, CMLR outperforms Naive Bayes, XGBoost, deep neural networks, and the leading CDE method FlexCode, while providing calibrated posteriors and interpretable structure. Comparisons with Naive Bayes demonstrate that explicitly modeling correlation structure substantially improves decoding performance across all brain areas.
To our knowledge, CMLR represents the first application of CDE to neural decoding, enabling flexible estimation of full posterior distributions over behaviorally or perceptually relevant variables. While modern high-capacity nonlinear models may achieve stronger predictive accuracy in very large-data regimes, they typically require extensive hyperparameter tuning and offer limited interpretability. CMLR plays complementary role: yields fully probabilistic and well-calibrated conditional densities, provides access to the full posterior rather than point estimates, is data-efficient with stable performance across datasets, and remains easy to train because it requires only a few design parameters and no dataset-specific hyperparameter tuning. In addition, CMLR produces transparent tuning functions that can be directly visualized and compared across neurons and conditions, in the same manner classical "sstems-neurosciens black-box box models. We therefore view CMLR as a practical and interpretable baseline or diagnostic model that complements more complex nonlinear approaches in modern neural decoding settings.
tuning curves, property not typically supported by

---

Published as a conference paper at ICLR 2026
Connections to Prior Work: CMLR builds rich body of research in nonparametric Bayesian inference using Gaussian processes (Williams & Barber, 1998; Girolami & Rogers 2006] Bishop
2006), but is conceptually distinct from standard GP regression (Rasmussen & Williams] [006]] Chan
2013) and GP classification (Liu et al. 2022). In GP regression, the output is modeled as GP function of the input, yielding Gaussian predictive distributions. In GP classification, latent GP is passed through a nonlinear link to produce class probabilities over discrete outputs. Both approaches lack the flexibility to represent rich or structured conditional densities over continuous outputs, and related methods that approximate uncertainty for point estimators, such as Laplace approximation (Daxberger et al., 2022), also do not provide full conditional densities. In contrast, CMLR places GP priors on weight functions defined over the output space, enabling flexible estimation of full conditional densities. CMLR is also closely related to logistic Gaussian processes Tokdar et al. 2004) and early unconditional density models such as the Gaussian process density sampler (Murray ct al. 2008), but generalizes them into a conditional, feature- desture-dccomposed framework that preserves the additive structure and interpretability of MLR. As a CDE model, CMLR complements several
existing approaches, including kernel-ratio estimators (Bashtannyk & Hyndman 2001 Holmes
et al. 2012 (Sugiyaaa et al. (2010), mixture density networks (Bishop] 1994), neural neural-kernel mixture
models (Ambrogioni et al. 2017). conditional normalizing flows (Papamakarios et al., 2017), logistic
Gaussian process partition models (Payne al. 2019), Lindsey's Method (Gao & Hastie 2022). and
histogram trees (Yang& & van Lecuwen (024)) However, many of these methods face limitations
in scalability, interpretability, or statistical robustness in high-ddmensional settings. CMLR offers a nonparametric and additive alternative that provides interpretable weight functions, supports multiple output dimensions, and achieves computational scalability through structured inference. CMLR also has conceptual links to neural operator models, particularly those defined over continuous input-output mappings (Kovachki et al. 2023 Li et al. 202I). To scale CMLR to large datasets, we leverage sparse Fourier-domain representations developed for variational GP inference (Hensman ct al. 2018 Keeley et al. 2020; Gondur et al. 2024). closely related in motivation and formulation to Fourier Neural Operators (Li et al. 202 l). Our model also builds directly on the discrete multinomial
logistic regression framework for neural decoding in Greenidge et al.] (2024). extending it to handle continuous outputs in principled and computationally efficient manner.
Limitations and Future Directions: While CMLR provides a unified and scalable framework for CDE, several directions remain open for future work. First, unlike XGBoost and DNNs that can implicitly capture temporal continuity by incorporating activity from multiple adjacent time bins (Glaser et al. 2020), CMLR does not yet include explicit output-space priors, such as temporal or spatial smoothness, which could be useful in navigation or motor decoding. Integrating structured priors into the conditional likelihood, or combining CMLR with latent dynamical models (Park) et 2015) Damianou et al. (011]. may enhance temporal generalization. Second, the current CMLR model treats each decoding weight independently; multivariate GP priors (Bonilla et al., 2007 Keeley et al., (220)) could share structure across neurons time and potentially improve generalization. However, such coupling would obscure neuron-specific tuning curves and thus compromise interpretability. Third, while fixed Fourier-domain bases and RBF kernels enable scalable inference, more flexible alternatives such as adaptive basis functions (Evans & Nair 2018) or spectral mixture kernels (Paciorek & Schervish] (203) could further enhance accuracy, though it remains to be tested whether such extensions retain the same computational scalability. Fourth, while additivity over inputs ensures scalability, convexity, and interpretability, incorporating low-rank or kernelized interaction models (Duvenaud et al., (201I) could enhance model expressiveness by capturing higher-order feature dependencies. Fifth, although we used Riemann integration to support tractable variational inference, advanced numerical techniques such as general quadrature methods (Hildebrand 1987), polynomial approximations (Trefethen] 2019). or adaptive binning (Wand) 1997) may offer improvements in accuracy or efficiency. Sixth, hybrid approaches that combine inducing- point methods with Fourier domain representations (Hensman et al. 2018) could further improve computational scalability. Finally, although we focused here sensory and motor decoding, CMLR is broadly applicable to other neural decoding tasks, including decoding motor intention (Tam et al.. 2019), speech (Chen et al. 2024).. spatial attention (Smith al.) 2013), and decision variables (Baeg et al. 2003). Beyond neuroscience, CMLR may also benefit other CDE applications such as head pose estimation MMurphy-Chutorice & Triiedi 2009). time-to-event modeling (Gensheimer & Narasimhan 2019), climate prediction (Rasp et al. 2018). photometric redshift estimation (Dalmasso
al. 2020), and geolocation and trajectory forecasting (Rhinehart et al. 2019). Together, these directions position CMLR as a versatile CDE framework for neural decoding and beyond.
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENTS
This work was supported by the NIH T32 institutional training grant (T32MH065214), the Simons Collaboration on the Global Brain (SCGB AWD543027), the NIH BRAIN Initiative (R01DA056404, RO1EB02694), the National Eye Institute (NEI) the National Institutes of Health (NIH) (R01EY033064), and U19 NIH-NINDS BRAIN Initiative Award (5U19NS 1237166.. We also thank C. Stringer, M. Michaelos, D. Tsyboulski, S. Lindo, and M. Pachitariu for providing the publicly available mouse V1 datasets; A. B. A. Graf, A. Kohn, M. Jazayeri, and J. A. Movshon for providing the monkey VI datasets; P. Jercog and M. Schnitzer for providing the publicly available mouse hippocampuu CAI datasets; and I. Glaser, M. G. Perich, P. Ramkumar, L. E. Miller, and K. Kording for providing the publicly available monkey motor cortex dataset.
REFERENCES
Douglas G Altman and Patrick Royston. The cost of dichotomising continuous variables. BMJ, 332 (7549):1080, May 2006.
Luca Ambrogioni, Umut Güçlü, Marcel van Gerven, and Eric Maris. The kernel mixture network: A nonparametric method for conditional density estimation of continuous random variables, 2017.
jRR_sstss://axxii/ org/ abs/1705. 07111
Mikio C. Aoi and Jonathan W. Pillow. Scalable bayesian inference for high-dimensional neural receptive fields. bioRxiv, 2017. doi: 10.1101/212217.
Baeg. B Kim, K Huh, Mook-Jung, H Kim, and M W Jung. Dynamics of population code for working memory in the prefrontal cortex. Neuron, 40(1):177-188, September 2003.
David M. Bashtannyk and Rob Hyndman. Bandwidth selection for kernel conditional density
estimation. Computationnl Statistics & Data Analysis, 36(3):279-298 2001. ISSN 0167-9473. doi: http::/wwuugg////0006666//777200 URL htt.tttts///// wwwp://wwwwinnncceee.....
com/ http://ccwe/aarrreee//////////1/
M. J. Beal. Variational algorithms for approximate Bayesian inference. PhD thesis. University of London, University College London (United Kingdom), 2003. URL hittps://discovte
uel ww..ac.kkkka
Christopher M. Bishop. Mixture density networks. Working paper, Aston University, 1994.
Christopher M Bishop. Patlern Recognition and Machine Learning. Springer, 2006.
David M Blei, Alp Kucukelbir, and Jon D McAuliffe. Variational inference: A review for statisticians. Journal of the American Statistical Association, 112(518):859-877, 2017. doi: https://do..or//10
1080/01621459.222 1087.7785773.
Edwin Bonilla, Kian Chai, and Christopher Williams. Multi-task gaussian pro- cess prediction. In J. Platt, D. Koller, Y. Singer, and S. Roweis (eds.), Ad- vances in Neural Information Processing Systems, volume 20. Curran Associates, Inc., 2007. URL
www.eripps.ccppppppeeeeeee Files/paper:2007/ Bec (tteee5588d33322991bb333 fi wtte:/6668220fFFF 1 pappeerrppbbppp..p......
E. Brown, L. Frank, D. Tang. M. Quirk, and M. Wilson. A statistical paradigm for neural spike train decoding applied to position prediction from ensemble firing patterns of rat hippocampal place cells. Journal of Neuroscience, 18:7411-7422 1998.
Antoni B. Chan. Multivariate generalized gaussian process models, 2013. URL hhttps:/aaat
htt//bb//3311. 0360
Tianqi Chen and Carlos Guestrin. Xgboost: A scalable tree boosting system. In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD '16, pp. 785-794, New York, NY, USA, 2016. Association for Computing Machinery. ISBN 9781450342322. doi: 1111445/2939672229997755
Utthttps: //odi////////./// wwt.do.o//...11.55551100
2939 08339672.255 2939778

---

Published as a conference paper at ICLR 2026
Xupeng Chen, Ran Wang, Amirhossein Khalilian-Gouus Leyao Yu, Patricia Dugan, Daniel Friedman, Werner Doyle Orrin Devinsky, Yao Wang, and Adeen Flinker. A neural speech decoding framework leveraging deep learning and speech synthesis. Nature Machine Intelligence. 6(4):467-488, Apr 2024. ISSN 2522-5839. doi: 10. 10100886422560444088223-3 URL https
7/doi.org C. T038 4256/842256-222-008832
N. Dalmasso, T. Pospisil, A.B. Lee, R. Izbicki, P.E. Freeman, and A.I. Malz. Conditional density estimation tools in python and with applications to photometric redshifts and likelihood cosmological inference. Astronomy and Computing, 30:100362, 2020. ISSN 2213-1337 doi: com/
https://doi..rr/10... https//wwooooggg//ooooo///aaaama asc.m. 2019..000221 URL htt.::::/www.sciennee Bciiencee/arriic.ee///////22...1
Andreas Damianou, Michalis Titsias, and Neil Lawrence. Variational gaussian process dy- namical systems. In J. Shawe-Taylor, R. Zemel, P. Bartlett, F. Pereira, and K.Q. Wein- berger (eds.), Advances Neural Information Processing Systems, volume 24. Curran Asso- ciates, Inc., 2011. URL httpst httptppss///rroreeeiins.sssse/ wwtpriipssccppppppeeeeeeee wc/paper.p........ cattser_fllssspppper Paperr.
A. Dawid. Present position and potential developments: Some personal views: Statistical theory: The prequential approach. Journal of the Royal Statistical Society. Series A (General), 147(2): 278 292, 2025/09/10/ 1984. doi: 10.2307/298 1683. URL http:::/tp////tttt http://oo...og//....0..03
2981 1 683
Full publication date: 1984.
Erik Daxberger, Agustinus Kristiadi, Alexander Immer, Runa Eschenhagen, Matthias Bauer, and Philipp Hennig, Laplace redux effortless bayesian deep learning, 2022. BtRhhtps::/ hattps://arxix:
orq/abs 2106. 14806
David Duvenaud, Hannes Nickisch, and Carl Edward Rasmussen. Additive gaussian processes, 2011. URL
wwwpxxiiv...g.a..a.. abs/ /e////111.... 4394
Trefor Evans and Prasanth Nair. Scalable Gaussian processes with grid-structured eigenfunctions (GP-GRIEF). In Jennifer Dy and Andreas Krause (eds.), Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 1417-1426. PMLR, 10-15 Jul 2018. URL
www.nngsmm..rr......m.m ml wtt.g/sss.......ss.......
evansl8a.
Zijun Gao and Trevor Hastie. Linede: conditional density estimation via lindsey's method. Journal of machine learning research, 23(52):1-55, 2022.
Jacob Gardner, Matt J Kusner, Zhixiang Eddie Xu, Kilian Q Weinberger, and John P Cunningham. Bayesian optimization with incquality constraints. In ICML, volume 2014, pp. 937-945, 2014.
Michael F Gensheimer and Balasubramanian Narasimhan. scalable discrete-time survival model for neural networks. PeerJ, 7:e6257, January 2019.
A. Georgopoulos, J. Kalaska, R. Caminiti, and J. Massey, On the relations between the direction of two-dimensional arm movements and cell discharge in primate motor cortex. Journal of Neuroscience, 2:1527-37, 1982.
Mark Girolami and Simon Rogers. Variational bayesian multinomial probit regression with gaussian process priors. Neural Computation, 1898):1790-1817, 2006. doi: 10. 10.11 1890
Joshua Glaser, Matthew G. Perich, Pavan Ramkumar, Lee Miller, and Konrad P. Kording. Population coding of conditional probability distributions in dorsal premotor cortex. Nature Communications,, 9(): May 2018. ISSN 2041-1723 doi: 10. 4117/444166--1---44444 URL httph:://///m.. /ittps://doi.. /ww.///oo///////////////o./ org 0. 1 038/038800 84144444677774 018 04062-6
Joshua I. Glaser, Ari S. Benjamin, Raeed H. Chowdhury, Matthew G. Perich, Lee E. Miller, and Kon- rad P. Kording, Machine learning for neural decoding. eNewro. 7(4), 2020. doi: 10.1UR//EEEURO. 0506- 19.2020. URL https: www. eneuro. 2020
0506-19
12

---

Published as a conference paper at ICLR 2026
Rabia Gondur, Usama Bin Sikandar, Evan Schaffer, Mikio Christian Aoi, and Stephen L. Keeley. Multi-modal gaussian process variational autoencoders for neural and behavioral data. In The Twelfth International Conference Learning Representations, ICLR 2024, Vienna, Austria, May 7-11, 2024. OpenReview.net, 2024. URL https:/ htttps://openrev /tpp:/vveeeevveeeee/ee www. ew. net. /t///foron forum?
aGH43r
lan Goodfellow, Yoshua Bengio, and Aaron Courville. Deep Learning. MIT Press, 2016. http:
www.deep// www.::ddeeeeerinnngiiinn eernniggooonooooooe
Arnulf B. A. Graf, Adam Kohn, Mehrdad Jazayeri, and J. Anthony Movshon. Decoding the activity of neuronal populations in macaque primary visual cortex. Nature Neuroscience, 14(2)233922) 245, Feb 2011. ISSN 1546-1722 doi: 10. 10.1038/mm.2733 URL
2733
C. Daniel Greenidge, Benjamin Scholl, Jacob L. Yates, and Jonathan W. Pillow. Efficient decoding of large-scale neural population responses with gaussian-process multiclass regression, Neural Computation, 36(2):175-22) 01 2024. ISSN 0899-7667. doi: 10. 162/neco a 01600 URL
https::/doi..oo http////o.or/..ooooo//o wtt./////................ 1162/neccocoo 0163C
Omer Hazon, Victor H. Minces, David P. Tomås, Surya Ganguli, Mark J. Schnitzer, and Pablo E. Jercog. Noise correlations in neural ensemble activity limit the accuracy of hippocampal spatial representations. Nature Communications, 13(1))4276, Jul 2022. ISSN 2041-1723, doi: 10.1038/ s41467--22-311544
URLhtttps://daa otg./oo...g.///oooo.o 1038/ 111111101/01117700-22021222
James Hensman, Nicolas Durrande, and Arno Solin. Variational fourier features for gaussian processes, Journal of Machine Learning Research, 18(151):1-52, 2018.
Francis Hildebrand. Introduction Numerical Analysis. Dover Publications, New York, 2nd edition, 1987.
Matt Hoffman, David M. Blei, Chong Wang, and John Paisley. Stochastic variational inference, 2013. URL
wwwp:xxwivoovgooooooo// abs 1206. 7051
Michael Holmes, Alexander G Gray, and Charles Lee Isbell. Fast nonparametric conditional density estimation. arXiv preprint arXiv: 2012.
Heikki Huttunen, Tapio Manninen, Jukka-Pekka Kauppi, and Jussi Tohka. Mind reading with regularized multinomial logistic regression. Machine vision and applications, 2416):131 2013.
Rob J. Hyndman, David M. Bashtannyk, and Gary K. Grunwald. Estimating and visualizing conditional densities. Journal of Computational and Graphical Statistics, 5(4):315-336, 1996. ISSN 10618600. URLhttp: httep://wwer //www//wwww/...../ jstor. org. wtg.:jstoooootttee.ooo.ooo htt.://orrrggstaabe/eee//ee
Rafael Izbicki and Ann B. Lee. Converting high-dimensional regressiooooiiirrreeeeesssooo hign-dimensional con- ditional density estimation. Electronic Journal of Statistics, 11(2):2800 2831, 2017. doi:
10.1214/1-EEJS138
/Rtpstts:://dol http:///ooooorg/1//// 1.c...21444 1ttp://11...11411-1111113330
Pablo Jercog and Mark Schnitzer. Functional calcium imaging recordings hippocampal area cal from mice performing linear track task. 2022. Available at CRCNS.org.
6080
Michael Jordan, Zoubin Ghahramani, Tommi S Jaakkola, and Lawrence K Saul. An introduction to variational methods for graphical models. Machine Learning, 37(2):183-33,, 1999. doi: https://doi.org/10.. http://oo..oooo..000022010000000
Stephen Keeley, Mikio Aoi, Yiyi Yu. Spencer Smith, and Jonathan W Pillow. Identifying signal and noise structure in neural population activity with gaussian process factor mod- els. In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (eds.). Advances in Neural Information Processing Systems, volume 33, pp. 13895.13805 Curran Asso- http:///....pppppaaaaaaaaeaa//// files//aper/
ciates, Inc., 2020. htL c
13

---

Published as a conference paper at ICLR 2026
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization, 2017. URL
pattps:/arxivy/ org/ wtt/abs//////oooooeo..ooo...esss abs/1412 6980
Nikola Kovachki, Zongyi Li. Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications to pdes. Journal of Machine Learning Research, 24(89): 2023.
Volodymyr Kuleshov, Nathan Fenner, and Stefano Ermon. Accurate uncertainties for deep learning using calibrated regression. In International conference on machine learning, pp. 2796-2804. PMLR, 2018.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar, Fourier neural operator for parametric partial differential equations, 2021. URL Attps: www.wiiwww..aaaaa..a
org
08895
Haitao Liu, Yew-Soon Ong, Ziwei Yu, Jianfei Cai, and Xiaobo Shen. Scalable gaussian process classification with additive noise for non-gaussian likelihoods. IEEE Transactions on Cybernetics, 52(7):54425884, 2022. doi: 10. 10.1100 (0.9/CCBB.2220.30111
Erik Murphy-Chutorian and Mohan Manubhai Trivedi. Head pose estimation computer vision: survey. IEEE Transactions Pattern Analysis and Machine Intelligence, 31(4):607-626, 2009.
doi: 10.1109/PAMM..200..UU
lain Murray, David MacKay, and Ryan P Adams. The gaussian process density sampler. In D. Koller, D. Schuurmans, Bengio, and Bottou (eds.), Ad- vances in Neural Information Processing Systems, volume 21. Curran Associates, Inc.. 2008. URL File/
files//pperr20008 Paper. .ppp
Nirenberg and P. Latham. Decoding neuronal spike trains: how important are correlations? PNAS, 1000:338-7353, 2003.
Farnaz Nojavan A., Song S. Qian, and Craig Stow. Comparative analysis of discretization methods in bayesian networks. Environmental Modelling & Software, 87:64-71, 2017. ISSN 1364-8152.
doi: http:///nnneaarreeee////////1113
10.007. URL hou.https::. /www. sc enced rect
com/ sc
Orbach. Principles of neurodynamics. perceptrons and the theory of brain mechanisms. Archives of General Psychiatry, 7(3):218-219, 09 1962. ISSN 0003-990X. doi: 10.1001/
(tttpsyy._9222001200300600600000 URL 01720030064010
Christopher Paciorek and Mark Schervish. Nonstationary covariance functions for gaus- sian process regression. In S. Thrun, L. Saul, and B. Schölkopf (eds.), Advances in Neural Information Processing Systems, volume 16. MIT Press, 2003. URL http://eeeeerrrose////ppeeeee//e 33580550004555544440333------003 pdff
Paninski and JP Cunningham. Neural data science: accelerating the experiment-analyste theory cycle in large-scale neuroscience. Current Opinion Neurobiology, 50:232-241, 2018. ISSN 0959-4388. doi: htt.a//cc.......a.......///////.
URL https:/ //www/
scciencdiiree wtteeeettcccccccieeeeeeeeeeeeeee nologies.
Neurotech-
George Papamakarios, Theo Pavlakou, and Iain Murray. Masked autoregressive flow for density estimation. In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 30. Curran Asso- ciates, Inc., 2017. URL
http://aaprfffflel.aeeeeeee files/paper/ pdf
2017 0077077000/001222222222222283353
Mijung Park, Gergo Bohner, and Jakob H Macke. Unlocking neural population non-stationarities using hierarchical dynamics models. In C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and
14

---

Published as a conference paper at ICLR 2026
R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 28. Curran Asso- htt.://add......e.o.o..o.ooo//po files/paper/
ciates, Inc., 2015. URL
Payne, N Guha, Y Ding, and B K Mallick. conditional density estimation partition model using logistic gaussian processes. Biometrika, 107(():173-19))) December 2019.
Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian processes for machine learning. Adaptive computation and machine learning. MIT Press, 2006. ISBN 026218253X.
Stephan Rasp. Michael S. Pritchard, and Pierre Gentine. Deep learning to represent subgrid pro- cesses in climate models. Proceedings of the National Academy of Sciences, 115(39)9988-- 9689, 2018. doi: 10. 10.1073/pnas. 1810286115. hallht:p:://ett wwtp pnas org/ dos//////l 1073/pnas 18 018880266600 6115
Nicholas Rhinehart, Rowan McAllister, and Sergey Levine. Deep imitative models for flexible inference, planning, and control,
DRLLhttps.//xxi.oor///////
06544
Jaakko Riihimäki and Aki Vehtari. Laplace Approximation for Logistic Gaussian Process Density Estimation and Regression. Bayesian Analysis, 9(2):425 448, 2014. doi: 10. 10.1214/14-B/87 URL
Christian P Robert, George Casella, and George Casella. Monte Carlo statistical methods, volume 2. Springer, 1999.
Srikanth Ryali, Kaustubh Supekar, Daniel Abrams, and Vinod Menon. Sparse logistic regression for whole-brain classification of fmni data. Neurolmage. 51(2):552-764 2010.
Bernhard Schölkopf and Alexander J. Smola. Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond. The MIT Press, 12 2001. ISBN 9780262256933.. doi: 00.7555//ipprres/411500000000 URL
001 00001
David V. Smith, John A. Clithero, Christopher Rorden, and Hans-Otto Karnath. Decoding the anatomical network of spatial attention. Proceedings of the National Academy of Sciences, 110 (4):1518-1523, 2013. doi: 10. 12tp://aas.111161111 abs/10.1073
httpstss: wwwww pnww..pns... pnas. org/doi/
12t.s.1111110 1210126110
Sutao Song, Gongxiang Chen, Yu Zhan, Jiacai Zhang, and Li Yao. Comparative study of two sparse multinomial logistic regression models in decoding visual stimuli from brain activity of fmri. In Medical Imaging 2014: Image Processing, volume 9034, pp. 447-452. SPIE, 2014.
Carsen Stringer. Michalis Michaelos, Dmitri Tsyboulski, Sarah E. Lindo, and Marius Pachitariu. High- precision coding visual cortex. Cell, 184(10)2767-2778.11) May 2021. ISSN 0092-8674. doi: 1110066/.ce11.2221.....0 URL http::p:://.....//11....o.
2021. 03. 042
Masashi Sugiyama, Ichiro Takeuchi, Taiji Suzuki, Takafumi Kanamori, Hirotaka Hachiya, and Daisuke Okanohara. Conditional density estimation via least-squares density ratio estimation. In Yee Whye Teh and Mike Titterington (eds.), Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, volume 9 of Proceedings of Machine Learning Research, pp. 781-788, Chia Laguna Resort, Sardinia, Italy, 13-15 May 2010. PMLR. URL http:ss/ppocceedingg..
www..prees//..mm.. /atress/v//eet sug yama 1 Ca html
Wing- kin Tam, Tong Wu, Qi Zhao, Edward Keefer. and Zhi Yang. Human motor decoding from neural signals: a review. BMC Biomedical Engineering, 1(1):22, Sep 2019. ISSN 2524-4426. doi: 10. 1186/s42490-01--00022- URL /tt::://ooi.o .rgg 1ww11118800000 010124222200000022-2-22022099920
Surya Tokdar, Yu M Zhu, and Jayanta K Ghosh. Bayesian density regression with logistic gaussian process and subspace projection. Bayesian Analysis, 1(1), 2004.
Lloyd N. Trefethen. Approximation Theory and Approximation Practice, Extended Edition. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2019. doi: 10.1137/1..7816111 1975949. 0T11:/1..711111999
siam. wtt..oog//oo//bs///////// URL htttts::///pppt /ttp://ppppp...aaa.oo///////////
15

---

Published as a conference paper at ICLR 2026
M. P. Wand. Data-based choice of histogram bin width. The American Statistician, 51(1)599 64. 1997. doi: 10. 1tt.://00000110.1199997777700510 URL 155//22//1..1300000011133333333
C.K.I. Williams and D. Barber. Bayesian classification with gaussian processes. IEEE Transactions on Puttern Analysis and Machine Intelligence, 20(12): 1998. doi: 10.1 10.1100/34.735807
Lincen Yang and Matthijs van Leeuwen. Conditional density estimation with histogram trees. Advances in Neural Information Processing Systems, 37:117315-117339 2024.
16

---

Published as a conference paper at ICLR 2026
A APPENDIX
This appendix provides supplementary material that supports, extends, and validates the core contri- butions presented in the main text.
Section A.
presents the simulation study, demonstrating that the proposed inference proce-
dure robustly recovers the decoding weight functions and hyperparameters.
Section A.2
(Fg.
S2
presents full decoding weight visualizations from the simulation study
in Section
reports additional simulation results evaluating inference performance under varying model dimensions, including the number of features (D), sample size (N), and Fourier resolution (M), complementing the analysis in Section A.1
Section A.3 Fig.
Section A.4 (Fig. S4 extends the simulation-based validation to two-dimensional output spaces, demonstrating accurate recovery of decoding weights and GP hyperparameters in multidimensional setting.
Section A.5 describes the CMLR design parameters, provides guidance on how to select them, and Table summarizes the settings used in each real-data application.
count and dataset size on the mouse V1 dataset, and Fig. SS compares its performance and training time with XGBoost and DNN.
Section presents scalability analysis of CMLR with respect
Section A.7 (Fig- 56) compares CMLR performance with other methods the monkey V1 dataset recorded under drifting grating stimuli (data from Graf et al. (201T)), extending the
results in Section 5.1
Section A.8 (Fig. S7 provides complete visualization of decoding weight functions inferred from the monkey motor cortex dataset (extending Fig. 4. highlighting a rich diversity of smooth and interpretable 2D tuning profiles across all recorded neurons.
Section A.9 (Fig. (S8) examines the interpretability and generalization of the inferred decod- ing weights by comparing CMLR weight functions to empirical spike-velocity maps and by analyzing how the inferred functions behave when trained on progressively restricted
subsets the velocity space for the monkey motor cortex dataset (extending Fig. 4
Section A.10 details the procedures used to assess calibration of CMLR and FlexCode posterior densities, including PIT histograms and quantile calibration curves, with results evaluated across datasets (Fig. S9)
Section A.1I declares the limited use of large language models (LLMs) for minor editing of the manuscript.
These results collectively reinforce the accuracy, scalability, and interpretability of the CMLR framework.
17

---

Published as a conference paper at ICLR 2026
A.1 SIMULATION STUDY
To assess the accuracy and robustness of our inference framework, we conducted simulations using synthetic data. We considered continuous output space of orientation angles 0 [0, (0,2 with N 4000 samples. The input consisted of D = 200 features, with decoding weights waly) drawn from Gaussian process (GP) priors. Each weight function was sampled on a fine grid of J points.
(M (Mgr J with GP hyperparameters drawn as la E [0.1, 1.5) (length scales) and pa E [0.5, 2.5] (variances). Feature vectors Xn were sampled independently from a standard Gaussian distribution. ( 1/ were then drawn from the conditional density:
Outputs y/ E {=g
exp w(ygrid) gridd
X
p ggr Xn, w(v)
Σ exp w (Ugrr)) Xn
using a finite grid approximation. IgorithmmsSmmmmrizze summarizes the procedure. Note that this procedure can be made arbitrarily accurate by refining the output grid, and exact sampling is also possible via inverse CDF sampling due the smoothness of the GP-drawn weight functions.
Algorithm 1 Simulation procedure for generating synthetic data ())
Define an output grid {Ugid 1==1 with Ugrid [0, 2Tr
2: for d
D do
3: Sample weight function waly)) from GP prior at (1) with ba [0.1, 1.5), Pa
[0.5, 2.5]
4: end for
5: for = 1 N do
6: Sample feature vector X, ~ (0, Ip)
7: Compute posterior p(y) Xn, w(y)) above
8: Sample observed output p(y Xx, w(y))
9: end for
For inference, we used T 100 bins, MM = 17 Fourier components, a mini-batch size N' 1500, learning rate 0.05, and 3 Monte Carlo samples to approximate the ELBO during training.
Fig. SI summarizes the results. Panel A shows the ELBO trajectory over 580 iterations (104 seconds 2.4GHz Intel i7 CPU with 16GB RAM), exhibiting smooth, monotonic convergence and confirming stable stochastic optimization. Panel compares ground truth and inferred GP hyperparameters across all D 200 features. Inferred values closely track the identity line, with normalized absolute errors of 11 11.35% + 12.14% for length scales and 48.56%+ 54.81% for variances. Larger errors arose for flatter functions (longer length scales larger variances), where identifiability is inherently limited, Panel C shows true vs. inferred decoding weights for 5 representative features (full results in Fig. S2). with mean error across all features of 7.46% + 5.11%. 1, confirming recovery of both coarse- and fine-scale structure.
To further assess robustness, we performed additional simulations under varying conditions. Fig. S3 reports results as we varied the number of input features (D), number of samples (N), and number of Fourier components (M). Performance improves systematically with larger N and M, while method remains computationally scalable and statistically accurate in high-dimensional regimes. Notably, strong performance was achieved with as few as M x 20 Fourier components, supporting the use of compact reqquencydoomain representations reduce dimensionality and prevent overfitting. Together, these findings demonstrate that our variational inference framework enables efficient and reliable recovery of both hyperparameters and decoding weights across diverse settings.
18

---

Published as a conference paper at ICLR 2026
ELBO vs iterations
B
la (length-scales)
(variances)
-10*
HT 35% 12.19%
cnr 56% 81.%
ELBO
inferred
-10
0
2
iteration
ground truth
Ground truth vs inferred decoding weights
6 0.1
0224
0.38
Egs 0.49
4yss 1.32
3.7%
3.0%
3.7%
12.2%
2r
y
ground truth
inferred
Figure S1: Simulation study results. (A) ELBO versus training iterations. (B) Scatter plots comparing ground truth versus inferred length scales (left) and variances (right) across features; insets report mean and standard deviation normalized absolute errors. (C) Decoding weights for 5 representative features: ground truth (black) vs. inferred weights (red), shaded regions show posterior standard deviation. True length scales are shown above each subplot; normalized mean absolute errors are shown below. Full results in Fig. $2
19

---

Published as a conference paper at ICLR 2026
A.2 COMPREHENSIVE VISUALIZATION OF TRUE INFERRED DECODING WEIGHTS (EXTENSION OF FIG. ST)
To provide complete view of model performance, present decoding weight recovery for all 200 simulated features used in the simulation study in Section A. This figure extends the subset shown in Fig. [T and illustrates how well the inferred weights (red) match the ground truth (black) across all features. Features with varying ground truth GP length scales are included, demonstraing that our inference method successfully recovers both broad and sharply tuned weight functions.
www
Figure S2: Full decoding weight recovery for all 200 simulated features, extending the results in Fig. SI Each subplot shows the ground truth decoding weights (black) and the inferred weights (red) across the continuous output space for one feature, with shaded regions showing posterior standard deviation. The true GP length scale is indicated in parentheses. The results demonstrate that the proposed inference procedure accurately captures both smooth and localized tuning profiles.
20

---

Published as a conference paper at ICLR 2026
A.3 ROBUSTNESS AND SCALABILITY OF CMLR INFERENCE ACROSS SIMULATION SETTINGS (EXTENSION OF SECTION A.T)
To assess the robustness and computational efficiency of the CMLR model, we systematically evaluated its performance across a range of controlled simulation scenarios, extending the results in Section A.1 Fig. S3 summarizes how inference accuracy and training time vary with the number of input features (D), number of observed samples (N), and number of Fourier components (M). These results extend the main text in Section and the simulation results in Fig. [yTby confirming that CMLR maintains high accuracy in hyperparameter and weight recovery while scaling gracefully with data size and model complexity.
Simulation performance comparison under different settings
A
Input Feauture (D)
B Observed samples (N)
C Fourier Components (M)
etT. in Pd
0-
CI.. in ld
0.5
CIT. mu(())
0.4
0.2-
Training time (s)
1000-
1000
2000
0
5000
10000
50
100
150
No. of input features (D)
No. of samples (N) No. of Fourier compoents (M)
Figure S3: Simulation performance comparison under varying conditions. Rows (top to bottom)
show: mean + standard deviation of the normalized absolute error in inferred variances (pa). length scales (L,), and decoding weights (w,(y)), followed by training time in seconds (Intel i7 2.4GHz CPU, 16GB RAM). (A) Varying the number of input features D shows that inference accuracy remains stable across all metrics, with only modest increases in training time. (B) Increasing the number of observed samples N substantially improves inference accuracy for all parameters, reflecting the benefit of additional data for GP hyperparameter and weight recovery. Notably, training time increases only marginally due to the use of stochastic variational inference. (C) Varying the number of Fourier components M reveals that performance improves with increasing M, particularly in recovering fine- scale structure in decoding weights. However, the computational cost does not increase substantially the frequency-domain formulation. Notably, near-asymptotic performance
due to the efficiency achieved with as few as M z 20 components. In practice, keeping M small is advisable to reduce dimensionality and mitigate the risk of overfitting. These results demonstrate that CMLR offers reliable and scalable inference across wide range of practical settings, with strong performance attainable using moderate data sizes and dimeinnnsi frequency representations.
21

---

Published as a conference paper at ICLR 2026
A.4 VALIDATION OF INFERENCE ACCURACY THE TWO-DIMENSIONAL SETTING (EXTENSION OF SECTION 3.4)
To confirm the accuracy and stability of CMLR in two-dimensional decoding tasks described in Section 3.4 we conducted simulation study with synthetic neurons tuned to 2D outputs. Fig. [S4] illustrates key results: convergence behavior of the variational inference algorithm, recovery of GP hyperparameters, and comparison of inferred vs. true decoding weights. Together, these analyses validate the model's s ability recover both the coarse and fine-grat ane-grained structure of 2D tuning profiles, demonstrating the reliability of CMLR in multidimensional output spaces.
A
FLBOO vs iterations
B & (lentth-saale dim 1)
(length-scale dim 2)
Pd Pdvvariance)
16% 32%
ELBO
inferred
40
20
2000
20
40
20
iteration
ground Iruth
Ground truth vs inferred decoding weights for representative neurons
7.0%
decode weights
wic 5.0%%
WTT 6.2%
4XX 4.7%
CTX. 4.5%
WY: 1 6.1%
-2 true weights
Figure S4: Simulation study validating inference accuracy in the 2D setting. (A) Evidence lower bound (ELBO) over training iterations, showing smooth and stable convergence of the variational inference algorithm. (completed in 469 seconds on a 2.4GHz Intel i7 CPU with 16GB RAM) (B) Scatter plots comparing ground truth and inferred GP hyperparameters (length scales and variances) across both output dimensions. Each point corresponds to one of the proximity the identity line indicates accurate recovery. (C) Comparison of inferred and ground truth decoding weights for subset of representative neurons. Each subplot shows a scatter plot of the true versus recovered weights, with the normalized mean absolute error shown inset. Points closely align with the identity line, indicating accurate recovery.
100 simulated neurons;
22

---

Published as a conference paper at ICLR 2026
A.5 DESIGN PARAMETERS AND PRACTICAL SETTINGS
This appendix provides practical guidance for selecting the design parameters of CMLR and outlines the specific settings used across all real-data experiments. A key advantage of CMLR is that it does not require extensive dataset-specific hyperparameter tuning. The main design parameters control numerical precision and computational efficiency rather than model capacity. These include:
1. the number of bins T used in the Riemann approximation of the normalization constant,
2. the number of Fourier basis components MM used to represent the GP prior,
3. the mini-batch size N' for stochastic optimization,
4. the learning rate a of the Adam optimizer, and
5. the decoding grid resolution J.
Table summarizes the design parameter settings used in each real dataset, along with the number of neurons (D), number of trials (N), and output range (S). These settings were fixed per dataset, but cross-validation can be incorporated into CMLR without modification. In practice, the two most important design parameters for numerical accuracy are the number Riemann bins and the number of Fourier components M.
determines the numerical resolution used to approximate the normalization constant and should match the desired output-grid resolution. In practice, T must be large enough that the density changes smoothly between adjacent grid points. For bounded or circular outputs, should scale with the width of the output range avoid discretization artifacts. As rule of thumb, for an output range of [0, 1]. choosing T Rs 100 provides adequate resolution, with larger values offering limited additional benefit.
Choice of T (Riemann bins):
Choice of M (Fourier components): M determines the spectral resolution the Gaussian process prior. For smooth kernels such the RBF, most of the spectral mass lies in low-frequency components, so only a small number of terms is required. Empirically, M E [15, 50] captures the relevant structure across all datasets while keeping computational cost low. As shown in Fig. S3 decoding performance is largely insensitive to M beyond this range. Choosing a small M is also desirable for reducing dimensionality and mitigating overfitting.
Choice of N' (mini-batch size): The mini-batch size affects optimization speed rather than model capacity. Smaller batches increase the variance of stochastic gradients (which can aid exploration), whereas larger batches reduce noise but increase memory cost. We selected N' to balance runtime and CPU memory constraints, with values between 1000 and 2000 working robustly.
Choice of (learning rate): The learning rate determines the step size in stochastic variational optimization. Unlike and M, which control numerical resolution, a governs the stability and convergence speed of the optimizer. In practice, moderately small value (e.g., E [10-3, 10 10-')) provides stable convergence across datasets. Larger values can lead to divergence, whereas smaller values slow training without improving accuracy. We selected for each dataset using small number preliminary runs and found that its performance was largely insensitive within this range.
Choice of J (decoding grid resolution): Finally, the number of decoding bins J must be specified when evaluating predictive densities from the CMLR model, as is also the case for FlexCode and Naive Bayes. Since we observed consistent trend of improved performance with finer decoding grids, we report all final results at the highest resolution tested, J = 20, 000.
Table 2: Design parameter settings used in the real-data studies
Dataset
D (neurons)
N (trials)
(range)
T
M
N
Mouse V1
11311-20616
4282-4469
[0, 2T)
100
17
1500
0.005
Monkey VI
113-148
3600
0, 2T))
100
1500
0.005
Mouse CA1
151-497
3600-5524
0,
100
20
1400
0.02
Monkey motor cortex
164
25299
-30,30
1000
27
2000
0.1
23

---

Published as a conference paper at ICLR 2026
A.6 SCALABILITY OF CMLR WITH RESPECT TO NEURON COUNT AND DATASET SIZE
Here, we assess the robustness of CMLR on the mouse VI dataset (data from Stringer et al. 2022 under systematic variation of the neuron count (D) and dataset size (N), comparing scaling behavior to XGBoost and DNN, Fig. SS characterizes how decoding accuracy and training time scale under these manipulations. As the number of neurons increases (Fig. (5(A), CMLR maintains stable absolute circular error and exhibits training times that grow approximately linearly with D, consistent with our simulation results in Appendix A..3 Across all values of D. CMLR achieves lower decoding error than XGBoost and DNN, and the performance gap widens at smaller D, indicating improved data efficiency. When the number of training samples is reduced (Fig. (SB), CMLR's performance degrades only slightly, whereas the data-driven baselines show more substantial deterioration. Training time for CMLR increases sublinearly with N due to mini-bate stochastic variational inference. Overall, CMLR shows favorable computational and statistical scaling, and critically, offers clear advantages over data-driven models in low-data regimes where principled structure and regularization matter most.
A
Neurons (D)
B
Training samples (N)
mean abs circular error
60
100
CMLR
XGBoost
DNN
0-
0-
training time (s)
8000
1000
0
0
102
100
1000
3000
no. of neurons (D)
no. of training samples (N)
Figure S5: Scalability of decoding performance and training time for CMLR, XGBoost, and DNN on the mouse V1 dataset. Top: mean + standard deviation of the absolute circular decoding error. Bottom: training time (seconds) measured on a standard laptop (Intel 2.4 GHz CPU. 16 GB RAM). Both metrics are shown function of (A) the number of neurons D and (B) the number of training samples N.
24

---

Published as a conference paper at ICLR 2026
A.7 MACAQUE VI DATA: DECODING DRIFTING GRATING ORIENTATIONS (EXTENSION OF SECTION 5.1)
Here, show the results of applying our method to five electrophysiological recording datasets from the monkey primary visual cortex (data from Graf et al. (20TT)) extending the results in Section 5.1 These datasets contained spiking activity recorded from between D 113 and D 148 neurons (treated as imput features), and included 72 discrete stimulus orientations (outputs) spaced 5 intervals, with 50 trials per orientation, resulting in a total of N 3600 samples per dataset. stimuli, the discrete case is a special instance
Although our method is designed for continuous-valle of the framework and fully supported without modification. We applied the CMLR model using the same inference settings in the mouse VI dataset.
Fig. S6 summarizes the results. panel A, CMLR and FlexCode consistently outperform Naive Bayes mean absolute circular enror across decoding resolutions, with all models plateauing beyond 500. Fig. SS6 shows true versus decoded orientations for all models; CMLR predictions align closely with the identity line, and most large errors occur near 180%, reflecting underlying bimodality, Panel C shows that FlexCode achieves the lowest decoding error (mean + SD: 13 36 median: 3", IQR: 5"), closely followed by CMLR (14 + 42", 3 49), with Naive Bayes performing moderately worse (21 + 50%, 4", 6"), and XGBoost (27 + 29", 17 28"), and DNN (26 + 26", 20", 26") showing higher variability and error.
These results confirm the robustness of CMLR across both continuous and discretized decoding settings. Consistent with findings in the mouse VI datasets, CMLR's correlation-awee -waare formulation yields improved accuracy over correlation-blind baselines, reinforcing the importance of modeling shared variability in cortical populations for reliable neural decoding.
A Error decoding resolution
B
True decoded
C
Error comparisons
CMLR
mean abs cireular error
25-
Naive Bayes (NB)
FlexCode
decoded orientation
300
abs citcular ciror (9)
100
36'
27
1
2
15.
250
true orientation (?)
CMLR FlexCode
NB
XGBoost
DNN
10a
10P
10
(number of decoding classes)
CMLR
FlexCode NB
XGBoost
DNN Ideal
Figure S6: Application of the proposed method macaque V1 data (data from Graf al. (2011) (A) Mean absolute circular error for CMLR, FlexCode, and Naive Bayes models as a function of the number of decoding classes J. (B) Scatter plots of decoded versus true stimulus orientations for
CMLR, FlexCode, Naive Bayes, XGBoost, and DNN models. (C) Box plots comparing the absolute
circular decoding error across all methods. Mean errors (#t standard deviation) indicated as insets.
25

---

Published as a conference paper at ICLR 2026
A.8 FULL DECODING WEIGHT MAPS FOR MONKEY MOTOR CORTEX DATA (EXTENSION OF FIG.
To provide complete visualization of the model's output, show the full set of 2D decoding weight functions inferred by CMLR for all 164 neurons in the monkey motor cortex dataset illustrated in Section 5.3 This Figure extends the subset shown in Fig. 4 revealing the rich diversity 2D velocity tuning profiles captured by the model. The smooth, continuous weight surfaces illustrate how CMLR uncovers interpretable tuning structure across the two-dimensional velocity space, consistent with known motor cortex encoding properties.
Figure S7: Full set of decoding weights inferred by CMLR for all D 164 neurons in the monkey motor cortex dataset (data from Glaser et al. (2018) 2020)), extending the results in Fig, 4 Each subplot shows the neuron's 's two-dimensional tuning surface across the and y velocity space. These results highlight the diversity and smoothness of velocity tuning captured by the model in real neural recordings.
26

---

Published as a conference paper at ICLR 2026
A.9 INTERPRETABILITY AND GENERALIZATION OF CMLR DECODING WEIGHT FUNCTIONS
In this Appendix, we further probe the CMLR decoding weight functions inferred for the motor-cortex velocity decoding task (data from Glaser et al. (2018 2020)) shown in Fig. 4 Our goal is to evaluate both the interpretability of the inferred weights and their ability to generalize beyond the range of velocities observed during training. Fig. S8 presents two complementary analyses that address these questions by comparing the inferred weight functions to empirical neural tuning structure and by examining how the inferred functions behave when trained on progressively restricted subsets of the data.
First, we compared the inferred weight functions wa(y) with empirical spike derived directly from the training data (Fig- (S8A). To obtain these empirical maps, we binned all training samples according to their two-dimensional velocity and, for each neuron d, accumulated that neuron's 's spike counts within each velocity bin. The resulting smoothed and normalized to form empirical firing-rate density maps over the velocity space. These maps reveal the raw statistical structure of each neuron's 's tuning. They are strongly concentrated within the region of velocities that were actually visited during the experiment, which is typically a compact region around the origin, and they are often noisy and irregular because of finite sampling. In contrast, the CMLR inferred decoding weights exhibit smooth and spatially coherent tuning patterns defined over the entire velocity space. Importantly, in regions where empirical data are dense, the inferred weights align well with the empirical structure, capturing the dominant preferred directions and suppressive regions observed in the raw densities. This agreement indicates that the CMLR model faithfully extracts meaningful tuning characteristics from the data while simultaneously regularizing them into smooth functional form that reflects the underlying relationship between velocity and firing activity.
spike-velocity velocity associations
weighted histograms were then
Second, we studied how the inferred decoding weights behave when the training data are progressively restricted to smaller velocity ranges (Fig. S8 By training the model using only samples within |/ |<| 25, 20, 15, 10, we artificially reduce the support of the observable data and thus increase the amount of extrapolation required outside the training domain. As the black boxes in Fig. B8B indicate, large portions of the velocity space are never observed when these restrictions are applied. Despite this, CMLR continues to produce well-structured and neuron-specific weight patterns that extend smoothly into unobserved regions. The fine details of the tuning curves become more dependent on the Gaussian process prior as the available data shrink, but the global structure remains consistent across all levels of data restriction. Preferred directionality, antagonistic regions, and smooth gradations remain stable as long as the model has at least partial coverage of the tuning landscape. This behavior reflects the inductive bias imposed by the smooth Gaussian process prior. When data are sparse or absent, the model defaults to the simplest and smoothest continuation that remains consistent with the observed samples.
Together, these analyses demonstrate two key properties of the CMLR framework:
1. Interpretability; the inferred decoding weights correspond closely to empirical neural firing statistics where data are available, offering intuitive and biologically meaningful descriptions of each neuron's 's response characteristics.
2. Generalization: the weight functions extend gracefully beyond the data supporedd regions, providing coherent and stable extrapolations driven by the model's smoothness prior rather than fitting noise or artifacts.
These findings show that CMLR provides not only accurate decoding performance but also structured and interpretable functional maps that remain reliable when parts of the stimulus space are sparsely sampled. This is important for real-world neural decoding settings. because behavioral or sensory spaces are rarely sampled uniformly in practice.
27

---

Published as a conference paper at ICLR 2026
Comparison of the inferred weights with the empirical densities
neuron 22
acuron 75
neuron 90
pcuron 94
ncuron 97
neuron 112
neuron 147
Inferred
weights
2
Empirical density
velocity
20-
0.004
>-200
velocity
B
Generalization of the inferred weights beyond the training set
neuron 22
neuron 75
neuron 90
neuron 94 neuron 97
neuron 112 neuron 147
all data
Iel 1y0%
25
[10% byt < 20
(ati 1y0%
< 15
velocity
20
[1% <
-20
25
X velocity
Figure S8: Interpretability and generalization of CMLR decoding weight functions for motor-cortex velocity decoding (data from Glaser al. (2018) 2020 (A) Inferred decoding weight functions for example neurons (top) compared with their empirical spike-weighted density maps (bottom). The empirical densities concentrate within the observed velocity range, while CMLR produces smooth, interpretable functions over the full space. (B) Generalization of inferred weights when training data are progressively restricted smaller velocity bounds. Top bottom: unconstrained. |, (g (2) 25, 20, 15. 10. Black boxes indicate the region of observed data. Even with limited training coverage, CMLR extrapolates smoothly outside the data-supported domain, demonstrating robustness and interpretability.
28

---

Published as a conference paper at ICLR 2026
A.10 CALIBRATION OF PREDICTIVE POSTERIORS
Here, we present the results of calibrating the posterior densities predicted by the CMLR model and FlexCode model held-out data. While the conditional density estimation methods trained to optimize the conditional likelihood, this does not necessarily guarantee that their predicted posteriors are well-calibrated; that is, that their uncertainty reflects empirical variability in the data. To assess calibration, we computed Probability Integral Transform (PIT) values (Dawid 1984) and quantile calibration curves (Kuleshov et al.) 2018), which diagnose whether predicted posterior distributions reflect the empirical distribution of true outputs. We also report the Expected Calibration Error (ECE) as a summary statistic.
A.10.1 PROBABILITY INTEGRAL TRANSFORM (PIT)
The Probability Integral Transform (PIT) provides diagnostic tool for assessing the calibration of continuous predictive distributions. For each test input Xn, let Fz (y) denote the predicted cumulative distribution function (CDF) for the output y. Given the true output Yn the PIT value is defined as:
Un Fn (vn)
which represents the cumulative probability mass assigned to values less than equal to yn under the predicted posterior. If the model is perfectly calibrated, the PIT values { {u, independently and identically distributed as Uniform(0, 1).
NN should be
Deviations from uniformity reveal miscalibration: U-shaped PIT histograms indicate overconfident predictions (posterior too narrow), hump-shaped histograms suggest underconfident predictions (posterior too wide), and asymmetric shapes reflect bias in the predictive distributions. Visual inspection of PIT histograms thus provides an interpretable diagnostic of calibration performance.
A.10.2 QUANTILE CALIBRATION CURVES
Quantile calibration curves assess whether predicted quantiles contain the correct proportion of ground truth outputs. For each test input x,, and each nominal quantile level a E (0, 1), we compute the a-quantile g (a) of the predicted posterior and then determine the empirical frequency with which the true output falls below this value. The empirical coverage at level a is defined as:
1 N
Coverage(a) (a)
96
where N is the number of test samples and TT (.) is the indicator function. Plotting empirical cover- age versus nominal yields the quantile calibration curve, with the diagonal representing perfect calibration.
To summarize calibration error across all quantile levels, we compute the expected calibration error (ECE) as:
1 Σ
ECE-
Coverage(a) - -al,
|A|
aEA
where A is finite set of quantile levels (e.g., (0.05, 0.10...... ..., 0.95)). A smaller ECE indicates better calibration.
A. 10.3 CALIBRATION ANALYSIS: CMLR vs FLEXCODE
We compared the uncertainty calibration of CMLR and FlexCode across three datasets (Mouse V1, Macaque VI, and Mouse CA1) using three diagnostics: PIT histograms, quantile calibration curves, and the relationship between PIT values and decoding error (Fig. S9) We restrict calibration analysis to CMLR and FlexCode because calibration diagnostics such as PIT histograms and quantile calibration curves require access to the full conditional density. Regression- models like XGBoost and DNN provide only point predictions rather than probability distributions, so PIT values, quantile calibration, and ECE cannot be computed for them.
PIT histograms (Fig, A))): PIT histograms for CMLR are close to uniform across datasets, in- dicating that the model allocates probability mass in way that matches the empirical distribution
29

---

Published as a conference paper at ICLR 2026
of the observed outcomes. This reflects well-calibrated posterior predictions with systematic concentration of mass in any region of the cumulative distribution. FlexCode, in contrast, shows clear deviations from uniformity, often exhibiting central peak or bimodal shape. These pattems indicate that FlexCode distributes probability mass unevenly and is therefore miscalibrated, assigning much weight to certain regions of the output space relative to the true distribution.
Quantile calibration curves (Fig, 59B): CMLR's calibration curves closely follow the identity line across all datasets, with low ECE values (Mouse VI: 0.02 + 0.02; Macaque VI: 0.03 + 0.01: CAI: 0.04 + 0.01), indicating accurate uncertainty quantification throughout the full quantile range. FlexCode shows larger deviations from the diagonal, with higher ECE values (0.09 0.03 0.05 = 0.01), consistent with its PIT histograms. Its empirical coverage typically falls below the ideal line at small quantiles and rises above it at large quantiles, indicating systematic miscalibration that produces under-coverage in the lower tail and over-coverage in the upper tail.
PIT versus decoding error (Fig. S9C): A well-calibrated model should show a clear relationship between PIT values and decoding error. In an ideal case, errors should be smallest near PIT R 0.5, where the model is most confident, and should increase smoothly as PIT approaches 0 or where the model expresses greater uncertainty. This pattern indicates that the model's predictive distribution correctly reflects which points are easy or difficult to decode. Across all datasets, CMLR closely matches this ideal behavior: errors are lowest near PIT x 0.5 and rise gradually toward the extremes, demonstrating that the model's uncertainty estimates align with actual decoding difficulty. FlexCode, in contrast, shows elevated errors near PIT 0.5, which is also where its PIT values are most concentrated, along with irregular fluctuations across the quantile range. This pattern indicates that FlexCode' 's predictive uncertainty does reliably correspond true error and that its diffuse predictions do not effectively communicate when mistakes are more likely.
Summary: Across all diagnostics and datasets, CMLR provides substantially betier-calibrated and more informative uncertainty estimates than FlexCode. CMLR's posterior distributions align closely with empirical data and express uncertainty coherently, whereas FlexCode tends to produce overly diffuse densities that obscure predictive reliability. This demonstrates that CMLR not only improves decoding accuracy but also offers superior uncertainty quantification, which is essential for scientific and neural decoding applications.
30

---

Published as a conference paper at ICLR 2026
A
Probability Integral Trasform (PIT)
B
Quantile Calibration C
PIT decoding error
Curves
samples
mean idcal
median
Mouse VI
ECE 0.02 0.02
0
20
Macaque VI
CMLR
ECE 0.03 0.02
0
Mouse CAI
ECE 0.03 0.01
20
Mouse VI
BCE 0.09 + 0.01
20
FlexCode
Macaque VI
ECE 0.01
empirical coverage
Mouse CAI density
decoding error
0
ECE 0.04 0.01
PIT value
Nominal quantile (a)
PIT value
Figure S9: Calibration analysis for CMLR (top) and FlexCode (bottom) across Mouse V1, Macaque
VI, and Mouse hippocampus CAI datasets. (A) Probability integral transform (PIT) histograms. CMLR produces ncar-uniform histograms, whereas FlexCode exhibits clear non-ume uniformiimmmmo distributions. (B) Quantile calibration curves comparing empirical coverage nominal quantile o; the diagonal indicates perfect calibration. Insets report the expected calibration error (ECE; mean + standard deviation across folds and datasets). (C) PIT value versus decoding error for individual test points, with median trends shown in red to illustrate how uncertainty relates prediction emror
31

---

Published as a conference paper at ICLR 2026
A.11 DECLARATION OF USAGE OF LARGE LANGUAGE MODELS (LLMS)
We used large language models (LLMs) only for minor editing tasks, such as polishing grammar and improving readability. LLMs were not used to generate content, perform analyses, design methods, or conduct experiments. All scientific contributions and results in this manuscript are solely the work of the authors.
32


Reviews:
[{"review_id": "9j3VjaPfgB", "reviewer": "Reviewer_LDws", "summary": "This paper develops a novel approach, Continuous Multinomial Logistic Regression (CMLR), that is an extension of MLR to predict continuous-valued outputs, and showcase this approach for neural decoding (predicting behaviors/stimuli from neural activity). The model enables interpretable understanding of smooth 'tuning curves' of neurons within the decoding model via GP priors. They develop an approach to fit the model with stochastic variational inference. The authors demonstrate excellent neural decoding performance, in additional to interpretability of the model across multiple datasets using their method.", "questions": "1. Fig 4C - Euclidean error is hard to interpret - Coefficient of determination (the standard r2_score in sklearn) would be easier to interpret, and more standard for velocity decoding\n\n2. Discussion: \"First, unlike models such as XGBoost or DNNs, CMLR does not incorporate priors over outputs\" - I don't understand this statement in terms of XGBoost and DNNs having priors on outputs.\n\n3. Line 263 - closing bracket ] missing\n\n4. I would put the methods in section 5.3 into an actual methods section. They feel very out of place when reading results\n\n5. How long does your approach take to run compared to others?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "This is an original, creative, and novel technical development for neural decoding (and prediction problems more broadly). The method has a good blend of accuracy and interpretability (via the learned tuning curves within the model). The paper is clearly written and the authors are very thorough. It is a very strong paper overall.", "weaknesses": "My one fundamental concern is with the comparisons in results, specifically in terms of hyperparameter selection. I might have missed it, but it's not clear to me how/if hyperparameter tuning was done both for the author's model and for comparison models. It seems odd how poor many of the results are from what should be close to state-of-the-art approaches (XGBoost and DNNs), which makes me suspect poor hyperparameters (causing either overfitting or underfitting). Proper hyperparameter tuning on a held-out validation set (within the training set) should be done, if it wasn't already.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "7vaX8gNOXA", "reviewer": "Reviewer_R8W8", "summary": "This paper introduces CMLR, a generalization of classical multinomial logistic regression to handle continuous output variables. Instead of discrete class weights, CMLR defines smooth, output-dependent weight functions with GP priors, allowing it to model conditional probability densities over continuous variables such as time, orientation, or spatial position. The authors derive a SVI algorithm in the Fourier domain for scalable learning on large datasets. They show CMLR’s performance on diverse neural decoding tasks across multiple brain regions (mouse and monkey V1, hippocampus CA1, motor cortex), and show that it outperforms Naive Bayes, FlexCode, XGBoost, and DNN in both accuracy and calibration.", "questions": "Despite using SVI for scalability, the method may still be computationally expensive for large-scale or high-dimensional datasets. A theoretical or empirical runtime comparison with baselines would strengthen the paper.", "limitations": "", "rating": 4, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. Although it relies on strong modeling assumptions than nonlinear methods, the approach is flexible enough to handle circular and multidimensional outputs.\n\n2. The method is interpretable, learning weight functions that correspond to tuning curves.\n\n3. This method offers an attractive alternative for researchers who prefer to avoid the complexity and hyperparameter tuning required by nonlinear models while still getting strong decoding performance.", "weaknesses": "The DNN results in Figs 2 and 3 show relatively poor performance with large variance. While I understand that DNNs typically have higher variance than the proposed method, the degree of underperformance here suggests that the architecture or hyperparameter tuning might not have been well optimized for these decoding tasks. This raises some concern about the fairness of the comparison. It would be helpful for the authors to acknowledge this limitation in the paper. If, on the other hand, the proposed method achieves strong performance without requiring extensive hyperparameter tuning, that could be an additional advantage of the proposed method.\n\nRelatedly, I would encourage the authors to clarify how they see the contribution of this work in the context of increasingly expressive modern models and large-scale neuroscience datasets. Although the proposed method is more interpretable, it may lag behind nonlinear models (e.g., transformers) in decoding performance. From a practical standpoint, the utility of the method might be limited. I would appreciate a clarification of how the authors envision their approach complementing or coexisting with more complex models.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "zrUAl98iJp", "reviewer": "Reviewer_4VSC", "summary": "The authors introduce continuous multinomial logistic regression (CMLR), a flexible nonparametric model that allows for both discrete- and continuous-valued outputs by mapping inputs to a full probability density using per-feature additive functions with Gaussian priors. The resulting model is applied to a wide variety of neural decoding tasks, both continuous and discrete, and shows impressive results compared to baselines.", "questions": "What is the computational complexity of CMLR, i.e. how do training and inference time scale with D, T, M, J? The authors have included training times in the Appendix, but it might be useful to mention these numbers (or at least their order of magnitude) more explicitly in the main text.\n\nHow much data do I actually need to train CMLR? Fig S3 is very cool and helpful, and it would be interesting to see something similar with real data. Do the GP priors/additive structure of CMLR make it more amenable to fitting models with less data, as compared to XGBoost or DNNs?\n\nFig 2E/3C/4C/S4C: what is the value of J used for CMLR/FlexCode/NB? How was this value chosen?\n\nI find the brief mention of uncertainty calibration very interesting; it is well-known that DNNs, for example, are very often poorly calibrated. It would be cool to see the PIT results (i.e. Fig S5) for some of the other methods. Is it possible to do this for DNNs? Better calibration could be as strong a selling point as better accuracy for scientific applications, and might be worth emphasizing this in the main text more.\n\nMore with uncertainty calibration: having never seen PIT histograms, I can believe they are a good approach for quantifying calibration, but they feel very disconnected from the actual datasets that are being analyzed. The mouse hippocampal data offers an interesting example: here CMLR tends to make mistakes when the true or decoded position is at one end of the track. What do uncertainty estimates look like for these mistakes? Does the uncertainty scale with the magnitude of the error? How do those compare to uncertainty estimates from the other methods?\n\nMinor: not having a background in this literature, the following sentences in the first Intro paragraph were a bit confusing to me: \"However, many neural decoding tasks involve continuous variables... . In such settings, researchers who wish to use MLR-like decoding methods are commonly forced to discretize the output variable into a finite number of classes.\" At this point I don't really know what \"MLR-like decoding methods\" are and the first thing I think of is \"sure but doesn't linear regression work just fine?\" Perhaps clarifying this point early will help convince readers unfamiliar with the literature.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "The paper is well-written and clearly articulates both the gap in the current literature and the exposition of the method. As applied to the neural decoding problem, the weight functions become per-neuron tuning curves that provide more interpretability than tree-based or neural network approaches. The method is applied to a range of neural decoding problems (sensory and motor) and shows strong performance across all datasets. The limitations and future directions section makes clear that this is a rich model class with many exciting directions to explore, and the paper is a strong foundation from which to start this work.", "weaknesses": "My main concern is the thoroughness of the baseline comparisons. The authors chose a different number of Fourier components M and mini-batch size N' for each of their datasets, but did not clearly state how these values were chosen. Furthermore, it seems no attempt was made to search hyperparameter space for XGBoost or the DNN. Better hyperparameter tuning (especially for the baselines) would therefore give me more confidence in the stated results. The authors could, for example, choose 2 or 3 important hyperparameters per method and search over these by performing a train/val split of the 80% of training data and selecting the best hyperparameters on the validation data.\n\nI am also concerned about the robustness of the results if they are indeed only reported on 20% of the data. My own experience in neural decoding has been that the train/test split can significantly affect model performance, and conclusions from one split may not hold with another split. Performing full k-fold cross validation (where every trial lands in the test set exactly once) would better control for noise introduced in the sampling process.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "TV1Qq5077K", "reviewer": "Reviewer_Z69w", "summary": "The paper introduces a combination of multinomial logistic regression with the gaussian processes defined on the weights.\nThe method is designed for neural decoding, where the goal is to estimate external variables from neural population activity, such as running speed or direction.\nTo handle computational complexity, the method assumed a univariate gaussian process per weight, uses standard radial basis function in the Fourier basis and trains the model as a stochastic variational inference.\nThen they test the model on several neural datasets, comparing with XGBoost and DNNs.", "questions": "1. How is CMLR related to the following works [1-3]?\n2. Why the XGBoost and DNNs lines are missing in Fig 2c and Fig 3A?\n3. I might have missed it in the text but which loss function do you use to train the CMLR? Is it MSE loss everywhere?\n4. Why there is no comparison with the standard multinomial logistic regression? You do not really analyse uncertainties in the main text, and even for uncertainties methods like Laplace Redux [3] could be used to derive uncertainties for a non-bayesian methods.\n5. Have you tuned the hyperparameters of DNN and XGBoost per dataset? As your datasets had different sizes the ratio of data to parameters might be crucial for performance. How DNN parameters compared to the CMLR learn parameters?\n6. What are the additional interpretability benefits provided by CMLR? XGBoost can also give an interpretable weighted impact of each neuron on the target variable.\n\nReferences:  \n[1] Chan, Antoni B. \"Multivariate generalized gaussian process models.\" arXiv preprint arXiv:1311.0360 (2013).  \n[2] Payne, Richard D., et al. \"A conditional density estimation partition model using logistic Gaussian processes.\" Biometrika 107.1 (2020): 173-190.  \n[3] Daxberger, Erik, et al. \"Laplace redux-effortless bayesian deep learning.\" Advances in neural information processing systems 34 (2021): 20089-20103.  \n[3] Murray, Iain, David MacKay, and Ryan P. Adams. \"The Gaussian process density sampler.\" Advances in neural information processing systems 21 (2008).", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. **Clarity**. The paper is clearly written and suggests extensive (theoretical) comparison with prior work (*\"Connections to Prior Work\"* section)\n2. **Various datasets**. The papers tests their method on various datasets from different species and data modalities (calcium imaging and electrophysiology)\n3. **Provided implementation**. The code is attached in the supplementary material.\n4. **Strong baselines**. The baselines are chosen to cover both bayesian and performance-drive methods (like DNN), covering the field.", "weaknesses": "1. **Unclear computational scalability**. The paper claims to provide a *\"scalable framework\"*, however, the practical aspects of scalability are under-explored. What are the computational restrictions? How many additional parameters the method provides and how longer does it take to train compared to modelling regression with uncertainty or the other baselines provided (such as FlexCode or XGBoost)?\n2. **Baselines**.The comparison with standard multinomial logistic regression is missing.  For the other strong baselines, it is not clear how they were optimized per dataset if any optimisations were present.\n3. **Lack of ablations for design choices justifications**.  While authors acknowledge in limitations that multivariate Gaussian processes could be used and fixed Fourier-domain bases and RBF kernels could be replaces by adaptive basis functions - all of these are not conceptual limitations of the methods but rather a list of sometimes very straightforward technical improvements (like multivariate Gaussian processes) and could be done within the current submission.\n4. **Unclear interpretability gains**. See Q6 -  What are the additional interpretability benefits provided by CMLR? XGBoost can also give an interpretable weighted impact of each neuron on the target variable. Multinomial logistic regression with uncertainty can also be close in terms of interpretations.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)