REVIEW 5 major objections 4 minor 29 references
A probabilistic Koopman model with inclusion flags can be shrunk to the observables that actually feed the outputs, without hurting prediction accuracy.
Reviewed by Pith at T0; open to challenge. T0 means a machine referee read the full paper against a public rubric. the ladder, T0–T4 →
T0 review · deepseek-v4-flash
2026-08-03 11:26 UTC pith:CAFH5YFN
load-bearing objection A genuine method paper: VB spike-and-slab plus SCC graph pruning for Koopman dictionaries; worth refereeing, but the output-unaffected guarantee needs validation against true supports and the per-system epsilon tuning documented better. the 5 major comments →
Koopman Model Dimension Reduction via Variational Bayesian Inference and Graph Search
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
Core claim
The paper's central claim is that the posterior inclusion matrix Γ, after thresholding at ε, encodes a directed graph whose nodes are dictionary observables, and that keeping only the strongly connected components that contain or feed the output observables yields a reduced dictionary that leaves the output dynamics unaffected. The authors derive variational Bayesian updates for a spike-and-slab model in which each Koopman regression coefficient is multiplied by a Bernoulli inclusion flag, so that weight magnitude and inclusion in the model are separated. The reduced dictionary is then obtained by graph condensation and ancestor search on the thresholded inclusion graph. The paper argues tha
What carries the argument
A spike-and-slab prior over Koopman regression coefficients, inferred by variational Bayesian updates, produces an inclusion matrix Γ whose entries are posterior probabilities of each coefficient being active. Thresholding Γ at ε turns it into the adjacency matrix of a directed graph over observables; strongly connected component condensation and ancestor search then identify the minimal dictionary needed to predict the outputs. This graph interpretation is what carries the reduction argument.
Load-bearing premise
The guarantee that outputs remain unaffected depends on the thresholded inclusion matrix being a complete map of every regression coefficient that can influence an output; if a coefficient that genuinely matters has posterior inclusion below epsilon, it gets silently pruned and the output guarantee fails.
What would settle it
Build a synthetic linear system with a known sparse dictionary where a small genuine coefficient connects a state to the output. Run the proposed VB inference with a modest threshold, then check whether that coefficient's posterior inclusion probability falls below epsilon. If it does and the reduced dictionary excludes that state, one-step prediction will degrade on test data, refuting the 'outputs unaffected' claim.
If this is right
- Overcomplete Koopman dictionaries can be pruned automatically without manually re-selecting observables, reducing model complexity for downstream control design.
- Ill-conditioned Koopman matrices, which often plague pseudoinverse and thresholded least-squares fits, become numerically stable after reduction, as shown by the Wiener–Hammerstein heatmaps.
- The reduced dictionary preserves one-step prediction NMSE for the systems tested, so the smaller model is not trading away predictive content.
- The reduction works as an aiding module: other identification methods benefit from the reduced dictionary even though the inclusion matrix was produced by the proposed Bayesian model.
- Threshold choice has little effect on the USV results, suggesting the inclusion probabilities converge to near-zero or near-one values rather than remaining indeterminate.
Where Pith is reading between the lines
- If the inclusion graph faithfully reflects true causal dependencies, the same graph-based pruning could be applied to other dictionary-learning or encoder-decoder Koopman pipelines, not just the authors' VB model.
- The method's reliance on a fixed threshold ε is a natural place to add uncertainty quantification; a data-driven ε or a threshold based on posterior credible intervals could make the reduction more robust.
- On systems with hidden or weakly excited modes, a coefficient with genuine but small influence could fall below ε and be silently dropped, so the 'outputs unaffected' guarantee is only as strong as the support-graph fidelity.
- A natural next test is a high-dimensional system with known sparse structure, where the reduced dictionary can be compared against the ground-truth generative support to measure how often true edges are discarded.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper proposes a hierarchical Bayesian spike-and-slab model for identifying finite-dimensional Koopman operator approximations. Each column of the EDMD regression is treated independently; the model includes inclusion flags gamma, coefficients beta, noise precision rho, and their hyperparameters, estimated by mean-field variational Bayes. The posterior inclusion matrix Gamma is thresholded at eps to produce a binary adjacency matrix over dictionary observables. Algorithm 2 condenses the resulting directed graph into SCCs and retains the SCCs that contain the output observables plus all their ancestors; the reduced dictionary is then used to re-estimate the Koopman model. Experiments on Lorenz, an underactuated USV, and the Wiener-Hammerstein benchmark compare the full and reduced dictionaries for four identification methods. The reduced dictionary improves conditioning of pseudoinverse and thresholded least-squares baselines and often preserves one-step NMSE; the proposed VB method itself matches SBL on the Lorenz and USV examples but underperforms SBL on the real Wiener-Hammerstein dataset.
Significance. The central idea using posterior inclusion probabilities to build a reachability graph over observables and reduce the dictionary while preserving output-reachable structure is interesting and potentially useful for Koopman model identification and control. The VB update derivations are transparent and the algorithms are specified in enough detail to be reimplemented. However, the main contribution as stated ('output variables remain unaffected by the reduction') is not established: it holds only if the thresholded inclusion matrix is a faithful support graph, which is never validated, and the paper's re-estimation protocol can change outputs (Table 1 shows a degradation for the proposed method on Wiener-Hammerstein). With additional validation and more honest framing, this could be a practical tool for dictionary reduction in EDMD.
major comments (5)
- [§4.2, Eq. (34); contribution bullet in §1] The guarantee that the reduced dictionary leaves output variables unaffected depends on Gamma_eps being a faithful support graph: every coefficient that can influence an output must have posterior inclusion probability >= eps, and no discarded coefficient may affect the output through an unmodeled path. The paper never validates Gamma_eps against ground-truth supports (e.g., synthetic regressions with known sparse coefficients) or against the full K_F. The threshold eps is set ad hoc (0.01 Lorenz, 0.25 USV, 0.1 Wiener-Hammerstein), and Figures 2 and 7 show the reduced size depends on eps. A false exclusion would silently remove an ancestor and change the outputs. Please add a validation of the inferred support (recovery of known nonzeros) and a sensitivity analysis over eps, or explicitly frame the guarantee as conditional.
- [§5.3, Table 1; §4.2, Algorithm 2] After computing the reduced dictionary, the paper re-estimates the Koopman model from the reduced dictionary rather than deriving the reduced model by restricting the full K_F to the kept observables. Even with a perfect support graph, re-estimation can change the outputs. Table 1 directly demonstrates this: for the proposed method, the test one-step NMSE on the Wiener-Hammerstein system increases from 4.6e-3 (full) to 7.3e-3 (reduced). This contradicts the contribution bullet stating the outputs 'remain unaffected by the reduction.' To support the graph guarantee, report the outputs of the reduced model obtained by projecting the full K_F onto the kept dictionary (submatrix), and compare those with the full model; present the re-estimated results as practical performance, not as the guarantee.
- [§4.1, Eq. (17d)] The update for the noise precision rho uses b_bar = 1/2 ||t - Phi(gamma_hat * beta_hat)||^2 + b. Under the mean-field posterior, the exact expectation of ||t - Phi(gamma * beta)||^2 contains the variance of each product gamma_i beta_i; because q factorizes, this adds sum_i gamma_hat_i[(1 - gamma_hat_i) beta_hat_i^2 + sigma_i^2] ||phi_i||^2, which is not accounted for. The update is therefore an approximation, not the exact VB update promised in the contributions. Please either derive the exact expectation or explicitly denote the point-estimate approximation and discuss its effect on the inferred inclusion probabilities.
- [§5, Eq. (35)] The performance evaluation uses only the one-step prediction NMSE, while the paper itself cites [8] to note that one-step performance can be misleading for Koopman models. Since the stated motivation is model predictive control and long-horizon prediction, please report multi-step or long-horizon NMSE for the full and reduced models, especially for the Wiener-Hammerstein experiment where the one-step metric already degrades for the proposed method. Without this, the claim that the reduction 'preserves or improves prediction performance' is not supported for the regimes of interest.
- [§5.1–§5.3] The inclusion-probability threshold eps is hand-tuned per system (0.01, 0.25, 0.1) and there is no criterion for its selection. Figures 2 and 7 show that the reduced dictionary size and hence the identified model depend strongly on eps and on the noise level. The paper calls the thresholding 'scale-free,' but the actual choice remains a user-defined free parameter. Please provide a principled selection rule (cross-validation, stability selection, or similar) and use it consistently, or explicitly list eps as a tuning parameter of the method.
minor comments (4)
- [General] No code or data repository is provided, which limits reproducibility despite the detailed algorithmic description.
- [§3, §4.1] There is a notation inconsistency: in Section 3, gamma_j is a vector of length L+l, but in Eq. (33) and Algorithm 1, gamma_hat_i is used for the scalar posterior inclusion probability. Please disambiguate.
- [§1, Abstract] Minor typos and grammar: 'Koopman operator recently gained' in the abstract should be 'has recently gained'; 'overparameterized' and 'overfitting' are used inconsistently. Please proofread.
- [Figures 8–9] The heatmaps use a logarithmic color scale spanning many orders of magnitude; consider adding annotations or a reproducible colorbar to make the comparison of sparsity patterns easier.
Circularity Check
No significant circularity: the central reduction is evaluated against external benchmarks; the only self-reference is a non-load-bearing baseline, and the output-preservation guarantee is conditional on the inferred support graph rather than a fitted prediction.
full rationale
The paper's core derivation chain is: define a spike-and-slab regression (Sec. 3), derive VB updates (Sec. 4.1), threshold the posterior inclusion matrix to form Γ_ε, and run SCC/ancestor search to select a sub-dictionary (Algorithm 2). None of these steps assumes the conclusion: Γ_ε is estimated from data by the VB mean-field equations, and the reduced dictionary is then tested on unseen data against independent baselines (pseudoinverse/EDMD, sequential thresholded least squares, SBL [26]) on Lorenz, USV, and Wiener-Hammerstein benchmarks. The 'outputs remain unaffected' claim is a graph-reachability property of Algorithm 2: by construction, only nodes that are outputs or ancestors of outputs in Γ_ε are retained, so relative to that thresholded graph no discarded node can reach an output. This is a conditional/definitional guarantee, not evidence that Γ_ε correctly identifies all true dynamical influences; but conditionality is a correctness risk, not circularity. The threshold ε is tuned per system (0.01/0.25/0.1), but the reported NMSE, conditioning, and dictionary-size outcomes are measured, not fitted. The only self-citation, [26], is used as a comparison baseline and in a concluding recommendation; it is not invoked to justify the VB derivation or the graph algorithm. The paper explicitly acknowledges tuning and mean-field limitations in Sec. 6. Therefore, no step in the derivation is equivalent to its inputs by construction; score 2 reflects the minor non-load-bearing self-citation rather than circularity.
Axiom & Free-Parameter Ledger
free parameters (5)
- ε (inclusion probability threshold) =
0.01 (Lorenz), 0.25 (USV), 0.1 (Wiener-Hammerstein)
- p_d (damping coefficient) =
not reported
- prior hyperparameters a,b,c_i,d_i,e_i,f_i =
not reported
- δ (clipping constant) =
1e-8
- RBF kernel centers, exponents, frequencies =
not reported exactly
axioms (6)
- domain assumption Mean-field factorization q(z)=∏_i q_i(z_i) (Eq. 6) is used for all posterior inference.
- domain assumption Gaussian likelihood for target variables, Eq. (8).
- domain assumption Spike-and-slab prior: β ~ N(0,α^{-1}), γ ~ Bernoulli(π).
- domain assumption Finite dictionary spans a near-invariant Koopman subspace (model (4) is valid).
- domain assumption Thresholded adjacency Γ_ε captures all true dependencies affecting outputs.
- standard math Strongly-connected-component condensation preserves reachability in directed graphs (Section 2.3).
read the original abstract
Koopman operator recently gained increasing attention in the control systems community for its abilities to bridge linear and nonlinear systems. Data driven Koopman operator approximations have established themselves as key enablers for system identification and model predictive control. Nonetheless, such methods commonly entail a preselected definition of states in the function space leading to high dimensional, overparameterized models that may suffer from poor numerical conditioning and degraded long term prediction performance. We address this problem by proposing a hierarchical probabilistic approach for the Koopman model identification problem. In our method, elements of the model are treated as random variables and the posterior estimates are found using variational Bayesian (VB) inference updates. Our model distinguishes from others in the integration of inclusion flags. By the help of the inclusion flags, we intuitively threshold the probability of each state in the model. We then propose a graph search based algorithm to reduce the preselected states of the Koopman model. We demonstrate that the proposed reduction improves numerical conditioning and can preserve or improve prediction performance while substantially reducing the dictionary size.
Figures
Reference graph
Works this paper leans on
-
[1]
L. Shi, M. Haseli, G. Mamakoukas, D. Bruder, I. Abraham, T. Murphey, J. Cortes, K. Kary- dis, Koopman operators in robot learning, arXiv preprint arXiv:2408.04200 (2024)
Pith/arXiv arXiv 2024
-
[2]
B. O. Koopman, Hamiltonian systems and transformation in hilbert space, Proceedings of theNationalAcademyofSciences17(5)(1931) 315–318
1931
-
[3]
Mauroy, Y
A. Mauroy, Y. Susuki, I. Mezic, Koopman operator in systems and control, Vol. 484, Springer, 2020
2020
-
[4]
P. J. Schmid, Dynamic mode decomposition of numerical and experimental data, Journal of fluid mechanics 656 (2010) 5–28
2010
-
[5]
M. O. Williams, I. G. Kevrekidis, C. W. Rowley, A data–driven approximation of the koopman operator: Extending dynamic mode decomposition, Journal of Nonlinear Science 25 (6) (2015) 1307–1346
2015
-
[6]
J. L. Proctor, S. L. Brunton, J. N. Kutz, Gen- eralizing Koopman theory to allow for inputs andcontrol, SIAMJournalonAppliedDynam- ical Systems 17 (1) (2018) 909–930
2018
-
[7]
I. Abraham, G. De La Torre, T. D. Murphey, Model-based control using Koopman opera- tors, arXiv preprint arXiv:1709.01568 (2017)
Pith/arXiv arXiv 2017
-
[8]
J. Pan, D. Li, J. Wang, P. Zhang, J. Shao, J. Yu, Autogeneration of mission-oriented robot controllers using bayesian-based koop- man operator, IEEE Transactions on Robotics 40 (2024) 903–918.doi:10.1109/TRO.2023. 3344033
-
[9]
M. Khosravi, Representer theorem for learn- ing Koopman operators, IEEE Transactions on Automatic Control 68 (5) (2023) 2995–3010. doi:10.1109/TAC.2023.3242325
arXiv 2023
-
[10]
Lusch, J
B. Lusch, J. N. Kutz, S. L. Brunton, Deep learning for universal linear embeddings of nonlinear dynamics, Nature Communications 9 (1) (2018) 4950
2018
-
[11]
Q. Li, F. Dietrich, E. M. Bollt, I. G. Kevrekidis, Extended dynamic mode decom- position with dictionary learning: A data- driven adaptive spectral decomposition of the Koopman operator, Chaos: An Interdisci- plinary Journal of Nonlinear Science 27 (10) (2017)
2017
-
[12]
N. C. Thompson, K. Greenewald, K. Lee, G. F. Manso, The computational limits of deep learning (2022).arXiv:2007.05558. URLhttps://arxiv.org/abs/2007.05558
Pith/arXiv arXiv 2022
-
[13]
D. G. Tzikas, A. C. Likas, N. P. Galatsanos, The variational approximation for Bayesian inference, IEEE Signal Processing Magazine 25 (6) (2008) 131–146.doi:10.1109/MSP. 2008.929620
arXiv 2008
-
[14]
S. L. Brunton, J. L. Proctor, J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the national academy of sciences 113 (15) (2016) 3932–3937
2016
-
[15]
S. Watanabe, Tree-structured parzen estima- tor: Understanding its algorithm components and their roles for better empirical perfor- mance (2025).arXiv:2304.11127. URLhttps://arxiv.org/abs/2304.11127
Pith/arXiv arXiv 2025
-
[16]
R. Nayek, R. Fuentes, K. Worden, E. Cross, On spike-and-slab priors for Bayesian equation discovery of nonlinear dynamical systems via sparse linear regression, Mechanical Systems and Signal Processing 161 (2021) 107986. doi:10.1016/j.ymssp.2021.107986. URLhttp://dx.doi.org/10.1016/j. ymssp.2021.107986
arXiv 2021
-
[17]
C. Schlosser, M. Korda, Sparse decompositions of nonlinear dynamical systems and appli- cations to moment-sum-of-squares relaxations (2024).arXiv:2012.05572. URLhttps://arxiv.org/abs/2012.05572
Pith/arXiv arXiv 2024
-
[18]
C. M. Bishop, N. M. Nasrabadi, Pattern recog- nition and machine learning, Vol. 4, Springer, 2006
2006
-
[19]
A. P. Dempster, N. M. Laird, D. B. Rubin, Maximum likelihood from incomplete data via the EM algorithm, Journal of the royal statis- tical society: series B (methodological) 39 (1) (1977) 1–22. 13
1977
-
[20]
M. I. Jordan, Z. Ghahramani, T. S. Jaakkola, L. K. Saul, An introduction to variational methods for graphical models, Machine learn- ing 37 (2) (1999) 183–233
1999
-
[21]
Nuutila, E
E. Nuutila, E. Soisalon-Soininen, On finding the strongly connected components in a di- rected graph, Information processing letters 49 (1) (1994) 9–14
1994
-
[22]
H. Ishwaran, J. S. Rao, Spike and slab variable selection: Frequentist and Bayesian strategies, The Annals of Statistics 33 (2) (2005) 730–773. doi:10.1214/009053604000001147
-
[23]
Saad, Iterative methods for sparse linear systems, SIAM, 2003
Y. Saad, Iterative methods for sparse linear systems, SIAM, 2003
2003
-
[24]
B. A. Vu, D. Gunawan, A. Zammit-Mangion, R-vgal: a sequential variational Bayes algo- rithm for generalised linear mixed models, Statistics and Computing 34 (3) (Apr. 2024). doi:10.1007/s11222-024-10422-8. URLhttp://dx.doi.org/10.1007/ s11222-024-10422-8
-
[25]
J. L. Spouge, Computation of the gamma, digamma, and trigamma functions, SIAM Journal on Numerical Analysis 31 (3) (1994) 931–944
1994
-
[26]
S. E. Özcan, M. M. Ankaralı, Sparse Bayesian learning for Koopman based system identifica- tion, in: 202511thInternationalConferenceon Control, Decision and Information Technolo- gies (CoDIT), IEEE, 2025
2025
-
[27]
E. N. Lorenz, Deterministic nonperiodic flow, Journal of the Atmospheric Sciences 20 (1963) 130–141. URLhttps://api.semanticscholar.org/ CorpusID:15359559
1963
- [28]
-
[29]
Schoukens, J
J. Schoukens, J. Suykens, L. Ljung, Wiener– Hammerstein benchmark, in: Proceedings of the 15th IFAC Symposium on System Identifi- cation (SYSID 2009), St. Malo, France, 2009, july 6–8, 2009. 14
2009
discussion (0)
Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.