Sparse Polynomial Regression Model

• SPRM is a class of regression models that blends polynomial expansions with sparsity constraints to achieve parsimonious and interpretable high-dimensional approximations.
• It leverages techniques such as ℓ1-regularization, hierarchical variable selection, and tensor decompositions to balance model complexity with computational tractability.
• The approach is applied in surrogate modeling, system identification, and genotype–phenotype analysis, illustrating trade-offs between recovery accuracy and sample complexity.

Searching arXiv for the provided SPRM-related papers and closely related work. arxiv_search.query({"search_query":"all:\"Sparse Polynomial Regression\" OR all:\"Sparse Hierarchical Regression with Polynomials\" OR id:(Bertsimas et al., 2017) OR id:(Abolpour et al., 25 Aug 2025) OR id:(Kekatos et al., 2011) OR id:(2002.01290) OR id:(Hibraj et al., 2020) OR id:(Götte et al., 2021) OR id:(Hatstatt et al., 23 Jan 2026) OR id:(Azmi et al., 2020) OR id:(Zhang et al., 2014)","start":0,"max_results":10}) Sparse Polynomial Regression Model (SPRM) denotes a class of regression models in which the predictor is expanded in a polynomial feature dictionary, while sparsity is imposed on the coefficient vector, on the admissible monomials, on the active input variables, or on a structured representation of the polynomial coefficients. Across the literature, this includes direct monomial-basis regression, sparse polynomial chaos expansions (PCE), hierarchical sparse polynomial regression, tensor-kernel $\ell^p$-regularized regression, and structured tensor-network formulations. The common objective is to approximate a target map with a polynomial surrogate that remains parsimonious enough to be estimable in high dimension, interpretable at the level of active terms, and computationally tractable under suitable algorithmic assumptions (Kekatos et al., 2011).

1. Formal model class

A standard multivariate polynomial regression formulation writes the response as

$$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$

where $x^{(k)}\in \mathbb{R}^n$, $y^{(k)}\in \mathbb{R}$, $d$ is the polynomial degree, $$\Gamma_d=\{\alpha\in\mathbb{N}^n : \sum_{i=1}^n \alpha_i \le d\}$$, $m_d = |\Gamma_d|$, and $$(x^{(k)})^\alpha = \prod_{i=1}^n (x_i^{(k)})^{\alpha_i}$$. In this form, SPRM is obtained by constraining or regularizing the coefficient family $\{c_\alpha\}$ so that only a small subset of polynomial terms remains active (Abolpour et al., 25 Aug 2025).

An equivalent linear-in-parameters view is central to the subject. If the polynomial basis functions are assembled into a regression matrix $\Psi\in\mathbb{R}^{n\times P}$, with $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$0, then the coefficient vector $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$1 is estimated from a regression problem in feature space. In full PCE this yields

$$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$2

whereas the sparse version is explicitly written as

$$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$3

This $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$4-regularized formulation is presented as the core sparse polynomial regression objective in sparse PCE-based surrogate modeling (Hatstatt et al., 23 Jan 2026).

A related general formulation in sparse PCE approximates a model output $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$5 by

$$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$6

with $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$7 a finite candidate multi-index set. In this setting, sparse polynomial regression is the problem of estimating a coefficient vector with many zeros from relatively few model evaluations, motivated by the sparsity-of-effects principle and the practical infeasibility of dense ordinary least squares when $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$8 is large relative to the experimental design size (2002.01290).

The same linearized viewpoint also underlies sparse Volterra and polynomial regression models for nonlinear system identification. There the data model is written as

$$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$9

after lifting the nonlinear input-output map into a monomial feature vector. The underlying nonlinearity is then represented by a sparse parameter vector $x^{(k)}\in \mathbb{R}^n$0 in a high-dimensional basis, making compressed-sensing-style estimation applicable (Kekatos et al., 2011).

2. Sparsity structures and polynomial dictionaries

The simplest SPRM uses unstructured coefficient sparsity: a large candidate monomial dictionary is generated, and only a small number of coefficients are allowed to remain nonzero. In sparse PCE, the candidate basis is commonly chosen through total-degree truncation,

$x^{(k)}\in \mathbb{R}^n$1

hyperbolic truncation,

$x^{(k)}\in \mathbb{R}^n$2

or interaction-order truncation,

$x^{(k)}\in \mathbb{R}^n$3

