Poly-Dual Polynomial Regression

• The paper presents a hybrid piecewise regression approach that combines GPU-accelerated pilot density estimation with dual polynomial fitting to accurately capture asymmetric, heavy-tailed unimodal distributions.
• It demonstrates superior computational efficiency and fitting accuracy compared to conventional KDE methods, significantly reducing memory usage and speeding up inference.
• The method offers practical advantages for large-scale datasets in clinical, financial, and reliability applications by adapting to skewed data with minimal tuning.

Poly-Dual Polynomial Regression (Poly-DPR) is a piecewise regression framework developed for accurate, efficient estimation of asymmetric, heavy-tailed, unimodal probability densities. Leveraging a hybrid procedure, Poly-DPR first constructs a smooth pilot density estimate using GPU-accelerated kernel or histogram-based surrogates, then fits two separate polynomials (of a chosen order) to the empirical density function on either side of the mode. This architecture enables low-cost evaluation and strong data-adaptive flexibility, especially for real-world distributions that exhibit skewness or exponential decay and are not well-modeled by classical parametric or generic nonparametric density estimators (Sarkar et al., 3 Dec 2025).

1. Motivation and Problem Setting

Standard parametric families, including the Normal and Pearson Type I distributions, impose restrictive symmetry and tail constraints, resulting in poor performance on skewed or heavy-tailed data. Conversely, generic nonparametric methods such as histograms or kernel density estimation (KDE) lack structure yet suffer bin-edge artifacts or become computationally intractable due to $O(N \cdot M)$ scaling in both time and space (where $N$ is the sample size and $M$ is the evaluation grid). Poly-DPR addresses this gap by utilizing the strengths of piecewise modeling—flexibility to accommodate empirical irregularities—while enforcing computational efficiency and numerical tractability, particularly for very large datasets and real-time inference scenarios (Sarkar et al., 3 Dec 2025).

2. Mathematical Formulation

Given one-dimensional observations $x_1, \ldots, x_N$ from an unknown unimodal density $f(x)$, Poly-DPR executes in two stages. First, a smooth, empirical pilot estimate $\hat f(x)$ is constructed using either:

• Tensor-based KDE (tKDE):

$$\hat f(x_k) = \frac{1}{N h} \sum_{i=1}^N K\left( \frac{x_k - x_i}{h} \right),$$

where $h$ is set by Scott’s rule $h = \sigma N^{-1/5}$.

• Tensor-based Histogram Density Estimation (tHDE), smoothed by a Gaussian filter.

The mode $m$ is determined as $N$0. The domain is split at $N$1, and two independent polynomials of order $N$2 are fitted on the left ($N$3) and right ($N$4) segments, in shifted coordinates $N$5:

• For $N$6:

$N$7

• For $N$8:

$N$9

Continuity is automatically enforced at the mode ($M$0) because both polynomials agree at $M$1 (i.e., $M$2). No derivative matching is imposed, ensuring numerical stability and simplicity. The polynomial coefficients $M$3 are obtained by minimizing mean squared error between the polynomial fit and the pilot estimate over all grid points:

$M$4

Normalization is enforced such that the fitted PDF integrates to unity over the grid:

$M$5

Values are floored at $M$6 and clipped to not exceed $M$7 to ensure numerical stability (Sarkar et al., 3 Dec 2025).

3. Implementation Aspects and Computational Complexity

tKDE and tHDE are implemented using TensorFlow 2.x, exploiting GPU batch operations and @tf.function(reduce_retracing=True). Least-squares polynomial fitting is performed using tf.linalg.lstsq (Cholesky-based solver) on independent Vandermonde matrices for each segment:

$M$8

The core workflow consists of generating the pilot estimate, mode finding, data partitioning, separate polynomial least-squares fits, clipping and normalization.

Theoretical computational complexity:

