Pith. sign in

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 →

arxiv 2601.06315 v2 pith:CAFH5YFN submitted 2026-01-09 eess.SY cs.SY

Koopman Model Dimension Reduction via Variational Bayesian Inference and Graph Search

classification eess.SY cs.SY MSC 93B3093C5562F1505C20
keywords Koopman operatorvariational Bayesian inferencespike-and-slab priordictionary reductiongraph condensationstrongly connected componentssystem identificationEDMD
verification ladder T0 review T1 audit T2 compute T3 formal T4 reserved

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper tries to establish that a Koopman model built from an overcomplete dictionary can be reduced to a much smaller dictionary without changing its output behavior. It does this by treating every regression coefficient as having an inclusion probability, estimated through variational Bayesian inference, and then using those probabilities as a directed graph. Nodes that cannot reach any output are removed. In three systems — a chaotic Lorenz attractor, an underactuated surface vehicle, and a real Wiener–Hammerstein circuit — the reduced models keep one-step prediction error at the same level while fixing ill-conditioned matrices. If true, this gives practitioners a principled way to prune dictionaries for EDMD-type models.

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.

Watch this falsifier — get emailed when new claim-graph text bears on it.

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

These are editorial extensions of the paper, not claims the author makes directly.

  • 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.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit.

Referee Report

