respolygp 0.1.2

ResPolyGP: Orthogonal Legendre Trajectory Gaussian Process Prior Model

ResPolyGP

ResPolyGP is a Python library for Gaussian Process (GP) trajectory modeling that resolves standard exact GP computational bottlenecks and non-convex initialization challenges on non-stationary, cyclical sequence extrapolation.

By formulating secular polynomial trends as a structured finite-rank prior over Legendre orthogonal basis functions ($K = \Phi A \Phi^\top + \sigma_n^2 I$), ResPolyGP utilizes the Woodbury Matrix Identity and Matrix Determinant Lemma to project the standard cubic exact GP complexity onto a tiny basis space.

This accelerates exact GP polynomial trend inference from $O(N^3)$ time and $O(N^2)$ memory to $O(ND^2 + D^3)$ time and $O(ND + D^2)$ memory (where $D$ is the active basis degree, typically $D \le 5 \ll N$), delivering up to a 1000x computational speedup without any approximation.

---

Key Features

1. Orthogonal Bounded Legendre Kernels: Replaces standard unstable monomial trend kernels with recursively computed Legendre orthogonal polynomials mapped onto $x \in [-1, 1]$ to prevent floating-point underflow/overflow and preserve conditioning.
2. True $O(ND^2 + D^3)$ Woodbury Solver: Employs Woodbury matrix inversion and matrix determinant lemma in low-rank space, bypassing the Cholesky bottleneck and completely preventing NotPSDError.
3. Sequential Significance Early Stopping: Filters polynomial degrees dynamically via a variance-explained significance threshold ($q_d < 0.01$) to prevent overfitting to local noise crests.
4. Simplex-Constrained Softmax Priors: Normalizes and refines components on the simplex ($\sum \alpha_k = 1.0$), ensuring structural regularization.
5. Shape Parameter Freezing / Weights-Only Tuning: Locks shape parameters (e.g. Lomb-Scargle periods and ACF stochastic decay lengths) at their analytical time-domain values, optimizing only mixture weights to eliminate optimization drift.
6. Spectral Healing: Detrends secular trends analytically to restore stationarity before Lomb-Scargle periodic detection, reducing periodic extrapolation error by over 80%.

---

Installation

You can install ResPolyGP directly from PyPI:

pip install respolygp

Dependencies

• Python >= 3.8
• numpy
• scipy
• torch
• gpytorch

---

Quickstart

Below is a simple example showing how to initialize, fit, and predict with ResPolyGPModel using the fast Woodbury Representation (default), the standard GPyTorch backend, or the Traditional GP baseline.

1. Legendre Trend Model (Woodbury Representation - Fast $O(N)$)

The Woodbury solver is the default backend. It reduces GP computation to linear complexity $O(ND^2 + D^3)$.

import numpy as np
from respolygp import ResPolyGPModel

# Generate a noisy cubic trajectory
np.random.seed(42)
t_raw = np.linspace(10, 50, 100)
t_norm = (t_raw - 10) / 40.0 * 2.0 - 1.0
y_raw = 5.0 + 3.0 * t_norm - 2.0 * t_norm**2 + 4.0 * t_norm**3 + np.random.normal(0, 0.1, 100)

# Instantiate Legendre model with Woodbury backend (use_woodbury=True is default)
model_woodbury = ResPolyGPModel(max_degree=5, sig_threshold=0.01, use_woodbury=True)

# Run Legendre decomposition and initialize hyperparameters
model_woodbury.initialize_hyperparameters(t_raw, y_raw)

# Refine simplex weights via Adam optimizer (only 15-20 epochs needed)
history = model_woodbury.fit_gp(t_raw, y_raw, lr=0.05, epochs=20)

# Predict out-of-sample in raw physical scale units
t_test = np.array([55.0, 60.0])
mean, var, lower, upper = model_woodbury.predict_raw(t_test)
print(f"Woodbury predictions at {t_test}: Mean={mean}, Var={var}")

2. Legendre Trend Model (Standard GPyTorch Backend)

If you prefer to use GPyTorch's exact Cholesky solver ($O(N^3)$), configure use_woodbury=False. Both backends are mathematically equivalent but scale differently:

# Instantiate Legendre model with standard GPyTorch backend
model_gp = ResPolyGPModel(max_degree=5, sig_threshold=0.01, use_woodbury=False)
model_gp.initialize_hyperparameters(t_raw, y_raw)
model_gp.fit_gp(t_raw, y_raw, lr=0.05, epochs=20)

mean_gp, var_gp, _, _ = model_gp.predict_raw(t_test)
print(f"GPyTorch predictions at {t_test}: Mean={mean_gp}, Var={var_gp}")

3. Additive Trend + Periodic System

For non-stationary, cyclical sequence extrapolation, use the MostlyFinalAdditiveGPModel whichDetrends secular trends analytically to restore stationarity before Lomb-Scargle periodic detection:

from respolygp import MostlyFinalAdditiveGPModel

model_additive = MostlyFinalAdditiveGPModel(max_degree=5, sig_threshold=0.01)
model_additive.initialize_hyperparameters(t_raw, y_raw)
model_additive.fit_gp(t_raw, y_raw, lr=0.05, epochs=20)

mean_add, var_add, _, _ = model_additive.predict_raw(t_test)

4. Traditional GP Baseline Model

To compare Legendre polynomials against standard RBF/Periodic covariance models, use TraditionalGPModel:

from respolygp import TraditionalGPModel

# Traditional GP baseline (Constant + Linear + Periodic + RBF)
model_trad = TraditionalGPModel()
model_trad.initialize_hyperparameters(t_raw, y_raw)
model_trad.fit_gp(t_raw, y_raw, lr=0.05, epochs=50)

mean_trad, var_trad, _, _ = model_trad.predict_raw(t_test)

---

License

This project is licensed under the GNU General Public License v3 (GPLv3) - see the LICENSE file for details.