ANN: Smolyak.jl
Spencer Lyon
Chase Coleman and I (Spencer Lyon) have been working on some routines to do Smolyak interpolation/function approximation in Julia and Python. We are still in the early stages of the project, but have some promising results so far.
Smolyak Basics
Smolyak interpolation consists of two objects: a grid and a interpolating polynomial. Typically, the sparse grid is generating using nested sets of the extrema of the Chebychev polynomials. Similarly, the polynomial is constructed using nested sets of uni-dimensional basis polynomials (generally Chebychev polynomials of the first kind). There have been numerous applications of this procedure to economic models. A few examples are:
- Krueger and Kubler (2004)
- Malin, Krueger, and Kubler (2007)
Malin, Krueger, and Kubler (2011)
Gordon (2011)
To make the algorithm more efficient, Judd, et. al. (2013) suggest grid and interpolating polynomial construction routines that rely on disjoint sets of points and basis elements instead of the traditional nested sets.
Smolyak.jl
The novel idea from Judd, el. al. (2013) is to not repeat calculations or function evaluations by having the grid and polynomial to be constructed using unique elements only. We have taken this idea and improved upon it using some combinatorics tricks. We’ll be sure to document them soon (sorry I don’t do it here: it takes some background, which I don’t have time to lay right now). So far, this has enabled us to build grids with many more dimensions or much greater “density” (which typically translates to greater accuracy in interpolation) than is possible with classical nested sets.
At this stage, the project is very young and we have some more work to go, but wanted to make everyone aware of it and invite all who are interested to contribute. The code is on github under the EconForge organization and can be found here https://github.com/EconForge/Smolyak.jl.
As this is a bi-lingual project, the structure of the repository is a little unique. We try to follow conventions for both Julia and Python packages as explained in the README.
We will be updating the issue list on github soon, which should give a good idea of where the project currently stands and what still needs to happen. We also wanted to open this thread so we can have another location where people can give input.
Let us know if you are interested in using/contributing to the project!
Pablo Winant
The old smolyak code in dolo needed a refresh anyway since it was consuming to much memory. It's probably good to take a fresh start here. Have you thought about implementing linear interpolation, using hierarchical sparse grids too ? (as in Matlab's sparse grid interpolation toolbox ?) It has automatic dimension refinement and there maybe some inspiration to take there.
Pablo
Spencer Lyon
These are actually near the top of our todo list. I have updated the github issue list to reflect this. I will be filling in the details on the issues soon. A summary of the TODO list is:
- Develop a more efficient algorithm for constructing the matrix
(matrix used for interpolating polynomial)
- Implement anisotropic grids. This allows different
parameters (density of grid) along each dimension. Should be fairly simple
- Implement adaptive domain. Right now our grids are defined on
. We will have two methods of generalizing this. I will give details on the GH issue soon
- Implement linear interpolation. All the pieces are there for this, we just need to design a good interface and write a few functions.
I am not sure what you mean by hierarchical sparse grids. Can you elaborate?
Pablo Winant
These are actually near the top of our todo list. I have updated the github issue list to reflect this. I will be filling in the details on the issues soon. A summary of the TODO list is:
- Develop a more efficient algorithm for constructing the matrix
- Implement anisotropic grids. This allows different
- Implement adaptive domain. Right now our grids are defined on
- Implement linear interpolation. All the pieces are there for this, we just need to design a good interface and write a few functions.
I am not sure what you mean by hierarchical sparse grids. Can you elaborate?
One of the advantages is that the coefficient of the element on top, is a good indication of whether the the grid needs to be refined or not. The interpolant can also be represented as a tree that can be refined locally, when needed. The keyword for that line of research is "adaptive sparse grids" (not "dimension adaptive"). I think that Brumm and Scheidegger have applied that approach in economics.
On Thursday, November 14, 2013 3:54:05 AM UTC-5, Pablo Winant wrote:
Excellent ! That's a very good project, I'm impatient to see its progress.
The old smolyak code in dolo needed a refresh anyway since it was consuming to much memory. It's probably good to take a fresh start here. Have you thought about implementing linear interpolation, using hierarchical sparse grids too ? (as in Matlab's sparse grid interpolation toolbox ?) It has automatic dimension refinement and there maybe some inspiration to take there.
Pablo
Le mercredi 13 novembre 2013 19:37:52 UTC+1, Spencer Lyon a écrit :Chase Coleman and I (Spencer Lyon) have been working on some routines to do Smolyak interpolation/function approximation in Julia and Python. We are still in the early stages of the project, but have some promising results so far.
Smolyak Basics
Smolyak interpolation consists of two objects: a grid and a interpolating polynomial. Typically, the sparse grid is generating using nested sets of the extrema of the Chebychev polynomials. Similarly, the polynomial is constructed using nested sets of uni-dimensional basis polynomials (generally Chebychev polynomials of the first kind). There have been numerous applications of this procedure to economic models. A few examples are:
- Krueger and Kubler (2004)
- Malin, Krueger, and Kubler (2007)
Malin, Krueger, and Kubler (2011)
Gordon (2011)
To make the algorithm more efficient, Judd, et. al. (2013) suggest grid and interpolating polynomial construction routines that rely on disjoint sets of points and basis elements instead of the traditional nested sets.
Smolyak.jl
The novel idea from Judd, el. al. (2013) is to not repeat calculations or function evaluations by having the grid and polynomial to be constructed using unique elements only. We have taken this idea and improved upon it using some combinatorics tricks. We’ll be sure to document them soon (sorry I don’t do it here: it takes some background, which I don’t have time to lay right now). So far, this has enabled us to build grids with many more dimensions or much greater “density” (which typically translates to greater accuracy in interpolation) than is possible with classical nested sets.
At this stage, the project is very young and we have some more work to go, but wanted to make everyone aware of it and invite all who are interested to contribute. The code is on github under the EconForge organization and can be found here https://github.com/EconForge/Smolyak.jl.
As this is a bi-lingual project, the structure of the repository is a little unique. We try to follow conventions for both Julia and Python packages as explained in the README.
We will be updating the issue list on github soon, which should give a good idea of where the project currently stands and what still needs to happen. We also wanted to open this thread so we can have another location where people can give input.
Let us know if you are interested in using/contributing to the project!
--
You received this message because you are subscribed to the Google Groups "econforge" group.
To unsubscribe from this group and stop receiving emails from it, send an email to econforge+...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.
Sébastien Villemot
> They explain it well in "T. Gerstner, M. Griebel. Sparse Grids." It
> can be understood as a way to implement piecewise linear interpolation
> on a sparse grid. The elements of various orders are added on top of
> each other (hence the "hierarchical").
> One of the advantages is that the coefficient of the element on top,
> is a good indication of whether the the grid needs to be refined or
> not. The interpolant can also be represented as a tree that can be
> refined locally, when needed. The keyword for that line of research is
> "adaptive sparse grids" (not "dimension adaptive"). I think that
> Brumm and Scheidegger have applied that approach in economics.
>
> The disadvantage is that even after many successive refinements, the
> function may still not be smooth enough.
grid" (I attach a paper about this). This method seems very appealing
for finite elements, because the refinement steps *removes* basis
functions with a large support (while traditional hierarchical
refinement only adds new basis functions). Hence the system to be solved
becomes more sparse, because there is less overlap between basis
functions domains. I don't know however if anybody has used it in high
dimension (i.e. > 3D).
Anyways, your package is a great step forward. Thanks for doing it! I am
however a bit skeptical about mixing python and julia code in the same
package. My fear is that this may become a hindrance when trying to push
the package in both communities. Also, it will force synchronized
numbered releases between both languages, which is not necessarily what
you may want.
On my side, I plan to work on a julia package for solving systems of
nonlinear equations. It seems that no such package exists in pure julia:
there are packages that depend on external libs, and there is also an
Optim package (in pure julia) but that does not solve systems of
equations. Please let me know if you are aware of such an effort in the
julia ecosystem, to avoid redundant effort.
Best,
--
Sébastien Villemot
Researcher in Economics
Dynare developer
http://www.dynare.org/sebastien
Spencer Lyon
Sébastien Villemot
> Sébastien, I appreciate the concern with having a bi-lingual
> repository and see that this may cause some issues down the road.
> Right now our design goals are to develop similar, if not identical,
> libraries in both Julia and Python. We would like these libraries to
> have a common interface and offer the same functionality. If we can
> keep up with that goal I don't think we will have issues with
> releasing new versions of the code to both communities at the same
> time. For that reason we chose to put both languages in the same
> repository. We are definitely not opposed splitting them up. If
> someone has a compelling reason for why they should be split
> immediately we will do it, otherwise for our own convenience and to
> enforce parallel development we will probably keep them in a single
> repository for now.
go for it. You can always do the split later, when the package is
mature.
Pablo Winant
I think it's actually a clever solution, which respects the convention of both languages. We can always differ the packaging problem until later, but right now you can install both version using a single command-line that will fetch the code from github. The only strange bit is the package name finishing by .jl . Is it necessary for Julia to consider it as a package repository ?
In my opinion, there are many things to learn from this kind of language comparison. The constraints induced by JIT in Python favours the use of simpler constructs, which produces similar code on both sides. Also, Python/numba is slowly evolving in its own flavour of Python with more advanced features, partly inspired by Julia (see http://continuumio.github.io/numba-nextgen-docs/)
One interesting example is the use of generators. Julia has made it clear that is is sometimes better to "devectorize" code. Python has had for a while a nice syntax for generators. Using such generators, the smolyak products may be constructed without ever constructing a big grid of points explicitly. I think that it would be very elegant, and less memory-intensive, if not faster. Here is the current state of generators:
- pure python: part of the language ( https://wiki.python.org/moin/Generators )
- cython: since 0.15 (https://github.com/cython/cython/wiki/ReleaseNotes-0.15)
- numba: not yet but should be in version 1.0
- julia: unclear to me. There are coroutines with a produce function, but I don't know if it's the right thing to use. There are also macros for fast multidimensional indexing in the "Cartesian" package. It would be worth checking whether they can be adapted. (maybe also for the splines package)
Pablo
Spencer Lyon
I was able to confirm from the Julia guys that the repositories do not need to end in .jl. I have changed the name of the repo to be just Smolyak now. (see my post for discussion).
Also we will have to think carefully and do some testing to come up with the best solution as far as the possible use of generators goes. Like you said Pablo, julia does have co-routines they cal l Tasks, which we are already using a bit in the code where we have generators on the Python side (see util.jl and check out the pmute function and the comb_with_replacement function). I did, however, read on a Julia google group post that the Task, produce, consume interface has a bit of overhead in Julia and we might be better off implementing an iterator by defining the start(), next(), and done() methods for a custom type.
Cartesian is fantastic. We are actually using using one of their macros @forcartesian in the cartesian function in util.jl. I think the @nloops and @nref are quite flexible and will help write concise, de-vectorized, performant code in these splines and interpolation routines.
zac
Pablo Winant
Pablo
Hi,I thought some of you guys might be interested in this. Its linear interpolation Clenshaw-curtis or no boundary grids but as it's based on some stuff I used in matlab the core of itis in a c++ library, this code is also dumped at the above link. It's got a fairly quite routine for generating the grids and aliows you to use different bounds for variables (not [0,1]) and has some rudimentary grid size reduction based on hierarchical surpluses.
--
zac
Spencer Lyon
Hi Zac,
Thanks for the post, this does seem like it could be quite useful! The things Chase and I are working on are still very early in development — as can be seen in the fact that our code doesn’t actually “do” anything besides build the grid and basis functions.
Lately, most of our time working on this project (which is quite limited because we are both first year PhD students) has been spent designing an algorithm that should drastically decrease the times you have reported. We will also have parallel implementations soon (our algorithm makes it trivial to parallelize). Hopefully then these raw benchmarks will be more meaningful.
Do you mind posting the code used to do the benchmark? I am what your object Q is and how that enters into our code.
Also, we are more interested making sure the code is efficient for high dimensionality cases instead of “high density” cases. Therefore, the current code scales much better in what we call d than it does in what we call mu.
Sébastien Villemot
> - is it possible to add a simple build system to the Juila package, in
> order to get rid of the binary ? That could be a simple and very
> useful improvement for other platforms than Linux. Sebastien, do you
> know how to do that ?
BinDeps package. Have a look at other packages (like Cairo) to see how
to build stuff from the deps/ subdirectory.
zac
Chase Coleman
Thanks for putting this up. This looks like a great resource and I think it will be helpful going forward. I'm going to look at your code a little more closely in next day or so, but two quick things:
1. SmolyakGrid(d, mu) does more than just build the grid. It builds the grid and builds the matrix that is used to get the interpolation coefficients; the latter takes significantly longer than building the actual grid. I don't know what is encompassed in yours, but if it is just building the grid then the function build_grid(d, mu) is the comparison that you want to use (I didn't run your code, but it made the times on ours a little closer to what yours is at).
2. I looked through your code somewhat briefly and saw that you have a function called nbgrid_size(q, d). Is that the number of points in the grid hard coded in? If so, then I'm not sure we're building the same grids. Yours seem to have more points than ours. Like I said, I'll look at this a little more closely tonight or tomorrow, but I wanted to let you know, in case you want to look as well.