5 major / 4 minor

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)
  1. [§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.
  2. [§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.
  3. [§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.
  4. [§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. [§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)
  1. [General] No code or data repository is provided, which limits reproducibility despite the detailed algorithmic description.
  2. [§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.
  3. [§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.
  4. [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

0 steps flagged

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

5 free parameters · 6 axioms · 0 invented entities

The contribution rests on standard Bayesian and graph-theoretic machinery; no physical entities are invented, but several hyperparameters (especially the threshold ε) are fitted by hand, and the central preservation guarantee depends on the support-graph assumption.

free parameters (5)
  • ε (inclusion probability threshold) = 0.01 (Lorenz), 0.25 (USV), 0.1 (Wiener-Hammerstein)
    Eq. (34): probabilities below ε are flagged 0. Chosen ad hoc per experiment; larger ε gives smaller dictionaries and can remove true edges.
  • p_d (damping coefficient) = not reported
    Eq. (24) damps β updates on rugged loss surfaces; no value given, so convergence behavior is not reproducible.
  • prior hyperparameters a,b,c_i,d_i,e_i,f_i = not reported
    Gamma/Beta priors, Eqs. (11)-(13), determine sparsity and inclusion behavior; paper says 'uninformative' but does not specify exact values.
  • δ (clipping constant) = 1e-8
    In Sec. 4.1.5, γ probabilities clipped to [δ,1-δ]; visible effect on tiny K̂_F entries in Figures 8-9.
  • RBF kernel centers, exponents, frequencies = not reported exactly
    Dictionary built in Sec. 5; exponents/frequencies set to 'different orders of magnitude' with no exact values, affecting every result.
axioms (6)
  • domain assumption Mean-field factorization q(z)=∏_i q_i(z_i) (Eq. 6) is used for all posterior inference.
    Used for all posterior updates; ignores dependencies between β, γ, α, π, ρ; acknowledged as a limitation in Section 6.
  • domain assumption Gaussian likelihood for target variables, Eq. (8).
    Measurement noise after observables is not necessarily Gaussian because dictionary functions are nonlinear; this is an unverified modeling assumption.
  • domain assumption Spike-and-slab prior: β ~ N(0,α^{-1}), γ ~ Bernoulli(π).
    A standard variable-selection prior, but its suitability for Koopman EDMD is assumed, not derived.
  • domain assumption Finite dictionary spans a near-invariant Koopman subspace (model (4) is valid).
    Standard EDMD assumption needed for the graph interpretation of K̂_F; otherwise the identified matrix does not represent the operator.
  • domain assumption Thresholded adjacency Γ_ε captures all true dependencies affecting outputs.
    Algorithm 2's preservation guarantee holds only if no discarded edge matters; never validated against ground-truth supports.
  • standard math Strongly-connected-component condensation preserves reachability in directed graphs (Section 2.3).
    Standard graph theory; not a physics or data assumption.

pith-pipeline@v1.3.0-alltime-deepseek · 13846 in / 15478 out tokens · 152579 ms · 2026-08-03T11:26:16.269586+00:00 · methodology

0 comments
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

Figures reproduced from arXiv: 2601.06315 by Mustafa Mert Ankarali, Selin Ezgi Ozcan.

Figure 1
Figure 1. Figure 1: The hierarchy of the random variables in the pro [PITH_FULL_IMAGE:figures/full_fig_p005_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: Reduced dictionary sizes of the Lorenz system with [PITH_FULL_IMAGE:figures/full_fig_p009_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: NMSE performance of full size models identified [PITH_FULL_IMAGE:figures/full_fig_p009_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: NMSE performance of reduced size models iden [PITH_FULL_IMAGE:figures/full_fig_p010_4.png] view at source ↗
Figure 6
Figure 6. Figure 6: NMSE performance of reduced size models iden [PITH_FULL_IMAGE:figures/full_fig_p011_6.png] view at source ↗
Figure 7
Figure 7. Figure 7: Reduced dictionary sizes of the USV with different [PITH_FULL_IMAGE:figures/full_fig_p011_7.png] view at source ↗
Figure 8
Figure 8. Figure 8: Heatmap of the recovered |Kˆ F | matrices of differ￾ent methods using the full sized dictionary. Gray indicates zero entries. weights, inclusion flags and resulting regression er￾ror are also treated randomly and inferred. For inference of the random variables, we use VB up￾dates. Finally, when the inclusion variables for all observables and in each observable regression prob￾lems are estimated, we get a m… view at source ↗

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

29 extracted references · 2 canonical work pages

  1. [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)

  2. [2]

    B. O. Koopman, Hamiltonian systems and transformation in hilbert space, Proceedings of theNationalAcademyofSciences17(5)(1931) 315–318

  3. [3]

    Mauroy, Y

    A. Mauroy, Y. Susuki, I. Mezic, Koopman operator in systems and control, Vol. 484, Springer, 2020

  4. [4]

    P. J. Schmid, Dynamic mode decomposition of numerical and experimental data, Journal of fluid mechanics 656 (2010) 5–28

  5. [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

  6. [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

  7. [7]

    Abraham, G

    I. Abraham, G. De La Torre, T. D. Murphey, Model-based control using Koopman opera- tors, arXiv preprint arXiv:1709.01568 (2017)

  8. [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. [9]

    Khosravi, Representer theorem for learn- ing Koopman operators, IEEE Transactions on Automatic Control 68 (5) (2023) 2995–3010

    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

  10. [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

  11. [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)

  12. [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

  13. [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

  14. [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

  15. [15]

    Watanabe, Tree-structured parzen estima- tor: Understanding its algorithm components and their roles for better empirical perfor- mance (2025).arXiv:2304.11127

    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

  16. [16]

    Nayek, R

    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

  17. [17]

    Schlosser, M

    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

  18. [18]

    C. M. Bishop, N. M. Nasrabadi, Pattern recog- nition and machine learning, Vol. 4, Springer, 2006

  19. [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

  20. [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

  21. [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

  22. [22]

    Ishwaran, J

    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. [23]

    Saad, Iterative methods for sparse linear systems, SIAM, 2003

    Y. Saad, Iterative methods for sparse linear systems, SIAM, 2003

  24. [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. [25]

    J. L. Spouge, Computation of the gamma, digamma, and trigamma functions, SIAM Journal on Numerical Analysis 31 (3) (1994) 931–944

  26. [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

  27. [27]

    E. N. Lorenz, Deterministic nonperiodic flow, Journal of the Atmospheric Sciences 20 (1963) 130–141. URLhttps://api.semanticscholar.org/ CorpusID:15359559

  28. [28]

    Atasoy, O

    S. Atasoy, O. K. Karagöz, M. M. Ankar- ali, Trajectory-free motion planning of an un- manned surface vehicle based on MPC and sparse neighborhood graph, IEEE Access 11 (2023) 47690–47700.doi:10.1109/ACCESS. 2023.3275433

  29. [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