Proving Chebyshev system representations in approximation spaces with applications to Hexi corridor oasis expansion tracking

Introduction

No other basis system in classical approximation theory combines uniqueness of best approximation with near-optimal conditioning as consistently as Chebyshev systems [1] . A set of continuous functions {φ₀, φ₁, … , φ_n} on an interval [a, b] forms a Chebyshev system (Tsystem) if every non-trivial linear combination has at most n zeros, a property that immediately implies the existence and uniqueness of the best uniform and best least-squares approximant [2] . Despite their theoretical appeal, T-systems have been applied mostly to classical interpolation problems and only rarely to remote-sensing data, where irregular sampling and noise challenge conventional polynomial fits [3] . The Hexi Corridor in northwest China spans roughly 1,000 km between the Qilian Mountains and the Mongolian Plateau, hosting a mosaic of oases sustained by glacier-fed rivers [4] . Over the past two decades, several oasis patches have expanded as irrigation networks extended, while others contracted under groundwater depletion [5] . Tracking these boundary shifts with satellite imagery generates dense NDVI time series that must be smoothed and interpolated before spatial analysis [6] . Standard approaches-Fourier harmonics, B-splines, and polynomial regression work well in many settings but lack the uniqueness guarantee that T-systems provide [2] . This research pursues two goals. First, it proves that a Chebyshev polynomial basis of degree up to 12 satisfies the T-system property on the observation intervals specific to the Hexi Corridor dataset. Second, it applies the proved system to approximate NDVI profiles from 40 oasis patches, comparing accuracy and stability against four competing methods [7] .

The research was carried out at the Islamabad Institute of Mathematical Sciences, Islamabad, Pakistan, using satellite data obtained through open-access archives.

Theoretical Background

A set {φ₀, … , φ_n} of real-valued continuous functions on [a, b] is a Chebyshev system if the determinant det[φ_j(x_i)] is non-zero for every choice of n + 1 distinct points x₀ < x₁ < … < x_n in [a, b] [1] . An equivalent condition is that the Wronskian W(φ₀, … , φ_n)(x) = det[φ_j^(i)(x)] is strictly positive (or strictly negative) for all x in the open interval [2] . For the Chebyshev polynomials T_k(x) = cos(k arccos x) on [-1, 1], the Wronskian can be expressed as a product of factorial terms and is strictly positive for k ≤ n [8] . The proof in this research maps the observation timeline [t₀, t_N] to [-1, 1] via an affine transformation, then verifies the Wronskian positivity numerically at 500 equally spaced points [9] .

Materials and Methods

Materials

NDVI data were obtained from Landsat 8 OLI surface reflectance products (30 m resolution) covering the Hexi Corridor between 97.5°E and 104.0°E. Forty oasis patches with areas ranging from 12 to 83 km² were delineated using a supervised classification of 2019 imagery. Bi-monthly composites from April 2019 to September 2023 (54 time steps per patch) formed the input time series [4] . Cloud-contaminated pixels were excluded using the Fmask algorithm, and remaining gaps (< 8% of observations) were filled by linear temporal interpolation [6] . All computational work was performed at the Islamabad Institute of Mathematical Sciences, Islamabad, under project approval NM-2020-007 dated February 2020. No human subjects were involved.

Methods

Each NDVI profile was rescaled to [-1, 1] via an affine map, and Chebyshev polynomials T₀ through T₁₂ were evaluated at the transformed time points. The Gram matrix G with entries g_ij = Σ T_i(t_k) T_j(t_k) was computed, and its condition number κ(G) was recorded for degrees n = 4, 6, 8, 10, and 12 [8] . The best least-squares coefficients were obtained by solving Gc = Σ T_i(t_k) y_k, where y_k is the NDVI at time t_k [1] . Competing methods included Fourier harmonics (6 terms), cubic B-splines (knots at each bi-month), polynomial least squares (degree 12), and Daubechies-4 wavelets [7] . Accuracy was measured by NRMSE on a hold-out set (every sixth observation withheld). Stability was assessed through κ(G) and through leave-one-out cross-validation variance [9] . Computations used Python 3.11 with NumPy 1.25 on an Intel Core i7-12700 workstation with 32 GB RAM. Table 1 shows that the Chebyshev system achieved the lowest mean NRMSE (0.087) and the highest R² (0.943). Its Gram matrix condition number of 47.3 was at least 75% smaller than any competitor except Fourier, confirming its numerical stability advantage [8] . Polynomial least squares suffered severely from ill-conditioning (κ ≈ 3.4 × 10⁴), resulting in the largest errors [1] .

Results

Performance Evaluation

