CaratheodoryFejerApprox.jl

Build Status Coverage

Introduction

CaratheodoryFejerApprox.jl offers robust near-minimax approximation of arbitrary smooth functions using the Carathéodory-Fejér method. This method approximates real functions with polynomials or rationals by transplanting the problem onto the unit disk in the complex plane, applying the Carathéodory-Fejér theorem, and constructing a near-best approximation from the eigenvalues and eigenvectors of a Hankel matrix containing coefficients from the Chebyshev series expansion of the function.

Much of the functionality in CaratheodoryFejerApprox.jl began as a translation of Matlab code from the fantastic Chebfun package, specifically chebfun/@chebfun/cf.m. Additionally, internally we make heavy use of the analagous Julia package ApproxFun.jl.

Package Features

Polynomial Approximation

polynomialcf(f, m::Int) -> RationalApproximant{Float64}
polynomialcf(f, dom::NTuple{2, T}, m::Int) -> RationalApproximant{T}

Approximate a function f with a degree m polynomial CF approximant on the interval dom. If not specified, dom defaults to (-1.0, 1.0).

Rational Approximation

rationalcf(f, m::Int, n::Int) -> RationalApproximant{Float64}
rationalcf(f, dom::NTuple{2, T}, m::Int, n::Int) -> RationalApproximant{T}

Approximate a function f with a type (m, n) rational CF approximant on the interval dom, where m is the numerator degree and n is the denominator degree. If not specified, dom defaults to (-1.0, 1.0).

Minimax Fine-Tuning

minimax(f, m::Int, n::Int) -> RationalApproximant{Float64}
minimax(f, dom::NTuple{2, T}, m::Int, n::Int) -> RationalApproximant{T}

Compute the type (m, n) CF approximant and then, if necessary, fine-tune the approximant to become a true minimax approximant using the Remez algorithm. If not specified, dom defaults to (-1.0, 1.0).

Polynomial Coefficient Basis

chebcoeffs(res::RationalApproximant{T}) -> NTuple{2, Vector{T}}
monocoeffs(res::RationalApproximant{T}; transplant = true) -> NTuple{2, Vector{T}}
monocoeffs(res::RationalApproximant{T1}, ::Type{T2} = BigFloat; transplant = true) -> NTuple{2, Vector{T1}}

Extract polynomial coefficients in the monomial basis via monocoeffs(res) or in the Chebyshev basis via chebcoeffs(res). As converting to monomial coefficients can be numerically unstable, optionally pass a higher precision type to monocoeffs for intermediate computations.

When transplant = true (the default), the monomial coefficients correspond to the original function f(x) on the interval dom. If transplant = false, they correspond to the linearly transplanted function g(t) = f((x - mid) / rad) where mid and rad are the midpoint and radius of dom and -1 <= t <= 1.

Note that, particularly when |mid| is large, it can be much more numerically stable to evaluate the approximant via evalpoly(t, p) using the non-transplanted coefficients p and t = (x - mid) / rad.

The Chebyshev coefficients always correspond to the linearly transplanted function g(t) = f((x - mid) / rad) used internally; they are not transplanted to dom.

Usage

julia> using CaratheodoryFejerApprox

Polynomial approximant

Compute a degree 5 polynomial CF approximant to cos(x) on the interval [-0.5, 0.5]:

julia> res = polynomialcf(cos, (-0.5, 0.5), 5) 
RationalApproximant{Float64}
  Type:   m / n = 4 / 0
  Domain: -0.5 ≤ x ≤ 0.5
  Error:  |f - p| ⪅ 6.721e-7
  Approximant:
      p(t) = 0.9385⋅T_0(t) - 0.06121⋅T_2(t) + 0.0003215⋅T_4(t)
  where: t = (x - 0.0) / 0.5

A few things to note:

  • While we requested a degree 5 approximant, polynomialcf automatically recognized that cos(x) is an even function and therefore truncated the approximant to degree 4
  • The polynomial approximant is displayed in the Chebyshev basis, transplanted to the standard interval [-1, 1]. In this case, we see that cos(x) ≈ p(t) where t = 2x
  • Despite being low degree the approximant is quite accurate, with an estimated 6-7 digits of worst-case accuracy

