[ANN] ArnoldiVandermonde.jl

Package Announcements
1

I am pleased to announce version 0.1.0 of ArnoldiVandermonde.jl, a lightweight package to numerically orthogonalize polynomials on discrete point sets in the complex plane.

A monomial basis 1, x, x^2, \ldots is disastrously ill-conditioned unless your domain is the unit circle in \mathbf{C}. The Chebyshev (and other Jacobi) polynomials are famously well-conditioned for real intervals, but not as useful on most other domains. The Arnoldi–Vandermonde iteration uses the Arnoldi iteration (i.e., modified Gram–Schmidt) to numerically orthogonalize a polynomial basis as it is grown. This basis can be used to give stable, spectrally accurate L_2 projections of smooth functions. See Brubeck & Trefethen (2021) for an introduction to the algorithm.

For instance:

julia> using ArnoldiVandermonde
julia> f(x) = sin(exp(x));

julia> @time p = project(f, 0, 1)
  0.000195 seconds (684 allocations: 1.636 MiB)
Arnoldi polynomial of degree 17 on 65 nodes

julia> @time p = project(f, BigFloat(0), 1)
  7.304747 seconds (171.66 M allocations: 15.354 GiB, 10.02% gc time)
Arnoldi polynomial of degree 77 on 257 nodes

julia> maximum(abs(f(x) - p(x)) for x in range(BigFloat(0), 1, 10000)) |> Float64
3.0407129046474806e-72
2

Great! Seems like it fully resolves this thread: Julia analog to R's poly function, orthogonal polynomials by QR decomposition

How does it compare to the ArnoldiFit object in Polynomials.jl?

3

Good question!

f(x) = sin(exp(x));
x = range(0, 1, 800);

function test1(f, x, n)
   B = ArnoldiBasis(x, n)
   return B \ f
end

function test2(f, x, n)
   return fit(ArnoldiFit, x, f.(x), n)
end

julia> @benchmark test1($f, $x, 16)
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max):  114.875 ΞΌs …  22.270 ms  β”Š GC (min … max):  0.00% … 99.25%
 Time  (median):     135.208 ΞΌs               β”Š GC (median):     0.00%
 Time  (mean Β± Οƒ):   170.619 ΞΌs Β± 366.259 ΞΌs  β”Š GC (mean Β± Οƒ):  20.17% Β± 12.29%

  β–ˆβ–                                                             
  β–ˆβ–ˆβ–‚β–‚β–‚β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–‚β–β–β–β–β–β–β–‚β–β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚β–‚ β–‚
  115 ΞΌs           Histogram: frequency by time         1.54 ms <

 Memory estimate: 631.92 KiB, allocs estimate: 53.

julia> @benchmark test2($f, $x, 16)
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max):  164.667 ΞΌs …  20.761 ms  β”Š GC (min … max):  0.00% … 98.29%
 Time  (median):     340.375 ΞΌs               β”Š GC (median):     0.00%
 Time  (mean Β± Οƒ):   475.733 ΞΌs Β± 604.246 ΞΌs  β”Š GC (mean Β± Οƒ):  30.12% Β± 22.17%

  β–„β–ƒβ–„β–ˆβ–ˆβ–†β–ƒ                                                       β–‚
  β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–‡β–ƒβ–β–β–β–ƒβ–β–β–β–β–β–β–β–β–β–ƒβ–β–β–β–ƒβ–„β–„β–†β–…β–†β–†β–‡β–†β–‡β–‡β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–‡β–‡β–‡β–‡β–‡β–‡β–‡β–†β–‡β–‡β–†β–†β–†β–‡β–† β–ˆ
  165 ΞΌs        Histogram: log(frequency) by time       2.69 ms <

 Memory estimate: 3.16 MiB, allocs estimate: 1360.

For evaluation, I did a separate implementation for array inputs that pays off:

y = range(0, 1, 2000)
p1 = test1(f, x, 16)
p2 = test2(f, x, 16)

julia> @benchmark $p1.($y)
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max):  225.042 ΞΌs …  1.282 ms  β”Š GC (min … max): 0.00% … 80.42%
 Time  (median):     236.709 ΞΌs              β”Š GC (median):    0.00%
 Time  (mean Β± Οƒ):   240.896 ΞΌs Β± 28.685 ΞΌs  β”Š GC (mean Β± Οƒ):  0.38% Β±  2.70%

  β–β–ˆ                                                            
  β–ˆβ–ˆβ–„β–„β–„β–‡β–…β–ƒβ–„β–„β–ƒβ–ƒβ–ƒβ–„β–ƒβ–ƒβ–ƒβ–ƒβ–ƒβ–ƒβ–„β–ƒβ–ƒβ–ƒβ–ƒβ–„β–„β–‚β–ƒβ–‚β–‚β–‚β–‚β–‚β–β–β–‚β–‚β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β– β–‚
  225 ΞΌs          Histogram: frequency by time          288 ΞΌs <

 Memory estimate: 16.06 KiB, allocs estimate: 3.

julia> @benchmark $p1($y)
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max):  60.541 ΞΌs …   4.412 ms  β”Š GC (min … max):  0.00% … 97.42%
 Time  (median):     74.958 ΞΌs               β”Š GC (median):     0.00%
 Time  (mean Β± Οƒ):   92.949 ΞΌs Β± 183.744 ΞΌs  β”Š GC (mean Β± Οƒ):  18.08% Β±  9.45%

  β–ˆβ–ˆ                                                            
  β–ˆβ–ˆβ–ƒβ–‚β–‚β–‚β–‚β–‚β–‚β–β–‚β–‚β–β–β–β–β–β–β–β–β–β–‚β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–‚β–β–β–β–β–‚β–β–β–β–β–β–‚ β–‚
  60.5 ΞΌs         Histogram: frequency by time          907 ΞΌs <

 Memory estimate: 304.14 KiB, allocs estimate: 6.

julia> p2 = test2(f, x, 16)
ArnoldiFit of degree 16

julia> @benchmark $p2.($y)
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max):  174.583 ΞΌs …   4.799 ms  β”Š GC (min … max):  0.00% … 95.91%
 Time  (median):     181.458 ΞΌs               β”Š GC (median):     0.00%
 Time  (mean Β± Οƒ):   246.001 ΞΌs Β± 493.594 ΞΌs  β”Š GC (mean Β± Οƒ):  24.63% Β± 11.59%

  β–ˆ                                                             ▁
  β–ˆβ–ˆβ–ƒβ–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–β–„β–†β–…β–‡β–… β–ˆ
  175 ΞΌs        Histogram: log(frequency) by time       4.06 ms <

 Memory estimate: 391.06 KiB, allocs estimate: 4003.

I’ll have to add a mention of ArnoldiFit to the readme; I completely forgot about it. But I need the basis as well as the polynomial for some other stuff I’m working on.

4

The reference should be Brubeck, Nakatsukasa, and Trefethen (2021).