Approximation

Overview

The sangi approximation module collects algorithms that approximate data or functions with polynomials and rational functions. It splits into the following five categories by purpose.

• Regression — fit a line $y = ax + b$ to observed points $(x_i, y_i)$ in the least-squares sense (the linearRegression family)
• Polynomial fitting — fit a degree-$d$ polynomial to observed points in the least-squares sense (polynomialFit)
• Chebyshev approximation — expand a known function $f$ on an interval $[a, b]$ as a Chebyshev polynomial series for a near-minimax approximation (ChebyshevApprox)
• Padé approximation — build a rational function $p(x)/q(x)$ from a Taylor coefficient sequence, improving accuracy over a wider range than a plain truncated Taylor series (the padeApprox family)
• Continued fractions and rational interpolation — the qd algorithm that converts a Taylor series into a continued fraction (taylorToSFraction / taylorToJFraction), and Thiele-type rational interpolation at distinct sample points (rationalInterpolation)

All templates live in namespace sangi, and the element type T must satisfy the concept concepts::OrderedField (e.g. double or the arbitrary-precision Float).

Build

#include <math/approx/approximation.hpp>

Header-only. No library linkage required. Internally it depends on the linear-algebra module (algorithms::solve / algorithms::expm) to solve linear systems, but those are likewise self-contained in the bundled headers.

Linear Regression

Fits a line $y = ax + b$ to observed points in the least-squares sense. The result is returned as a LinearRegressionResult<T> holding the slope, intercept, and coefficient of determination $R^2$.

LinearRegressionResult<T>

template<concepts::OrderedField T>
struct LinearRegressionResult {
    T slope;       // slope a
    T intercept;   // intercept b   (y = a*x + b)
    T r_squared;   // coefficient of determination R²
};

Function reference

linearRegression

// (1) Both x and y supplied
template<concepts::OrderedField T>
LinearRegressionResult<T> linearRegression(std::span<const T> x, std::span<const T> y);
template<concepts::OrderedField T>
LinearRegressionResult<T> linearRegression(const std::vector<T>& x, const std::vector<T>& y);

// (2) Only y supplied (x = 0, 1, ..., n-1 generated automatically)
template<concepts::OrderedField T>
LinearRegressionResult<T> linearRegression(std::span<const T> y);
template<concepts::OrderedField T>
LinearRegressionResult<T> linearRegression(const std::vector<T>& y);

Behaviour: solves for $a, b$ in $y = a x + b$ by standard least squares. After subtracting the means $\bar x, \bar y$, it computes $S_{xx} = \sum (x_i - \bar x)^2$ and $S_{xy} = \sum (x_i - \bar x)(y_i - \bar y)$, then sets $a = S_{xy}/S_{xx}$ and $b = \bar y - a\,\bar x$. The $y$-only overload treats the independent variable as the index $i = 0, 1, \dots, n-1$, useful for estimating the slope of equally spaced data.

Parameters:

Note: throws std::invalid_argument if the sizes of $x$ and $y$ differ, if fewer than 2 points are given, or if all $x$ are identical ($S_{xx} = 0$).

// Example: fit the perfect line y = 2x + 1
std::vector<double> x = {0, 1, 2, 3, 4};
std::vector<double> y = {1, 3, 5, 7, 9};
auto r = sangi::linearRegression(x, y);
std::cout << r.slope << ", " << r.intercept << ", " << r.r_squared;
// Result: 2, 1, 1   (slope=2, intercept=1, R²=1)

weightedLinearRegression

template<concepts::OrderedField T>
LinearRegressionResult<T> weightedLinearRegression(
    std::span<const T> x, std::span<const T> y, std::span<const T> weights);
template<concepts::OrderedField T>
LinearRegressionResult<T> weightedLinearRegression(
    const std::vector<T>& x, const std::vector<T>& y, const std::vector<T>& weights);

Behaviour: least squares with a weight $w_i$ on each point, $\min \sum w_i (y_i - a x_i - b)^2$. The normal equations $\begin{bmatrix} \sum w x^2 & \sum w x \\ \sum w x & \sum w \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} \sum w x y \\ \sum w y \end{bmatrix}$ are solved stably by a direct $2 \times 2$ method (LU). $R^2$ is computed from the weighted residual sum of squares, $1 - \mathrm{SS}_{\mathrm{res}} / \mathrm{SS}_{\mathrm{tot}}$.