Next, let's extract the coefficients of the rational approximant in the monomial basis:

julia> p, q = monocoeffs(res) 
([0.9999993278622336, 0.0, -0.49995153387633173, 0.0, 0.04114863415981116], [1.0])

The resulting monomial coefficients are close to (but not equal to) the Taylor expansion of cos(x) at the origin: cos(x) = 1 - x^2/2 + x^4/24 + O(x^6) ≈ 1.0 - 0.5 x^2 + 0.04166 x^4. This is a common feature of minimax approximants; they are often similar to Taylor series expansions, but adjusted to tradeoff accuracy near the Taylor expansion point for accuracy over the whole interval.

Rational approximant

We can simlarly compute a type (4, 4) rational approximant to exp(x) on [-1, 1] as follows:

julia> res = rationalcf(exp, 4, 4) 
RationalApproximant{Float64}
  Type:   m / n = 4 / 4
  Domain: -1.0 ≤ x ≤ 1.0
  Error:  |f - p / q| ⪅ 1.538e-10
  Approximant:
      p(x) = 1.054⋅T_0(x) + 0.511⋅T_1(x) + 0.05434⋅T_2(x) + 0.003018⋅T_3(x) + 7.582e-5⋅T_4(x)
      q(x) = 1.053⋅T_0(x) - 0.5068⋅T_1(x) + 0.05333⋅T_2(x) - 0.002919⋅T_3(x) + 7.173e-5⋅T_4(x)

Note again the low error - approximately 10 digits of accuracy - despite the small degrees of the numerator and denominator.

Minimax approximant

Now compare the type (4, 4) CF approximant with the corresponding type (4, 4) minimax approximant:

julia> res = minimax(exp, 4, 4) 
RationalApproximant{Float64}
  Type:   m / n = 4 / 4
  Domain: -1.0 ≤ x ≤ 1.0
  Error:  |f - p / q| ⪅ 1.538e-10
  Approximant:
      p(x) = 1.054⋅T_0(x) + 0.511⋅T_1(x) + 0.05434⋅T_2(x) + 0.003018⋅T_3(x) + 7.582e-5⋅T_4(x)
      q(x) = 1.053⋅T_0(x) - 0.5068⋅T_1(x) + 0.05333⋅T_2(x) - 0.002919⋅T_3(x) + 7.173e-5⋅T_4(x)

We see that the coefficients are identical, and in fact the error bound is sharp over the whole interval. This is the magic of CF approximants: they are often very nearly true minimax approximants.

High-precision rational approximant

Care has been taken to ensure all internal computations work for generic float types, so we can easily compute a high degree minimax approximant in arbitrary precision:

julia> res = minimax(sin, BigFloat.((-1, 1)), 20, 16) 
RationalApproximant{BigFloat}
  Type:   m / n = 19 / 16
  Domain: -1.0 ≤ x ≤ 1.0
  Error:  |f - p / q| ⪅ 8.282e-56
  Approximant:
      p(x) = 0.8859⋅T_1(x) - 0.03733⋅T_3(x) + 0.0004165⋅T_5(x) - 1.997e-6⋅T_7(x) + 4.995e-9⋅T_9(x) - 7.17e-12⋅T_11(x) + 6.155e-15⋅T_13(x) - 3.147e-18⋅T_15(x) + 8.94e-22⋅T_17(x) - 1.102e-25⋅T_19(x)
      q(x) = 1.004⋅T_0(x) + 0.00447⋅T_2(x) + 4.94e-6⋅T_4(x) + 3.57e-9⋅T_6(x) + 1.861e-12⋅T_8(x) + 7.234e-16⋅T_10(x) + 2.072e-19⋅T_12(x) + 4.0509999999999997e-23⋅T_14(x) + 4.201e-27⋅T_16(x)

Here we have computed a type (20, 16) minimax approximation to sin(x) on [-1, 1] accurate to about 55 digits. Note again that internally sin(x) was recognized to be an odd function of x, and therefore the numerator degree was truncated from 20 to 19.

Application: implementing the modified Bessel function of the first kind of order zero

