CoordinateSplittingPTrees
Note: this package is under active development and is subject to change.
Interpolation, accuracy, and the curse of dimensionality
Interpolation is the process of building a model from data available at discrete locations. Fulldegree multidimensional polynomial interpolation differs from most other interpolation schemes in that it is not, in general, continuous over space. Instead, the goal is to find a local polynomial fit of full degree and use this for detailed investigations of local structurein essence, prioritizing accuracy of local interpolation over any global property like continuity or smoothness.
An important application fulldegree interpolation is derivativefree optimization (DFO), where many algorithms work by constructing a quadratic fit to the local function values and then minimizing the quadratic over some trust regionsince the goal is to find the location of the minimum (or minima, when there are several), global continuity of the interpolation scheme is not an urgent priority. For use in DFO (and presumably most other interpolation applications), one therefore requires a rule to select a "complete" subset of evaluation positions to determine the model coefficients. If different point subsets are selected in different regions of space, the polynomial will change discontinuously as one moves around the space; but for purposes of optimization, this is an acceptable tradeoff that enables one to quickly find the exact minimum via quasiNewton methods.
Unfortunately, when n
is large the computational burden of building
the model can be high. For example, the number of parameters in a
fulldegree quadratic model in n
dimensions is (n+1)*(n+2)/2
.
This implies that data from (at least) an equal number of evaluation
points must be used; moreover, solving for the model via linear
regression is in the general case an O(n^6)
problem. This becomes
prohibitive for n
in the hundreds and extremely costly even for much
more modest problems. Powell
described
a method for updating a quadratic model with a new data point; this
method has a cost in the range O(n^2)O(n^4)
per update depending
on the amount of history maintained, and for local DFO this represents
a significant improvement. However, this method is less attractive in
the context of global optimization, where multiple regions in space
may be pursued simultaneously but for which one may not wish to
maintain storage of O(n^2)
parameters for many different local
models.
Reducing the burden with CSptrees
This repository implements a new (?) data structure, the coordinate
splitting tree of degree p or CSptrees for short. As shown
here, this data structure supports efficient and
accurate fulldegree multidimensional polynomial interpolation. For
example, a CS1tree reduces the burden of fitting a quadratic model to
~O(n^4)
, and a CS2tree reduces it to O(n^2)
. Moreover, the
resulting coefficients are determined with significantly higher
precision than by traditional approaches. CSp trees naturally
implement adaptive mesh refinement (AMR), as adding new evaluation
points to the interpolation scheme corresponds to adding new nodes to
the tree. This package is written in the
Julia programming language.
CSptrees have a close relationship with
kd trees, and indeed CS1
trees can be implemented essentially as a kd tree with only one
evaluation point (and substantial restrictions on its position) per
box. Specifically, a kd tree "splits" space along one dimension at a
time; a CS1 tree selects a new evaluation point that differs from its
"parent" position in only one coordinate, thus creating a new box that
splits the parent along this coordinate. A general CSp tree splits p
dimensions simultaneously, creating a set of 2^p1
new boxes each
associated with a unique position and new function evaluation.
For now, to learn how to use this pacakge just type
?CoordinateSplittingPTrees
at the Julia prompt, and then read the
help associated with each function. For convenience, that overall
summary helptext is reproduced below.
Usage summary
Note that more detail is available in the help for each function.
The world bounds are specified by creating a World
object (see ?World
).
The tree is initialized with root = Box{p}(world)
, and new leaves are
added with Box(parent, dims, xs, metas)
(see ?Box
).
You can interact with boxes using the following API (again, this is currently under development and may be unimplemented, outdated, or in flux):
Polynomial fits
polynomial_full
: construct a full polynomial of the specified degreepolynomial_minimal
: construct a minimal polynomial of the specified degree
Tree API

position
: retrieve the position of aBox

meta
: retrieve the metadata associated with aBox

value
: retrieve just the function value associated with aBox

boxbounds
: retrieve the edges of aBox

isroot
: determine whether a givenBox
is the rootnode 
isleaf
: determine whether a givenBox
is a leafnode 
getroot
: get the root node associated with a tree 
getleaf
: get the leaf node associated with a (possiblyparent) box 
find_leaf_at
: find the leaf node containing a particular position 
addpoint!
: add a new ndimensional evaluation point 
leaves
: returns an iterator for visiting all leafnodes
Display
print_tree
: display the tree structure (from AbstractTrees.jl)splitprint
: a more compact display of the tree structuresplitprint_colored
: similar tosplitprint
but with highlights