These truncation schemes govern not only expressiveness but also conditioning and the ratio $x^{(k)}\in \mathbb{R}^n$4, both of which materially affect sparse recovery (2002.01290).

A more structured variant is hierarchical sparse polynomial regression. For polynomials $x^{(k)}\in \mathbb{R}^n$5 of total degree at most $x^{(k)}\in \mathbb{R}^n$6, expressed in a monomial basis $x^{(k)}\in \mathbb{R}^n$7 with

$x^{(k)}\in \mathbb{R}^n$8

the class $x^{(k)}\in \mathbb{R}^n$9 restricts the regressor to depend on at most $y^{(k)}\in \mathbb{R}$0 input variables and to contain at most $y^{(k)}\in \mathbb{R}$1 monomial terms. Its defining hierarchy is encoded by

$y^{(k)}\in \mathbb{R}$2

so a monomial may be selected only if every variable appearing in that monomial is itself selected. This yields a disciplined sparse model in which variable selection occurs “above” monomial selection (Bertsimas et al., 2017).

Another structured design arises in high-dimensional polynomial regression via homogeneous polynomial subspaces. With a one-dimensional basis dictionary $y^{(k)}\in \mathbb{R}$3, the full tensor-product space is

$y^{(k)}\in \mathbb{R}$4

whereas the homogeneous subspace of degree $y^{(k)}\in \mathbb{R}$5 is

$y^{(k)}\in \mathbb{R}$6

Here sparsity is induced by degree structure rather than by an arbitrary support pattern: the admissible coefficients occupy only a restricted subspace of the full tensorized polynomial space (Götte et al., 2021).

Tensor-kernel formulations provide a different representation of the same broad idea. For a linear tensor kernel of order $y^{(k)}\in \mathbb{R}$7,

$y^{(k)}\in \mathbb{R}$8

and for the polynomial tensor kernel,

$y^{(k)}\in \mathbb{R}$9

These kernels encode higher-order monomial interactions implicitly. When combined with $d$0-regularization with

$d$1

for even integer $d$2, the penalty lies in $d$3 and can be made arbitrarily close to $d$4, yielding a lasso-like sparse behavior in a kernelized polynomial-interaction model (Hibraj et al., 2020).

3. Estimation, optimization, and computational formulations

A core estimation family is based on $d$5-regularization. In sparse polynomial and Volterra regression, the noisy estimation problem is written as

$d$6

with $d$7 for Lasso and

$d$8

for weighted Lasso. The same work also develops recursive sparse estimation through an exponentially weighted objective,

$d$9

leading to the RL and RWL variants and a CCD-R(W)L algorithm with recursive sufficient statistics (Kekatos et al., 2011).

Sparse PCE adopts a pathwise model-selection workflow centered on LARS. The procedure described for sparse polynomial chaos consists of: constructing a candidate polynomial basis, running LARS to generate a sparse solution path, selecting the best model by leave-one-out cross-validation, and refitting the coefficients by ordinary least squares on the selected active set. In the cited implementation, this is specifically the Hybrid LARS procedure (Hatstatt et al., 23 Jan 2026).

Hierarchical sparse polynomial regression departs from convex relaxation and solves a reduced exact combinatorial problem. It first ranks inputs using ordinary polynomial kernel regression with

$$\Gamma_d=\{\alpha\in\mathbb{N}^n : \sum_{i=1}^n \alpha_i \le d\}$$0

then retains only the top $$\Gamma_d=\{\alpha\in\mathbb{N}^n : \sum_{i=1}^n \alpha_i \le d\}$$1 variables, and finally solves the reduced $$\Gamma_d=\{\alpha\in\mathbb{N}^n : \sum_{i=1}^n \alpha_i \le d\}$$2-sparse problem exactly by a cutting-plane / outer-approximation algorithm implemented with lazy constraints in a branch-and-bound solver. The exact reduced formulation is expressed as a pure binary optimization problem after analytically eliminating the continuous coefficients (Bertsimas et al., 2017).