Parameters:

Note: std::invalid_argument on a size mismatch among $x, y, \mathrm{weights}$, fewer than 2 points, or a zero weight sum. Use it when observations have different reliabilities (e.g. variance weighting).

Polynomial Fitting

Fits a degree-$d$ polynomial $p(x) = a_0 + a_1 x + \dots + a_d x^d$ to observed points $(x_i, y_i)$ in the least-squares sense.

PolynomialFitResult<T>

template<concepts::OrderedField T>
struct PolynomialFitResult {
    std::vector<T> coefficients;   // polynomial coefficients [a0, a1, ..., aN] (ascending)
    T             residual;       // residual sum of squares Σ (y_i - p(x_i))²
};

Function reference

polynomialFit

template<concepts::OrderedField T>
PolynomialFitResult<T> polynomialFit(
    std::span<const T> x, std::span<const T> y, std::size_t degree);
template<concepts::OrderedField T>
PolynomialFitResult<T> polynomialFit(
    const std::vector<T>& x, const std::vector<T>& y, std::size_t degree);

Behaviour: builds the normal equations $A c = b$ ($A_{ij} = \sum_k x_k^{i+j}$, $b_i = \sum_k x_k^i y_k$) and solves them by direct LU. For numerical stability it fits in a shifted space after subtracting the means of $x, y$, then converts the coefficients back to the original space via the binomial expansion. The residual sum of squares is returned alongside.

Parameters:

Note: std::invalid_argument on a size mismatch or fewer than $\mathrm{degree} + 1$ data points. Raising the degree worsens the condition number of the normal equations (a Vandermonde effect), so when a high degree is needed consider approximation via ChebyshevApprox.

// Example: fit a parabola y = x^2 (degree=2)
std::vector<double> x = {-2, -1, 0, 1, 2};
std::vector<double> y = { 4,  1, 0, 1, 4};
auto r = sangi::polynomialFit(x, y, 2);
// r.coefficients ≈ {0, 0, 1}  (a0=0, a1=0, a2=1 → p(x)=x²)
// r.residual ≈ 0

Chebyshev Approximation

Expands a known function $f$ on an interval $[a, b]$ into a degree-$n$ Chebyshev polynomial series $f(x) \approx \tfrac{c_0}{2} + \sum_{k=1}^{n} c_k\, T_k(y)$ (where $y$ is the linear map $[a,b] \to [-1,1]$). Because the function is sampled at the Chebyshev nodes (extrema) and the coefficients are obtained by a DCT, the error is distributed evenly across the whole interval — making it a near-minimax approximation among polynomial approximations. The maximum error is far smaller than a plain truncated Taylor series.

class ChebyshevApprox<T>

template<concepts::OrderedField T>
class ChebyshevApprox {
public:
    ChebyshevApprox(const std::function<T(T)>& f, T a, T b, std::size_t n);

    T operator()(T x) const;                       // evaluate with all terms
    T evaluate(T x, std::size_t m) const;          // evaluate with the first m terms only (truncated)

    const std::vector<T>& coefficients() const;     // Chebyshev coefficients c_k
    std::size_t          size() const;             // number of terms (= n + 1)
    T                    lower() const;            // interval lower bound a
    T                    upper() const;            // interval upper bound b

    std::vector<T> toPolynomialCoefficients() const; // convert to ordinary polynomial coefficients [a0,...,aN]
};

Constructor ChebyshevApprox(f, a, b, n)

Behaviour: at construction, evaluates $f$ at the Chebyshev nodes $y_k = \cos\!\big(\pi (k + \tfrac12)/N\big)$ ($N = n + 1$) and computes the Chebyshev coefficients $c_0, \dots, c_n$ in one pass via DCT-II. Subsequent evaluations are fast because they reuse the coefficients.

Parameters:

Note: std::invalid_argument if $a \ge b$.

operator() / evaluate

Behaviour: operator()(x) evaluates all $n + 1$ terms with the Clenshaw recurrence. evaluate(x, m) performs a truncated evaluation using only the first $m$ terms ($m \le n + 1$), trading accuracy for speed (returns $0$ when $m = 0$). The evaluation point $x$ is assumed to lie within $[a, b]$ (extrapolation makes the error grow rapidly).

Parameters:

toPolynomialCoefficients

