Memory-free multivariate polynomials

Package Announcements
,
20

SpectralKit.jl is not really intended to work with coefficients of the polynomials directly (it uses Chebyshev bases). It is possible but you have to iterate. Its intended uses case is getting a linear combination with a given coefficients, allocation-free. Eg

using SpectralKit, StaticArrays
dom = BoundedLinear(0, 200)
basis = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(4, 2), Val(2)) ∘
    coordinate_transformations(dom, dom)
θ = randn(dimension(basis));     # random coefficients
linear_combination(basis, θ, (1, 1))

The last line should not allocate when inside a function. You can iterate the last argument on CartesianIndices, as in (continuing the example)

julia> using BenchmarkTools

julia> cids = CartesianIndices((200, 200))
CartesianIndices((200, 200))

julia> buffer = Vector{Float64}(undef, length(cids));

julia> function eval_at_θ!(buffer, basis, θ, cids)
       for (i, x) in enumerate(cids)
       buffer[i] = linear_combination(basis, θ, Tuple(x))
       end
       buffer
       end
eval_at_θ! (generic function with 2 methods)

julia> @ballocated eval_at_θ!($buffer, $basis, $θ, $cids)
0

Note that I am never collecting the CartesianIndices, since they iterate just fine online.

One caveat: using Smolyak basis as tensor basis (as above) is not ideal and is not what this package has been optimized for. If there is a demand, I could introduce a TensorBasis, it is much simple as it does not have to jump through all the hoops of Smolyak indexing.

1 Like