Robust SPRM under anomalous data is formulated as a mixed-integer model. Binary variables $$\Gamma_d=\{\alpha\in\mathbb{N}^n : \sum_{i=1}^n \alpha_i \le d\}$$3 select monomials and $$\Gamma_d=\{\alpha\in\mathbb{N}^n : \sum_{i=1}^n \alpha_i \le d\}$$4 select non-anomalous samples: $$\Gamma_d=\{\alpha\in\mathbb{N}^n : \sum_{i=1}^n \alpha_i \le d\}$$5 The resulting MILP minimizes the maximum absolute residual $$\Gamma_d=\{\alpha\in\mathbb{N}^n : \sum_{i=1}^n \alpha_i \le d\}$$6, converts to a nonconvex QCQP, and is then mapped to a fractional program (FP). The proposed practical solver is the Two-Step Convex Relaxation and Recovery (TS-CRR) algorithm: first solve the linear-based relaxation of the FP to infer the sparse monomial set and retained samples, then solve a recovery LP for the coefficients on the selected support (Abolpour et al., 25 Aug 2025).

High-dimensional structured SPRM can also be solved in tensor-network form. In the block-sparse Tensor Train (TT) approach, the coefficient tensor is restricted by the homogeneous polynomial condition $$\Gamma_d=\{\alpha\in\mathbb{N}^n : \sum_{i=1}^n \alpha_i \le d\}$$7, which induces a block sparsity pattern in the TT cores. Optimization proceeds by an ALS-type scheme using measurement matrices $$\Gamma_d=\{\alpha\in\mathbb{N}^n : \sum_{i=1}^n \alpha_i \le d\}$$8 and local contractions $$\Gamma_d=\{\alpha\in\mathbb{N}^n : \sum_{i=1}^n \alpha_i \le d\}$$9, so that each core update reduces to a restricted least-squares subproblem on the nonzero block pattern (Götte et al., 2021).

On the implementation side, tensor-kernel sparse regression has been improved through two practical devices: a storage layout that retains only the unique entries of symmetric tensors, and a Nyström-type subsampling strategy that solves the tensor-kernel problem on a subset $m_d = |\Gamma_d|$0 with $m_d = |\Gamma_d|$1. For a fourth-order symmetric tensor, the number of stored entries becomes

$m_d = |\Gamma_d|$2

rather than $m_d = |\Gamma_d|$3, and the code is reported as implemented in C++ (Hibraj et al., 2020).

4. Recoverability, sample complexity, and computational barriers

The statistical motivation for SPRM is that polynomial feature spaces are combinatorially large. In sparse Volterra and polynomial regression, the number of distinct coefficients is

$m_d = |\Gamma_d|$4

which makes dense least squares impractical unless the sample size is correspondingly large. Sparse recovery theory is therefore framed through restricted isometry properties (RIP). For sparse linear-quadratic polynomial regression with independent inputs, one result states that

$m_d = |\Gamma_d|$5

suffices for RIP with high probability, while for second-order non-homogeneous Volterra filtering the sufficient scaling is

$m_d = |\Gamma_d|$6

These results quantify how sparsity can reduce measurement requirements relative to the ambient polynomial dimension (Kekatos et al., 2011).

Sample complexity is also shaped by the choice of polynomial subspace. In the TT-based high-dimensional formulation, the variation constant of the full Legendre product space satisfies

$m_d = |\Gamma_d|$7

whereas for homogeneous degree-$m_d = |\Gamma_d|$8 spaces the bound is much smaller: $m_d = |\Gamma_d|$9 With optimally weighted sampling, the improved scaling

$$(x^{(k)})^\alpha = \prod_{i=1}^n (x_i^{(k)})^{\alpha_i}$$0

is reported. This indicates that homogeneous polynomial ansatz spaces are statistically easier to learn than the full tensor-product polynomial space (Götte et al., 2021).

Sparse PCE benchmarking reaches a closely related conclusion from the empirical side. Over 11 models, 30–50 replications per setting, multiple experimental design sizes, and several solver/sampling combinations, the choice of sparse regression solver and sampling scheme can change the resulting relative mean-squared error by orders of magnitude. The benchmark reports regime dependence rather than a universal winner: BCS is strongest in high-dimensional or low-information settings, SP often becomes best in low-dimensional problems with larger experimental designs, and advanced sampling schemes help mainly in low-dimensional basis-rich regimes (2002.01290).

A separate line of theory identifies a computational-statistical gap. For sparse linear regression,

$$(x^{(k)})^\alpha = \prod_{i=1}^n (x_i^{(k)})^{\alpha_i}$$1

the $$(x^{(k)})^\alpha = \prod_{i=1}^n (x_i^{(k)})^{\alpha_i}$$2-constrained estimator achieves