Behaviour: builds ordinary monomial-basis polynomial coefficients $[a_0, a_1, \dots, a_n]$ ($p(x) = \sum_j a_j x^j$, ascending order) from the Chebyshev coefficients $c_k$. The $[a, b] \to [-1, 1]$ change of variable is folded into the expansion, so the result is a polynomial in the original $x$. Use it when you want to hand the result to other modules (Polynomial, root finding, etc.).

Note: at high degree the Chebyshev-to-monomial conversion worsens the condition number, so if you only need to evaluate, operator() (Clenshaw) is more numerically stable.

// Example: 8th-degree Chebyshev approximation of exp(x) on [-1, 1]
sangi::ChebyshevApprox<double> approx(
    [](double x){ return std::exp(x); }, -1.0, 1.0, 8);
double v = approx(0.5);             // ≈ exp(0.5) = 1.6487212707...
std::cout << v;
// Result: 1.64872127  (matches the true exp(0.5) to high precision)

Padé Approximation

Builds the $[M/N]$ Padé approximation $p(x)/q(x)$ ($\deg p = M$, $\deg q = N$, normalized to $q(0) = 1$) from a Taylor coefficient sequence. It extends accuracy beyond the radius of convergence farther than a Taylor truncation of the same number of terms, and is strong at approximating functions with poles.

PadeResult<T>

template<concepts::OrderedField T>
struct PadeResult {
    std::vector<T> numerator;     // numerator coefficients [p0, p1, ..., pM] (ascending)
    std::vector<T> denominator;   // denominator coefficients [1, q1, ..., qN]  (normalized to q0 = 1)
    bool          valid;         // whether a solution was obtained
};

Function reference

padeApprox

template<concepts::OrderedField T>
PadeResult<T> padeApprox(std::span<const T> taylor, std::size_t M, std::size_t N);
template<concepts::OrderedField T>
PadeResult<T> padeApprox(const std::vector<T>& taylor, std::size_t M, std::size_t N);

Behaviour: builds the $[M/N]$ Padé from Taylor coefficients $a_0, a_1, \dots, a_{M+N}$. It first solves for the denominator coefficients $q_1, \dots, q_N$ via an $N \times N$ linear system (direct LU), then computes the numerator $p_m = a_m + \sum_{l=1}^{\min(m,N)} q_l\, a_{m-l}$. When $N = 0$, the denominator is $1$ and the Taylor coefficients become the numerator directly.

Parameters:

Note: std::invalid_argument if fewer than $M + N + 1$ coefficients are given. Returns valid = false when the linear system cannot be solved (ill-conditioned).

// Example: [2/2] Padé from the Taylor series of exp(x)
std::vector<double> taylor = {1.0, 1.0, 0.5, 1.0/6, 1.0/24};  // 1, x, x²/2, ...
auto pa = sangi::padeApprox(taylor, 2, 2);
double v = sangi::evaluatePade(pa, 0.5);   // ≈ exp(0.5)
std::cout << v;
// Result: 1.64872...  (close to the true exp(0.5)=1.6487212707)

evaluatePade

template<concepts::OrderedField T>
T evaluatePade(const PadeResult<T>& result, T x);

Behaviour: for the result of padeApprox, evaluates the numerator and denominator each by Horner's method and returns $p(x)/q(x)$.

Parameters:

padeTable

template<concepts::OrderedField T>
std::vector<std::vector<PadeResult<T>>>
padeTable(std::span<const T> taylor, std::size_t Mmax, std::size_t Nmax);
template<concepts::OrderedField T>
std::vector<std::vector<PadeResult<T>>>
padeTable(const std::vector<T>& taylor, std::size_t Mmax, std::size_t Nmax);

Behaviour: computes the $[M/N]$ Padé for every $(M, N)$ with $0 \le M \le \mathrm{Mmax}$ and $0 \le N \le \mathrm{Nmax}$, returning a two-dimensional table indexed by row $M$ and column $N$ ($\mathrm{table}[M][N]$). Use it to compare the diagonal $[M/M]$ or anti-diagonal $[M + N = \text{const}]$ entries and pick the best approximation. Degenerate cells are stored as-is with valid = false.

Parameters:

Note: std::invalid_argument if there are too few coefficients. Since each cell solves an $N \times N$ system, the cost is $O(\mathrm{Mmax} \cdot \mathrm{Nmax} \cdot N^2)$. If you only want the diagonal, calling padeApprox directly is cheaper.