To show how minimax approximants may be used in practice, let's derive polynomial approximants for the modified Bessel function of the first kind of order zero I₀(x) using polynomialcf. We'll compare the approximants to those provided by the Bessels.jl package.

Following Bessels.jl, we'll use the following piecewise approximant for I₀(x):

  • For x < 7.75 we let I₀(x) = 1 + (x/2)^2 P₁((x/2)^2)
  • For x ≥ 7.75 we let I₀(x) = exp(x) P₂(1/x) / sqrt(x)

where P₁ and P₂ are polynomial approximants. Note that I₀(x) is an even function, so we need only consider positive x.

First, we need the true I₀(x). We will compute this via the integral representation I₀(x) = (1 / π) ∫_0^π exp(x cos(t)) dt, which can be evaluated precisely using QuadGK.jl with BigFloats:

julia> using QuadGK, Bessels 

julia> function besseli0_quadgk(x::BigFloat; expscale = false) 
           I, E = quadgk(BigFloat(0.0), BigFloat(π); order = 21, rtol = 1e-50) do t
               # I₀(x) diverges exponentially for large x; optionally compute exp(-x) I₀(x) instead
               return expscale ? exp(x * (cos(t) - 1)) : exp(x * cos(t))
           end
           return I / π
       end; 

julia> besseli0_quadgk(x; kwargs...) = oftype(x, besseli0_quadgk(BigFloat(x); kwargs...));

Let's first check that this matches Bessels.jls implementation:

julia> @assert Bessels.besseli0(0.5) ≈ besseli0_quadgk(0.5; expscale = false) 

julia> @assert Bessels.besseli0x(100.0) ≈ besseli0_quadgk(100.0; expscale = true)

Now, we use polynomialcf to compute the polynomial P₁ for x < 7.75:

# We are aiming to find P₁(t) such that I₀(x) = 1 + (x/2)^2 P₁((x/2)^2)
julia> xdom = Double64.((0.0, 7.75)); # compute coefficients in higher precision 

julia> tdom = (xdom ./ 2) .^ 2; # t = (x/2)^2 

julia> res = polynomialcf(tdom, 13) do t 
           x = 2√t
           I₀ = besseli0_quadgk(x; expscale = false)
           return t == 0 ? one(t) : (I₀ - 1) / t # P₁(t) = (I₀(x) - 1) / (x/2)^2 = (I₀(x) - 1) / t
       end 
RationalApproximant{Double64}
  Type:   m / n = 13 / 0
  Domain: 0.0 ≤ x ≤ 15.02
  Error:  |f - p| ⪅ 2.0e-16
  Approximant:
      p(t) = 8.518⋅T_0(t) + 10.13⋅T_1(t) + 3.143⋅T_2(t) + 0.5982⋅T_3(t) + 0.0769⋅T_4(t) + 0.007112⋅T_5(t) + 0.0004954⋅T_6(t) + 2.689e-5⋅T_7(t) + 1.169e-6⋅T_8(t) + 4.158e-8⋅T_9(t) + 1.232e-9⋅T_10(t) + 3.087e-11⋅T_11(t) + 6.625e-13⋅T_12(t) + 1.231e-14⋅T_13(t)
  where: t = (x - 7.508) / 7.508

We see that the approximant is accurate to about eps(Float64) over the interval. Now, let's compare these coefficients to those from Bessels.jl:

julia> const P₁_vec = monocoeffs(res)[1] .|> Float64; # coefficients for P₁(t) 

julia> const P₁_tup = (P₁_vec...,); # `evalpoly` is faster with tuples of coefficients 

julia> maximum(abs, (P₁_tup .- Bessels.besseli0_small_coefs(Float64)) ./ P₁_tup) # maximum relative difference 
0.0

Our derived coefficients are identical. Similarly, let's now use polynomialcf to compute the polynomial approximant P₂(t) for x ≥ 7.75:

# We are aiming to find P₂(t) such that I₀(x) = exp(x) P₂(1/x) / sqrt(x)`
julia> tdom = Double64.((1 / 1e6, 1 / 7.75)); # t = 1/x; domain bounds copied from `Bessels.jl` 

