Proposal: Hierarchical Basis functions

[floswald]

hi there

I have been reading this paper recently: http://johannesbrumm.com/wp-content/uploads/2017/09/Brumm-Scheidegger-2017-ECTA.pdf and they use adaptive sparse grids to compute a huge economic model. I thought that would be very cool to have in julia.

• Here, I started a mini-prototype of their basis functions, which are just simple hat functions on [0,1]. they are called hierarchical basis functions because each consecutive order nests the previous orders.
• the main functionality is largely missing in this PR. this is just to get a discussion going and see if you are interested.
• I would definitely need help bringing this into the BasisMatrices.jl format.
• Here is all that this PR does:

julia> p = BasisMatrices.HierarcParams(-3,2,4)
Hierarchical Basis parameters from -3.0 to 2.0 
refinement level: 4
nodes:17

julia> p.nodes
17-element Array{Float64,1}:
 0.0   
 0.0625
 0.125 
 0.1875
 0.25  
 0.3125
 0.375 
 0.4375
 0.5   
 0.5625
 0.625 
 0.6875
 0.75  
 0.8125
 0.875 
 0.9375
 1.0   

julia> p.l
4

julia> BasisMatrices.evalbase(p,[0.0,0.01,0.02])
3×17 SparseMatrixCSC{Float64,Int64} with 5 stored entries:
  [1 ,  1]  =  1.0
  [2 ,  1]  =  0.84
  [3 ,  1]  =  0.68
  [2 ,  2]  =  0.16
  [3 ,  2]  =  0.32

What is so cool about this?

• there is a grid refinement step, that adds more points in "difficult" regions. (not in my code!)
• there is a fully fledged C++ library which does all of that, in fact I wrapped that last week for julia: https://github.com/floswald/Tasmanian.jl
• why have a julia version of this C++ library? no dependencies, greater flexibility. but also quite some work, so wanted so what you think.

cheers

[sglyon]

Hey @floswald this looks very cool and it something I'd love to see in Julia. Thanks for kicking the tires here

I've studied this paper before and written my own "toy" implementations in Julia before.

Do you have an idea for how much effort it would be to get a workable version in Julia? Also an estimate of how much of that work you think you'd be able to do would be super helpful.

[floswald]

I looked a bit into TASMANIAN, referenced above. This does a ton of stuff, so it makes no sense at all to redo all of this. So I would focus just on the simplest case for now, i.e. the proposed linear basis (hat) function. to get a simple version of the above paper, I think one would need

• a struct that carries current coefs and basis functions.
• a method that figures out the set of parent and offspring nodes (see page 14 in https://tasmanian.ornl.gov/documents/UserManual_5_1.pdf)
• a method setRefinement(g::ASG,tol) that will add points to the current g wherever the surplus is lager than tol. surplus is another name for coefficient. here, this is simply f(x) - \hat{f}_{l}(x), where l is the current level of refinement. in other words the surplus relates to the approximation error at point x. much better explained on page 14 of said pdf.
• extend this into multidimensions

It's low prioritiy for me after I wrapped Tasmanian, but I'd be happy to contribute a couple of hours! can you post your toy example here, maybe that would be useful?

[sglyon]

Ok, I think for now we will probably be content with using Tasmanian.jl (thanks for wrapping, but the way!)

I looked for my toy implementations of these routines and couldn't find them... Bad file management on my part. To be honest, I don't think we have lost too much