matrixExpPade

template<typename T>
Matrix<T> matrixExpPade(const Matrix<T>& A, int m = 6);

Behaviour: computes the matrix exponential $e^A$ by scaling-and-squaring plus a diagonal Padé kernel (Higham 2005). It is a thin wrapper over the library's matrix-exponential routine algorithms::expm; the Padé degree is selected automatically internally.

Parameters:

Note: std::invalid_argument for a non-square matrix. The argument m does not affect accuracy.

Continued Fractions (qd algorithm)

Converts a Taylor series $f(z) = c_0 + c_1 z + c_2 z^2 + \dots$ into a continued fraction. The coefficients are obtained by Rutishauser's qd (quotient-difference) algorithm. A continued fraction is equivalent to the Padé diagonal of the same number of terms; since evaluation can be truncated stage by stage, it is well suited to stable numerical evaluation of function values. Two forms are provided: the Stieltjes type (S-fraction) and the Jacobi type (J-fraction).

StieltjesFraction<T> / JacobiFraction<T>

// Stieltjes continued fraction:  f(z) = b0 / (1 + a1 z / (1 + a2 z / (1 + a3 z / ...)))
template<concepts::OrderedField T>
struct StieltjesFraction {
    T             b0;          // c_0 (leading constant, usually taylor[0])
    std::vector<T> a;          // a_1, a_2, a_3, ...
    bool          terminated;  // whether the qd table hit 0 and completed at finite length
    int           order;       // number of input Taylor coefficients - 1
};

// Jacobi continued fraction:  f(z) = c0 / (1 - β0 z - α1² z² / (1 - β1 z - α2² z² / ...))
template<concepts::OrderedField T>
struct JacobiFraction {
    T             c0;          // leading constant
    std::vector<T> beta;       // β_0, β_1, β_2, ...
    std::vector<T> alpha2;     // α_1², α_2², ...
    bool          terminated;
    int           order;
};

Function reference

taylorToSFraction / taylorToJFraction

template<concepts::OrderedField T>
StieltjesFraction<T> taylorToSFraction(std::span<const T> taylor, int kmax = -1);
template<concepts::OrderedField T>
StieltjesFraction<T> taylorToSFraction(const std::vector<T>& taylor, int kmax = -1);

template<concepts::OrderedField T>
JacobiFraction<T> taylorToJFraction(std::span<const T> taylor, int kmax = -1);
template<concepts::OrderedField T>
JacobiFraction<T> taylorToJFraction(const std::vector<T>& taylor, int kmax = -1);

Behaviour: runs the Taylor coefficient sequence through the qd algorithm and produces the S-fraction coefficients $a_k$ or the J-fraction coefficients $\beta_k, \alpha_k^2$. It truncates at the point where the qd table hits $0$ (Padé degeneracy) and sets terminated = true. The J-fraction is a form that merges two S-fraction stages into one, corresponding to the three-term recurrence of orthogonal polynomials.

Parameters:

Note: std::invalid_argument if the input is empty. When $c_0 = 0$ (leading term zero) the continued fraction is undefined, so it returns empty with terminated = true (the caller must shift the variable or similar).

evalSFraction / evalJFraction

template<concepts::OrderedField T>
T evalSFraction(const StieltjesFraction<T>& sf, T x);
template<concepts::OrderedField T>
T evalJFraction(const JacobiFraction<T>& jf, T x);

Behaviour: evaluates the continued fraction by building up from the bottom stage (bottom-up). evalSFraction computes $b_0 / (1 + a_1 x / (1 + a_2 x / \dots))$ and evalJFraction computes $c_0 / (1 - \beta_0 x - \alpha_1^2 x^2 / (1 - \beta_1 x - \dots))$. If the coefficient sequence is empty, the leading constant ($b_0$ / $c_0$) is returned directly.

Parameters:

// Example: convert the Taylor series of log(1+x) to an S-fraction and evaluate
std::vector<double> taylor = {0.0, 1.0, -0.5, 1.0/3, -0.25, 0.2};  // log(1+x)
auto sf = sangi::taylorToSFraction(taylor);
double v = sangi::evalSFraction(sf, 1.0);   // ≈ log(2) = 0.6931...
std::cout << v;
// Result: 0.6931...  (stable over a wider range of x than the truncated Taylor series)