julia> res = polynomialcf(tdom, 21) do t 
           x = inv(t)
           I₀ₓ = besseli0_quadgk(x; expscale = true) # exp(-x) * I₀(x)
           return √x * I₀ₓ # P₂(t) = √x * exp(-x) * I₀(x)
       end 
RationalApproximant{Double64}
  Type:   m / n = 21 / 0
  Domain: 1.0e-6 ≤ x ≤ 0.129
  Error:  |f - p| ⪅ 6.292e-17
  Approximant:
      p(t) = 0.4024⋅T_0(t) + 0.003488⋅T_1(t) + 7.408e-5⋅T_2(t) + 3.25e-6⋅T_3(t) + 2.43e-7⋅T_4(t) + 2.837e-8⋅T_5(t) + 4.332e-9⋅T_6(t) + 5.778e-10⋅T_7(t) - 2.278e-11⋅T_8(t) - 5.335e-11⋅T_9(t) - 1.818e-11⋅T_10(t) - 1.24e-12⋅T_11(t) + 1.425e-12⋅T_12(t) + 5.084e-13⋅T_13(t) - 4.223e-14⋅T_14(t) - 6.879e-14⋅T_15(t) - 6.505e-15⋅T_16(t) + 8.044e-15⋅T_17(t) + 1.82e-15⋅T_18(t) - 9.757e-16⋅T_19(t) - 3.463e-16⋅T_20(t) + 1.191e-16⋅T_21(t)
  where: t = (x - 0.06452) / 0.06452

The approximant is again accurate to about eps(Float64) over the interval. Let's again compare these coefficients to those from Bessels.jl:

julia> const P₂_vec = monocoeffs(res)[1] .|> Float64; # coefficients for P₂(t) 

julia> const P₂_tup = (P₂_vec...,); # `evalpoly` is faster with tuples of coefficients 

julia> maximum(abs, (P₂_tup .- Bessels.besseli0_med_coefs(Float64)) ./ P₂_tup) # maximum relative difference 
1.1172124075993363e-15

Again, we see that the derived coefficients are in good agreement.

Finally, we can now define our implementation of I₀(x) using the approximants P₁ and P₂. Modifying the code from Bessels.jl, we have:

julia> function besseli0_cf(x::Float64) 
           x = abs(x)
           if x < 7.75
               a = x * x / 4
               return muladd(a, evalpoly(a, P₁_tup), 1)
           else
               a = exp(x / 2)
               s = a * evalpoly(inv(x), P₂_tup) / sqrt(x)
               return a * s
           end
       end; 

julia> @assert Bessels.besseli0(0.5) ≈ besseli0_cf(0.5) 

julia> @assert Bessels.besseli0(100.0) ≈ besseli0_cf(100.0)

References

[1] Real Polynomial Chebyshev Approximation by the Carathéodory-Fejér Method (Gutknecht and Trefethen, 1982)

This foundational work introduces a novel method for near-best approximation of real functions by polynomials. The approach involves transforming the problem onto the unit disk and applying the Carathéodory-Fejér theorem. The resulting approximation is constructed from the principal eigenvalue and eigenvector of a Hankel matrix of Chebyshev coefficients. The method offers high-order agreement with the best approximation, making it of both practical and theoretical importance. This package implements this work in polynomialcf.

[2] The Carathéodory-Fejér Method for Real Rational Approximation (Trefethen and Gutknecht, 1983)

This work extends the Carathéodory-Fejér method for polynomial minimax approximants to real rational approximation, similarly leveraging eigenvalue analysis of a Hankel matrix of Chebyshev coefficients to achieve this. The rational CF approximants frequently approach the best rational approximation with high accuracy. This package implements this work in rationalcf.

[3] A Robust Implementation of the Carathéodory-Fejér Method for Rational Approximation (Van Deun and Trefethen, 2011)

This work details a robust implementation of the Carathéodory-Fejér method for both polynomial and rational approximation in Matlab within the Chebfun package ecosystem. CaratheodoryFejerApprox.jl is based largely on this reference implementation, providing users with a powerful tool for rational approximation in Julia which additionally works for generic float types beyond Float64, such as BigFloat and Double64.

This page was generated using Literate.jl.