Procedure	Training Time	Training Memory	Evaluation
tKDE + Poly-DPR	$M$9	$x_1, \ldots, x_N$0	$x_1, \ldots, x_N$1
tHDE + Poly-DPR	$x_1, \ldots, x_N$2	$x_1, \ldots, x_N$3	$x_1, \ldots, x_N$4
SciPy KDE	$x_1, \ldots, x_N$5	$x_1, \ldots, x_N$6	$x_1, \ldots, x_N$7

Poly-DPR yields major efficiency advantages both in memory and evaluation latency, supporting deployment on datasets for which full KDEs are infeasible (Sarkar et al., 3 Dec 2025).

4. Model Order Selection and Hyperparameters

Polynomial order $x_1, \ldots, x_N$8 sets the balance between fit complexity and stability. Empirical benchmarking across six synthetic unimodal distributions (including symmetric/skewed Gaussians, AMW-I, AMW-II, Asymmetric Laplace) indicates:

• Orders $x_1, \ldots, x_N$9 underfit asymmetric shapes.
• Orders $f(x)$0 often trigger Cholesky instability due to ill-conditioning in Vandermonde matrices.
• $f(x)$1 achieves optimal balance: low mean squared error, Jensen-Shannon divergence (JSD $f(x)$2), and Pearson correlation coefficients near unity ($f(x)$3).

A plausible implication is that, in most practical cases, choosing $f(x)$4 suffices for accuracy and stability; higher orders may require numerically robust solvers (e.g., QR decomposition) (Sarkar et al., 3 Dec 2025).

5. Empirical Evaluation and Benchmarking

Poly-DPR was benchmarked using synthetic and real-world unimodal datasets. Metrics include mean squared error (MSE), Jensen-Shannon divergence (JSD), Pearson $f(x)$5, and area under the curve (AUC). On six synthetic datasets, Poly-DPR ($f(x)$6, tKDE backend) achieved:

• Mean JSD $f(x)$7,
• MSE $f(x)$8,
• $f(x)$9,
• Inference time orders of magnitude faster than SciPy KDE.

On 300,000-patient systolic and diastolic blood pressure data, full-scale SciPy KDE failed due to $\hat f(x)$0 GB memory use, whereas Poly-DPR ($\hat f(x)$1, tKDE) trained in $\hat f(x)$2 seconds (for 1,000 grid points) and inferred in $\hat f(x)$3 seconds for 15,000 points. For systolic data, AUC = 1.00, JSD = 0.005, MSE $\hat f(x)$4, $\hat f(x)$5, and $\hat f(x)$6 faster inference than SciPy KDE. Diastolic data showed similar results. A sliding window analysis confirmed significantly lower median JSD for $\hat f(x)$7 vs.\ $\hat f(x)$8 (Mann–Whitney $\hat f(x)$9, $$\hat f(x_k) = \frac{1}{N h} \sum_{i=1}^N K\left( \frac{x_k - x_i}{h} \right),$$0), with batch inference latency $$\hat f(x_k) = \frac{1}{N h} \sum_{i=1}^N K\left( \frac{x_k - x_i}{h} \right),$$1 ms (Sarkar et al., 3 Dec 2025).

6. Applications, Limitations, and Extensions

Poly-DPR is particularly advantageous for large-scale clinical, financial, or reliability datasets characterized by unimodality, skewness, or heavy tails. Its strengths are efficiency, interpretability, and resilience to deviations from Gaussian assumptions. The main limitation is numerical instability for high polynomial orders ($$\hat f(x_k) = \frac{1}{N h} \sum_{i=1}^N K\left( \frac{x_k - x_i}{h} \right),$$2), beyond which ill-conditioned least-squares systems impede robust fitting. tKDE is preferable if sufficient GPU memory is available; tHDE is the recommended fallback for constrained environments. A plausible implication is that future developments could involve robustification of fitting routines (e.g., QR-based solvers) and auto-adaptive selection of polynomial order to maximize fit quality without sacrificing stability (Sarkar et al., 3 Dec 2025).