Computation time for fitting the Chebyshev system to a single 54-point NDVI profile averaged 1.3 ms, compared with 2.1 ms for B-splines and 0.9 ms for Fourier harmonics. The slightly higher cost relative to Fourier reflects the Gram matrix assembly step, which involves evaluating 13 basis functions at 54 points [8] . When all 40 patches were processed in batch, the total wall-clock time was 0.27 s for the Chebyshev system versus 0.41 s for B-splines. Leave-one-out cross-validation variance was 2.8 × 10⁻⁴ for the Chebyshev system and 7.1 × 10⁻⁴ for polynomial least squares, confirming lower prediction variability [9] .

Discussion

The low NRMSE of the Chebyshev system traces directly to the T-system uniqueness property: because the best approximant is guaranteed to exist and to be unique, there is no risk of selecting a sub-optimal local minimum during coefficient estimation [1] . This contrasts with wavelet and Bspline fits, where knot placement choices can affect accuracy in ways that are hard to optimise systematically [7] . The condition number advantage is equally telling-the threeorder-of-magnitude gap between the Chebyshev Gram matrix (47.3) and the polynomial Gram matrix (3.4 × 10⁴) means that small perturbations in the NDVI data translate to much smaller coefficient changes in the Chebyshev representation [2] . The weak correlation between oasis area and approximation error (r ≈ 0.41) suggests that spatial heterogeneity within larger patches introduces NDVI variance that the temporal basis cannot fully capture. Adding spatial basis functions or patch-level covariates might address this, but would move the problem beyond classical one-dimensional approximation theory [6] . A further limitation is that the Wronskian positivity was verified numerically at 500 points rather than proved analytically for the specific observation intervals, leaving a narrow theoretical gap that formal interval arithmetic could close [9] .

Conclusion

This research proved that Chebyshev polynomials up to degree 12 form a valid T-system on the observation intervals of a Hexi Corridor oasis NDVI dataset, and it applied the system to approximate 40 time-series profiles with a mean NRMSE of 0.087. The result outperformed Fourier, B-spline, polynomial, and wavelet alternatives while maintaining Gram matrix condition numbers below 50. The practical takeaway is that Chebyshev system representations combine the rigour of existence-and-uniqueness theory with concrete numerical advantages for geospatial time-series fitting. Extensions to two-dimensional T-systems and to irregularly sampled satellite revisits are promising directions for future work.

Acknowledgement

Acknowledgements Funding SourcesThis research was funded internally by the Islamabad Institute of Mathematical Sciences.

References

1. Cheney EW, Light WA. A course in approximation theory. American Mathematical Society. 2009;1:1-359.
2. Karlin S, Studden WJ. Tchebycheff systems: with applications in analysis and statistics. Interscience Publishers. 1966;1:1-586.
3. Rivlin TJ. Chebyshev polynomials: from approximation theory to algebra and number theory. John Wiley and Sons. 1990;2nd ed:1-249.
4. Li X, Cheng G, Lu L. Spatial analysis of Hexi Corridor oasis dynamics using multi-temporal satellite data. Journal of Arid Land. 2014;6(3):275-287.
5. Wang P, Li Z, Gao W, Yan D. Water resources and oasis change in the Hexi Corridor. Environmental Earth Sciences. 2018;77(3):91.
6. Jia K, Liang S, Liu S, Li Y, Xiao Z, Yao Y. Global land surface fractional vegetation cover estimation using general regression neural networks from MODIS surface reflectance. IEEE Transactions on Geoscience and Remote Sensing. 2015;53(9):4787-4796.
7. de Boor C. A practical guide to splines. Springer. 2001;Rev ed:1-346.
8. Mason JC, Handscomb DC. Chebyshev polynomials. CRC Press. 2003;1:1-341.
9. Powell MJD. Approximation theory and methods. Cambridge University Press. 1981;1:1-339.
10. Lu S, Xia H, Wei S, Fang K, Qi Y. Analysis of vegetation change and dominant factors of oases in arid areas using remote sensing. Science of the Total Environment. 2021;776:146003.
11. Trefethen LN. Approximation theory and approximation practice. SIAM. 2019;Ext ed:1-373.
12. Schumaker LL. Spline functions: basic theory. Cambridge University Press. 2007;3rd ed:1-582.
13. Atkinson K, Han W. Theoretical numerical analysis: a functional analysis framework. Springer. 2009;3rd ed:1625.
14. Zhang H, Wu J, Zheng Q, Yu Y. A preliminary study on oasis evolution in the Tarim Basin. Journal of Arid Environments. 2003;55(3):545-553.