{"id":"2f767c84-d368-4869-b4dc-1d6bc499950d","arxiv_id":"1908.03913","paper_version":2,"verdict":"ACCEPT","confidence":"HIGH","novelty_score":6.0,"correctness_risk":"low","formal_verification":"none","parameter_count":3,"one_line_summary":"Stable spline kernels induce RBF kernels for disturbance processes, enabling closed-form Gaussian process imputation of missing outputs in linear system identification.","lead":"The paper derives a new family of radial basis function kernels from stable spline kernels and uses them to impute missing output data in linear system identification. This gives a nonparametric replacement for parametric prediction error methods when some measurements are missing.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"No significant objection identified: the RBF kernel derivation, the discrete-time specializations, and the closed-form imputation are internally consistent, and the acknowledged Gaussian approximation is a non-fatal limitation.","rationale":"I read the paper as making two connected claims: (i) a new stationary covariance family for the disturbance is obtained from stable spline priors, and (ii) this enables closed-form linear minimum variance imputation. Claim (i) is a direct convolution calculation; I re-derived it and it is correct, including the discrete-time formulas. Claim (ii) is a standard kriging/Gaussian process regression result, and Proposition 3's algebra is consistent. The paper is honest about the Gaussian approximation in the marginal likelihood and does not claim exact non-Gaussian optimality. The residual risk is prior misspecification: if the true noise model is far from the stable spline prior, the covariance used in (24) is not the true one. This is a general property of any Bayesian or regularized estimator and does not undermine the central contribution, especially since the experiments randomize ARMAX orders and pole locations and still show strong performance. The reader's weakest_assumption points to exactly this misspecification/Gaussian-approximation risk, so I partially agree, but I do not consider it load-bearing enough to change the accept verdict.","tokens_in":12926,"tokens_out":29600,"duration_ms":295199,"concrete_test":"Run the ARMAX Monte Carlo of Sec. 6.1 with Student-t innovations (3-5 degrees of freedom) and compare SS imputation+SS against PEM+Oracle (missing) on COD_miss and COD_k. If the performance gap remains small, the Gaussian approximation in Eq. (26) is not load-bearing; if it degrades sharply, hyperparameter estimation via the Gaussian marginal likelihood becomes a concern that the current experiments (Gaussian e) do not address.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The load-bearing steps check out. Proposition 1 (Eq. 13) correctly computes the covariance of xi = e⊗f_p from the stable spline prior on f_p: for lag x, the integral or sum of K(y, y+|x|) follows directly from the whiteness of e and the independence of e and f_p. I verified the q=1 and q=2 closed forms, including the discrete-time formulas (18)-(19), using the corrected reading of Eq. (11) in which the second term is e^{-3β max}/6; the formulas match a direct summation. The enriched-kernel statement in Sec. 4.4 is consistent, since replacing K by the covariance of g_p⊗h in Eqs. 13/17 is exactly what the convolution covariance requires. Proposition 3's estimator (24) and posterior variance (25) are the standard linear minimum variance expressions, and the cross-covariance V(y(t), y_o) = λ(P(t,·)+R(t,·)) follows from independence of the input-driven part and the disturbance. The acknowledged Gaussian approximation in the marginal likelihood (Sec. 5.3, Eq. 26) is a real limitation, and prior misspecification of the noise covariance would remove the 'minimum variance' guarantee, but these are normal Bayesian modeling caveats rather than internal inconsistencies. The experiments, including the oracle-PEM baseline that selects order on test data, provide credible evidence for the competitive performance claim.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes a nonparametric Bayesian approach to identifying linear systems from output measurements with missing samples. The authors model the unknown impulse responses as zero-mean Gaussian processes with stable spline covariances and show that the disturbance component xi(t) = (e \\otimes f_p)(t) is stationary with a covariance that is a new RBF kernel obtained by integrating the stable spline kernel along its diagonal (Proposition 1, Eq. (13), with discrete-time versions in Eqs. (17)-(19)). Using this kernel together with the stable spline kernel, they derive a closed-form linear minimum variance estimator for missing outputs (Proposition 3, Eq. (24)) and an algorithm, stable spline imputation (Algorithm 4), in which hyperparameters are estimated by marginal likelihood under a Gaussian approximation (Eq. (26)). Numerical experiments on ARMAX models compare the method against PEM baselines, including an oracle that selects model order on the test target, and report competitive or better predictive performance even when the benchmarks use the full data set.","tokens_in":13227,"tokens_out":10593,"duration_ms":108219,"significance":"If the results are correct, the paper makes a useful contribution to kernel-based system identification. The main theoretical content, the derivation of RBF kernels from stable spline priors, is transparent and checkable: the q=1 case reduces to the Laplace kernel, and the q=2 closed form is consistent with direct summation. The estimator in Proposition 3 is the standard linear minimum variance expression, and the independence structure of the Bayesian network is clearly specified. The experimental design is strong: the PEM baselines include oracle order selection on the test target, which gives the parametric competitors an advantage, and the conclusion that the nonparametric imputation method is competitive is carefully worded. The paper is also clearly organized, and its scope covers both discrete- and continuous-time settings, which is broader than most existing missing-data solutions. The acknowledged Gaussian approximation in the marginal likelihood is a limitation but not an internal inconsistency. Overall this is a solid, publishable contribution.","major_comments":[],"minor_comments":[{"comment":"The notation for R_m is inconsistent: the text calls it a covariance matrix, but the displayed definition in Eq. (6) is E[\\xi_m \\xi^T], which is a cross-covariance, and Eq. (7) uses it in exactly that cross-covariance role. Please define R_{mo} = E[\\xi_m \\xi^T] explicitly and use a separate symbol for the auto-covariance of \\xi_m, or state clearly that R_m denotes the cross-covariance.","section":"Section 2.3, Eqs. (6)-(7)"},{"comment":"The proof of Proposition 1 is omitted with the note that it relies on simple computations. Since this proposition is the main theoretical novelty of the paper, please include a derivation in the appendix; this would also clarify the integration over R_+ and the role of |x|.","section":"Section 4.2, Proposition 1"},{"comment":"The Gaussian approximation for the disturbance is introduced but not discussed. Please add a sentence or two explaining the nature of the approximation and its expected effect on hyperparameter estimates; a brief empirical check would strengthen the paper, although the current experiments do not indicate a problem.","section":"Section 5.3, Eq. (26)"},{"comment":"The text says 'boxplots of the 1000 values of COD5' but the Monte Carlo study has 500 runs; this appears to be a typo and should be corrected to 500.","section":"Section 6.5, bottom panel"},{"comment":"The sentence reporting 'almost 87% of the cases' should clarify whether 'cases' refers to the 2000 individual COD_k values from the first 100 runs or to the number of Monte Carlo runs; the current wording is ambiguous.","section":"Section 6.4, top panel"},{"comment":"The definition of h_x(t) in the output kernel P(x,y) is hard to parse. Please rewrite it with a more explicit convolution formula, such as P(x,y) = \\sum_l \\int\\int u_l(y-s) u_l(x-r) K(s,r) ds dr, or state clearly that h_x(t) is the function t \\mapsto (u_l(\\cdot) \\otimes K(t,\\cdot))(x).","section":"Section 5.1, Eq. (21)"}],"recommendation":"minor_revision","confidential_remarks":"In my view the paper is essentially ready; the issues are local and can be fixed without further experiments. I would ask the authors to add the proof of Proposition 1 and to correct the notation in Eq. (6), as these would improve the paper's rigor and readability. No concerns about novelty or fit with the journal's scope."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Dear [Name],\n\nRead the Pillonetto–Chiuso–De Nicolao paper on stable spline imputation with missing data. The headline: this is a genuine extension of the stable spline framework to missing outputs, and the central mechanism—deriving a new RBF kernel family from stable spline priors—checks out. The math is clean; I verified the q=1 and q=2 discrete-time formulas in (18)–(19), and the cross-covariance argument behind the closed-form estimator (24) is standard and correct.\n\nWhat is new: the family of RBF kernels in (13)/(17), the discrete-time specializations, and the application to imputation via a minimum-variance linear estimator with hyperparameters estimated by marginal likelihood. The q=1 case reduces to the Laplace kernel, a useful sanity check, but the general family is not in prior work. The experiments are thorough: 500 Monte Carlo runs, missing-data PEM with oracle order selection, and full-data baselines. The claim that imputation plus a stable spline estimator is competitive with, and often better than, oracle PEM on missing data is credible—the oracle has genuine extra information, so the reported 87% figure is not a straw man.\n\nSoft spots, in proportion. The Gaussian approximation in the marginal likelihood (Eq. 26) is real and untested. It is a Bayesian modeling choice, not an error, but it means the hyperparameter estimate is not the exact marginal likelihood, and the minimum-variance guarantee is conditional on the prior being right. The MCAR assumption is implicit throughout; if missingness depends on the output, the method would need a different justification. No code or data are shipped, so the numerical results are not independently reproducible yet. The proof of Proposition 1 is omitted with a note that it uses simple computations—I re-derived it and it is fine, but a referee should ask for the details.\n\nNone of these are load-bearing. The central argument holds. This is a solid within-subfield contribution that people working on nonparametric system identification will want to read and cite. The self-citations are appropriate given that the stable spline framework is the authors' own. I would send this to peer review without hesitation, and would accept after minor revisions if the authors supply the missing proof details and comment on the Gaussian approximation.","headline":"Solid extension of stable spline identification to missing data; the kernel derivation is clean and the central claims hold, with a real but non-fatal caveat about the Gaussian approximation in the marginal likelihood.","tokens_in":13741,"tokens_out":2243,"would_cite":true,"duration_ms":24226,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["93E12"],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper establishes that stable spline kernels induce a new family of radial basis function kernels describing filtered white noise, so missing output data can be imputed in closed form and the resulting models can outperform full-data…","keywords":["linear system identification","missing data","Gaussian processes","kernel-based regularization","stable spline kernels","radial basis functions kernels","stable spline imputation"],"falsifier":"Run a Monte Carlo with an ARMAX system whose noise-path pole is much faster than the system's least-damped pole, so the shared-$\\beta$ RBF kernel is misspecified; if the closed-form imputer's mean-squared error on missing outputs is not close to that of the best linear estimator computed with the true noise covariance, the shared-$\\beta$ assumption is load-bearing. A simpler check is to estimate the residual autocovariance from a long simulation and compare its shape with Eq. (15).","tokens_in":12723,"feed_emoji":"🧩","tokens_out":5575,"duration_ms":53823,"temperature":0.7,"pith_summary":"This paper tries to establish a fully nonparametric cure for missing output data in linear system identification. Its key move is to show that the stable spline kernel, previously used to model impulse responses as Gaussian processes, also yields a new family of radial basis function kernels that describe the disturbance term of a linear model, the convolution of white noise with a stable impulse response. With both kernels in hand, the paper derives the minimum variance linear estimator of the missing outputs in closed form, parameterized only by a few hyperparameters estimated by marginal likelihood. Numerical experiments on ARMAX systems indicate that this stable spline imputation procedure is more predictive than parametric Prediction Error Methods, even when the latter see the full data set and use an oracle for model order selection.","feed_headline":"Stable spline imputation fills missing outputs, beats full-data PEM","feed_subtitle":"A two-parameter Bayesian kernel closes output gaps and predicts better than oracle-selected parametric models.","key_machinery":"The load-bearing object is the RBF kernel family defined by Eq. (13), $h(x) = \\int_0^\\infty K(y, y+|x|)\\,dy$, obtained by integrating the stable spline kernel $K$ along the diagonal; it expresses the second-order statistics of the noise component $\\xi$ as a function of time lag only. For $q=1$ it reduces to the Laplace kernel $h(x) = e^{-\\beta|x|}/\\beta$, and for $q=2$ to $h(x) = (3e^{-2\\beta|x|} - e^{-3\\beta|x|})/(18\\beta)$. This kernel, together with the stable spline output kernel $P$, enters the minimum variance linear estimator (24) and the marginal-likelihood criterion (26), which is optimized over only two or four hyperparameters before all missing outputs are returned in closed form.","core_discovery":"The central discovery is the covariance identity in Eq. (13): if the noise-path impulse response $f_p$ is a zero-mean Gaussian process with stable spline covariance $\\lambda K$, then the disturbance $\\xi(t) = (e \\otimes f_p)(t)$ is a zero-mean stationary process whose covariance is the translation-invariant radial basis function $R(s,t) = h(s-t)$, where $h(x) = \\int_0^\\infty K(y, y+|x|)\\,dy$. For the classic first- and second-order stable spline kernels this yields the explicit Laplace-like forms (14) and (15). Combined with the stable spline kernel for the system impulse responses, this makes the best linear predictor of the missing outputs in Eq. (24) computable in closed form; the paper calls the resulting procedure stable spline imputation. In Monte Carlo tests on ARMAX models with 25% output data missing, the imputation-then-identification pipeline produces models whose $k$-step-ahead prediction accuracy is better than that of PEM with oracle order selection in almost 87% of runs, and comparable to or better than full-data estimators.","pith_inferences":["The covariance identity suggests a direct diagnostic: given a long record from a known system, the empirical autocorrelation of the residuals after removing the deterministic input could be compared against the predicted $h(x)$; a mismatch would indicate either a wrong $\\beta$ or a noise path not well modeled by the stable spline prior.","The same closed-form imputer could be adapted to missing inputs rather than outputs by symmetry of the linear model, though the paper only treats output gaps.","Because the RBF kernel arises from a convolution, a similar construction may yield kernels for other Gaussian-process priors, such as Mat\\'ern or periodic kernels, extending the family of tractable imputation problems."],"forward_implications":["If the covariance identity holds, missing outputs can be imputed in closed form with no model-order search, just a two- or four-dimensional optimization for the hyperparameters.","The estimator extends to continuous-time systems without re-derivation, since the kernel identities are stated for both discrete and continuous time.","The approach can enrich the stable spline kernel with a low-order parametric component, so it can represent ARMAX systems with a small parametric part while retaining the closed-form imputer.","In the experiments, imputation plus stable spline estimation matches or beats the full-data stable spline estimator and full-data PEM, suggesting that the imputed values do not degrade subsequent model quality.","Because the estimate does not depend on the scale factor $\\lambda$, imputation only requires relative hyperparameters, reducing the tuning burden."],"supporting_citations":[{"why":"Defines the stable spline kernel $K(s,t)$ of Eq. (11) that serves as the prior covariance for impulse responses.","marker":"[3]"},{"why":"Supplies the nonparametric Gaussian regression predictor and the marginal-likelihood machinery that Algorithm 4 inherits after imputation.","marker":"[4]"},{"why":"Provides the Gaussian process regression framework and kernel covariance language used throughout the paper.","marker":"[5]"},{"why":"Gives the parametric maximum-likelihood ARMA estimation under missing observations that the PEM+Oracle (missing) baseline uses.","marker":"[9]"},{"why":"Introduces the EM algorithm, the standard imputation baseline for missing-data problems that the paper compares against conceptually.","marker":"[16]"},{"why":"Supplies the standard formula for the best linear mean squared estimator used in Proposition 3 and the appendix proof.","marker":"[20]"}],"fun_headline_variants":["Stable spline kernels fill missing outputs, beat full-data PEM","Closed-form imputation via stable splines tops oracle PEM","New Bayesian kernel solves missing data in linear systems","Imputation from stable splines outperforms PEM on gaps","Kernel identity yields imputation that beats parametric rivals"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The estimator is the minimum variance linear one only when the disturbance process $\\xi(t)$ is zero-mean stationary with covariance exactly $\\lambda R$, which holds if the noise-path impulse response $f_p$ is drawn from a Gaussian process with stable spline covariance using the same stability decay $\\beta$ as the system's impulse responses.","fun_headline_variants_meta":{"raw":{"variants":["Stable spline kernels fill missing outputs, beat full-data PEM","Closed-form imputation via stable splines tops oracle PEM","New Bayesian kernel solves missing data in linear systems","Imputation from stable splines outperforms PEM on gaps","Kernel identity yields imputation that beats parametric rivals"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000118,"raw_usage":{"total_tokens":1068,"prompt_tokens":918,"completion_tokens":150,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":534,"completion_tokens_details":{"reasoning_tokens":70}},"tokens_in":534,"tokens_out":150,"duration_ms":3060,"temperature":1.0,"reasoning_tokens":70,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-14T13:58:55.647829+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run a Monte Carlo with an ARMAX system whose noise-path pole is much faster than the system's least-damped pole, so the shared-$\\beta$ RBF kernel is misspecified; if the closed-form imputer's mean-squared error on missing outputs is not close to that of the best linear estimator computed with the true noise covariance, the shared-$\\beta$ assumption is load-bearing. A simpler check is to estimate the residual autocovariance from a long simulation and compare its shape with Eq. (15).","supporting_citations":[{"cited_title":"A new kernel-based approach for linear system identification,","cited_arxiv_id":null,"evidence_quote":"Defines the stable spline kernel $K(s,t)$ of Eq. (11) that serves as the prior covariance for impulse responses."},{"cited_title":"Prediction error identification of linear systems: A nonparametric Gaussian regression approach,","cited_arxiv_id":null,"evidence_quote":"Supplies the nonparametric Gaussian regression predictor and the marginal-likelihood machinery that Algorithm 4 inherits after imputation."},{"cited_title":"Rasmussen and C","cited_arxiv_id":null,"evidence_quote":"Provides the Gaussian process regression framework and kernel covariance language used throughout the paper."},{"cited_title":"Maximum likelihood fitting of arma models to time series with observations,","cited_arxiv_id":null,"evidence_quote":"Gives the parametric maximum-likelihood ARMA estimation under missing observations that the PEM+Oracle (missing) baseline uses."},{"cited_title":"Maximum likelihood from incomplete data via the em algorithm,","cited_arxiv_id":null,"evidence_quote":"Introduces the EM algorithm, the standard imputation baseline for missing-data problems that the paper compares against conceptually."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the standard formula for the best linear mean squared estimator used in Proposition 3 and the appendix proof."}],"review_version":1}