$$(x^{(k)})^\alpha = \prod_{i=1}^n (x_i^{(k)})^{\alpha_i}$$3

whereas thresholded Lasso under the restricted eigenvalue condition incurs

$$(x^{(k)})^\alpha = \prod_{i=1}^n (x_i^{(k)})^{\alpha_i}$$4

with high probability. Under the assumption

$$(x^{(k)})^\alpha = \prod_{i=1}^n (x_i^{(k)})^{\alpha_i}$$5

the lower bound shows that polynomial-time-efficient estimators can be forced to suffer a prediction risk of order

$$(x^{(k)})^\alpha = \prod_{i=1}^n (x_i^{(k)})^{\alpha_i}$$6

on hard ill-conditioned designs, while the optimal combinatorial estimator avoids the $$(x^{(k)})^\alpha = \prod_{i=1}^n (x_i^{(k)})^{\alpha_i}$$7 penalty (Zhang et al., 2014). The cited work is stated to be directly relevant to sparse polynomial regression more broadly; this suggests that SPRM-like methods based on polynomial-time optimization may inherit a fundamental limitation on ill-conditioned polynomial design matrices.

5. Empirical regimes and application domains

SPRM has been used in nonlinear system identification, genotype–phenotype analysis, surrogate modeling, optimal control, and data-driven forecasting. In sparse Volterra and polynomial regression, reported applications include neuroscience spike-train connectivity models, genotype–phenotype or epistasis models, and nonlinear system-identification settings such as LNL/Wiener/Hammerstein-type structures. In a synthetic genotype–phenotype experiment with $$(x^{(k)})^\alpha = \prod_{i=1}^n (x_i^{(k)})^{\alpha_i}$$8, $$(x^{(k)})^\alpha = \prod_{i=1}^n (x_i^{(k)})^{\alpha_i}$$9, and $\{c_\alpha\}$0, the reported predictive error, MSE, and number of nonzeros were $\{c_\alpha\}$1, $\{c_\alpha\}$2, and $\{c_\alpha\}$3 for ridge; $\{c_\alpha\}$4, $\{c_\alpha\}$5, and $\{c_\alpha\}$6 for Lasso; and $\{c_\alpha\}$7, $\{c_\alpha\}$8, and $\{c_\alpha\}$9 for weighted Lasso. On a real barley height dataset, the reported prediction errors were $\Psi\in\mathbb{R}^{n\times P}$0 for ridge, $\Psi\in\mathbb{R}^{n\times P}$1 for Lasso, and $\Psi\in\mathbb{R}^{n\times P}$2 for weighted Lasso (Kekatos et al., 2011).

Exact hierarchical sparse polynomial regression reports a phase transition in recoverability. In one synthetic cubic example with $\Psi\in\mathbb{R}^{n\times P}$3, $\Psi\in\mathbb{R}^{n\times P}$4, and $\Psi\in\mathbb{R}^{n\times P}$5, recovery becomes essentially perfect around $\Psi\in\mathbb{R}^{n\times P}$6, and the same transition is reflected in accuracy $\Psi\in\mathbb{R}^{n\times P}$7, false alarms $\Psi\in\mathbb{R}^{n\times P}$8, runtime, and the number of cutting planes needed. With the ranking heuristic, the method is reported to scale to problems with $\Psi\in\mathbb{R}^{n\times P}$9 observations and $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$00 inputs (Bertsimas et al., 2017).

In surrogate modeling, sparse PCE is explicitly treated as a sparse polynomial regression model. The benchmark literature emphasizes that different methods excel in different regimes: BCS performs best most often on small experimental designs in low dimension and is clearly strongest in high-dimensional models, while SP often becomes the best overall solver in large-design low-dimensional settings. LHS is usually slightly better than MC, and near-optimal sampling can be best when feasible (2002.01290).

Gradient-augmented sparse polynomial regression has been used for optimal feedback law recovery. There, a polynomial surrogate $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$01 for the value function is fitted from offline state-value and state-gradient data, and the feedback is recovered via

$$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$02

The reported numerical tests include a Van der Pol oscillator, a controlled Allen–Cahn equation, and a Cucker–Smale consensus problem in state dimension $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$03. In the Van der Pol example, adding the $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$04 penalty reduced the number of active coefficients from $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$05 to $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$06, with a reported online floating-point-operation reduction of about $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$07. In the Cucker–Smale test, gradient-augmented sparse regression required roughly two orders of magnitude fewer samples than gradient-free sparse regression for comparable error levels (Azmi et al., 2020).