Rational Interpolation

A multipoint Padé (Cauchy / Newton-Padé type) that fits a rational function $p(x)/q(x)$ to values $y_i$ at distinct sample points $x_i$. Whereas padeApprox is built from Taylor coefficients at a single point, rationalInterpolation is built from multiple sample points $(x_i, y_i)$ — that is the difference (if you need the Taylor-coefficient version, use padeApprox).

RationalInterpResult<T>

template<concepts::OrderedField T>
struct RationalInterpResult {
    std::vector<T> numerator;     // [p0, p1, ..., pM] (ascending)
    std::vector<T> denominator;   // [1, q1, ..., qN]  (normalized to q0 = 1)
    bool          valid;
};

Function reference

rationalInterpolation

template<concepts::OrderedField T>
RationalInterpResult<T> rationalInterpolation(
    const std::vector<T>& x, const std::vector<T>& y,
    std::size_t M, std::size_t N);

Behaviour: builds the linear system ($M + N + 1$ unknowns) with $p(x_i) - y_i\, q(x_i) = 0$ at each sample point and solves for the numerator and denominator coefficients at once by direct LU. The $[M/N]$ rational function satisfies $f(x_i) = y_i$ (with $M + N + 1$ equal to the number of data points). In floating point the interpolation conditions are not exact, so check the residual on the caller side if needed.

Parameters:

Note: std::invalid_argument when the sizes of $x$ and $y$ differ, or when the number of points is not $M + N + 1$. Returns valid = false when the linear system cannot be solved. Reference: Stoer & Bulirsch, "Introduction to Numerical Analysis" §2.2.4.

evalRational

template<concepts::OrderedField T>
T evalRational(const RationalInterpResult<T>& r, T x);

Behaviour: evaluates the numerator and denominator each by Horner's method and returns $p(x)/q(x)$.

Parameters:

// Example: [2/2] rational interpolation from 5 points
std::vector<double> x = {-2, -1, 0, 1, 2};
std::vector<double> y = {/* f(x_i) */};
auto r = sangi::rationalInterpolation(x, y, 2, 2);
if (r.valid) {
    double v = sangi::evalRational(r, 0.5);
    std::cout << v;
}

Example

#include <math/approx/approximation.hpp>
#include <cmath>
#include <iostream>
#include <vector>
using namespace sangi;

int main() {
    // (1) Chebyshev approximation: degree-10 approximation of cos(x) on [0, π], then evaluate
    ChebyshevApprox<double> cheb(
        [](double x){ return std::cos(x); }, 0.0, std::acos(-1.0), 10);
    std::cout << "cos(1.0) ≈ " << cheb(1.0) << '\n';
    // Output: cos(1.0) ≈ 0.540302...  (true value cos(1) = 0.5403023059)

    // (2) Padé approximation: [3/3] rational approximation from the Taylor series of exp(x)
    std::vector<double> taylor = {
        1.0, 1.0, 0.5, 1.0/6, 1.0/24, 1.0/120, 1.0/720
    };  // Maclaurin coefficients of exp(x) (7 = 3+3+1)
    auto pade = padeApprox(taylor, 3, 3);
    if (pade.valid) {
        std::cout << "exp(0.5) ≈ " << evaluatePade(pade, 0.5) << '\n';
        // Output: exp(0.5) ≈ 1.648721...  (true value 1.6487212707)
    }

    // (3) Linear regression: fit to a line
    std::vector<double> xs = {0, 1, 2, 3, 4};
    std::vector<double> ys = {1.1, 2.9, 5.2, 6.8, 9.1};
    auto fit = linearRegression(xs, ys);
    std::cout << "slope=" << fit.slope
              << " intercept=" << fit.intercept
              << " R²=" << fit.r_squared << '\n';
    // Sample output: slope≈2.0 intercept≈1.0 R²≈0.99
}

Related Mathematical Background

The following articles explain the mathematical concepts underlying the Approximation module.

• Approximation Theory — The framework of function approximation and best approximation
• Least Squares — Foundation of linear regression and polynomial fitting
• Chebyshev Approximation — Near-minimax approximation via Chebyshev series
• Minimax Approximation — Uniform-norm minimization and the equioscillation theorem
• Padé Approximation — Rational approximation and continued fractions