Robust SPRM with anomalous-data filtering has been evaluated on an electricity price dataset and a temperature forecasting dataset. For the electricity data, the setup uses polynomial degree $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$08, $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$09 training points, $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$10 candidate monomials, $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$11 anomalous points removed, and $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$12 selected monomials. For the temperature dataset, the setup uses $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$13, $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$14, $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$15, $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$16, and $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$17. In both cases TS-CRR is reported as giving stable interpolation and extrapolation performance, and the paper emphasizes its balanced generalization relative to linear regression, polynomial regression, and decision-tree-type baselines (Abolpour et al., 25 Aug 2025).

The following summary organizes the main empirical problem settings.

Domain	SPRM formulation	Reported empirical theme
System identification and genotype–phenotype analysis	Sparse polynomial / Volterra regression with Lasso or weighted Lasso	Weighted Lasso yields sparser models and improved estimation accuracy (Kekatos et al., 2011)
High-dimensional nonlinear regression	Exact hierarchical $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$18-sparse polynomial regression	Phase transition in support recovery; scalability to $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$19, $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$20 (Bertsimas et al., 2017)
Surrogate modeling and control	Sparse PCE, gradient-augmented sparse polynomial regression	Strong regime dependence; gradient data reduce sample complexity (2002.01290)

6. Robustness, uncertainty quantification, and methodological cautions

A recurrent robustness issue is that sparse polynomial structure is data-dependent. In anomalous-data SPRM, this is handled explicitly by binary sample-selection variables $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$21, so anomaly filtering is not a preprocessing step but part of the optimization model itself. The retained sample budget $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$22 and monomial budget $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$23 jointly determine the final robust sparse polynomial surrogate (Abolpour et al., 25 Aug 2025).

Uncertainty quantification for sparse polynomial surrogates has been developed through conformal prediction. For PCE-based surrogate models, two conformal approaches are integrated into both full and sparse PCEs: full conformal and Jackknife+. For full PCE, the fixed-basis OLS structure permits computational shortcuts through Sherman–Morrison updates and hat-matrix formulas, with residuals affine in the trial value. For sparse PCE, however, a naive approach that selects a sparse basis once and then freezes it is not symmetric with respect to the data and therefore breaks the conformal validity assumptions. The cited work states that this naive shortcut leads to poor empirical coverage (Hatstatt et al., 23 Jan 2026).

The same study distinguishes the valid and approximate sparse conformal strategies. Jackknife+ for sparse PCE remains valid when Hybrid LARS is rerun separately in each leave-one-out fit. Full conformal for sparse PCE is computationally harder, and the reported tractable implementation fixes the number of active regressors, maps the LARS solution to a pseudo-regularization parameter $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$24, uses the LASSO homotopy path as a surrogate, and then applies Brent’s method in a one-dimensional optimization. The paper is explicit that fixing $$y^{(k)}=\sum_{\alpha \in \Gamma_d} c_\alpha (x^{(k)})^\alpha + \epsilon^{(k)}, \qquad k=1,\dots,N,$$25 sacrifices strict finite-sample guarantees, although empirical coverage remains close to nominal, and both conformal approaches are reported as achieving better-calibrated prediction intervals and superior coverage relative to bootstrap while maintaining moderate computational cost (Hatstatt et al., 23 Jan 2026).

Several methodological cautions follow from the combined literature. First, sparsity does not by itself remove computational hardness: exact support selection can require mixed-integer optimization, cutting planes, or combinatorial search, even when the final model is small (Bertsimas et al., 2017). Second, polynomial-time sparse estimators can be fundamentally suboptimal on ill-conditioned designs, so poor performance is not always attributable to a weak implementation; this suggests a genuine computational-statistical barrier rather than merely an algorithm-engineering deficit (Zhang et al., 2014). Third, solver quality depends strongly on basis truncation, sampling scheme, and the structural assumptions encoded in the polynomial space, rather than on sparsity regularization alone (2002.01290).

Taken together, the literature portrays SPRM not as a single algorithm but as a modeling paradigm: polynomial approximation constrained by sparsity, hierarchy, robustness, or low-rank structure, with performance governed jointly by basis design, optimization strategy, sample regime, and the conditioning of the induced design matrix.