CompEcon function approximation needed for Python

[sebgraves]

Python currently lacks the following CompEcon functions, which are useful for approximating functions using Chebyshev polynomials and splines:

fundef
funnode
gridmake
funfitxy

There may be other related CompEcon functions that are important, but these seem to be the critical ones.

[davidrpugh]

@sebgraves What is the dimension of the functions that these routines can be used to approximate? Many (possibly all) of these routines will exist in numpy.polynomial package or scipy.interpolate (possibly other places as well) I have used NumPy and SciPy to solve functional equations using both splines and Cheyshev polynomials in 2-3 dimensions.

Can you check and see if any of the routines from CompEcon are not covered by either of these packages.

[jstac]

@davidrpugh Good point. However, some potential users are being discouraged by the lack of obvious equivalents. So it will still be valuable to have matching functions that do the same job as found in the CompEcon toolkit. Hopefully those routines you mention can form the back end and do the heavy lifting.

[davidrpugh]

@jstac I see. So wrapping standard numerical routines for functional approximation, interpolation, etc that already exist in other packages is definitely something that we want to do in quantecon? I assume that we are not trying to mimic the CompEcon MatLab API, but rather have obvious equivalents available via a more Pythonic API?

[thomassargent30]

I think that David's suggestion about merely having a clear and prominently available statement of readily available substitutes for key compecon programs could be all that we want.

…

On Wed, Feb 8, 2017 at 1:24 AM, David R. Pugh ***@***.***> wrote: @jstac <https://github.com/jstac> I see. So wrapping standard numerical routines for functional approximation, interpolation, etc that already exist in other packages is definitely something that we want to do in quantecon? I assume that we are not trying to mimic the CompEcon MatLab API, but rather have obvious equivalents available via a more Pythonic API? — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#281 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AD22TnJwCVHRLaAkW04NlEjwqHlJS3BGks5raV-DgaJpZM4L6ORh> .

[jstac]

Thanks for the input. It's hard to know the best way forward. The key point here is that someone was blocked from working in Python due to these missing packages (the back story behind OP's message). It's part of our mission to avoid such outcomes. But this can be addressed in different ways. Options are:

(a) just document alternatives
(b) provide similar, more Pythonic functions, or
(c) provide identical functions

I can see the benefit of (a) and (b), if the documentation and examples are clear. That said, if we have the manpower, I would still lean towards (c). It's easy work for RAs since it's just a straight port, and the whole task is easily quantified. Furthermore, it's easier on the demand side for MATLAB refugees.

But that's just my opinion...

[davidrpugh]

We definitely do not want people to feel blocked from using QuantEcon because their favorite CompEcon MATLAB routine appears to be unavailable.

1. Have RA's to do a first-pass port of MATLAB routines into Python/Julia for reasons that @jstac mentions.
2. Have more experienced devs refactor the output of step 1 into more idiomatic Python/Julia public API and make it clear in the documentation the mapping between Python/Julia API and the CompEcon MATLAB routines. Refactoring would probably have both usability and performance benefits.

I think we want to avoid terse, non-descriptive MATLAB-like public APIs...

[jstac]

Sounds like a good plan to me.

[cc7768]

A few thoughts:

• In case others haven't seen already: @sglyon has already written routines for Julia that are a superset of the things done in the function approximation side of CompEcon (see QuantEcon/BasisMatrices.jl). Much of the original code in this library is based off of CompEcon (with Miranda and Fackler's very generous permission). However, he has put time into writing a nice "Julian" API for interacting with it. There is a companion package CompEcon.jl that emulates the original Matlab API, but I think it is made clear that people should try to avoid that unless necessary (This is in the same flavor of what @davidrpugh was suggesting about avoiding "terse, non-descriptive MATLAB-like public APIs"). This code is prettly well documented and, in my opinion, if we did embark on writing these routines in Python, this would be the right place to start instead of the original Matlab code.
• I agree that we don't people want to feel prevented from using QuantEcon (in any language). A first step, instead of writing a CompEcon "emulation" that just wraps a bunch of scipy functions, might be to put together a tutorial on how to do function approximation in Python. This could take the form of a notebook in the notebook gallery.
• There are a few related projects to getting function approximation routines in Python in addition to scipy.
  • interpolation: This is a package which @albop has developed which has support for splines and complete polynomials. It would be nice not to duplicate his efforts (he has been careful about writing things in numba etc...). I think it would be natural to make contributions to this package rather than start from scratch.
  • Python version of CompEcon: @sglyon pointed this package out to me the other day. It is written by one of Miranda and Fackler's students. I don't know how efficient (or fast) it is, but if someone really would like to use the CompEcon API then this mimics it pretty closely.

Ultimately, some of these suggestions may require a work schedule that is a bit more front loaded than they would experience if they were able to just copy and paste their previous commands, but I think encouraging people to be idiomatically Pythonic (or Julian) will pay large returns in the medium to long run.

[sglyon]

Thanks for those comments everyone.

Here's my take.

Why?

I think there are two main reasons to want a compecon style interpolation library in python (e.g. this is my motivation for creating BasisMatrices.jl that @cc7768 linked to above):

1. To give compecon-Matlab users a familiar starting place when coming to Python. I've heard from many economists that if they could use their favorite compecon routines inside language XX they'd be happy to give XX a shot instead of Matlab. As has been pointed out here, I don't think the answer is to expose the exact same API compecon-Matlab has, but rather (1) make sure the functionality is covered and (2) we have clear instructions (probably in a notebook as @cc7768 suggests) showing how to translate compecon-Matlab based code into whatever solution we come up with.
2. Some algorithms can be built on top of the building blocks provided by a compecon style function approximation library. Without getting too much into the math, compecon-Matlab does function approximation by constructing basis matrices. The final approximate is computed as the kronecker (or row-wise kronecker) product of 1d basis matrices times a coefficient vector. Some algorithms can leverage these basis matrices directly to make algorithms more efficient.

How?

If those reasons are convincing enough that we want to build this library in Python, here's how I would proceed if I were doing it.

I'd start by surveying what's already out there

• There is a branch named albop/linear_basis in the interpolation.py library @cc7768 linked to. This is the start of constructing routines for building up basis matrices for piecewise linear, arbitrary degree B-spline, and Chebyshev interpolation. Performance critical routines are implemented in numba and can be very fast. The developer who takes on this project would be able to work with @albop (and potentially myself) to tidy up these routines and get them ready for mass consumption.
• The compecon-python library by @randall-romero. He's done great work there and has implemented most of the underlying machinery behind making the functions @sebgraves listed "just work". I haven't looked closely at the implementation, but assuming it is efficiently implemented, this might be all we need to do. If performance is lacking in any particular areas, perhaps we could work with @randall-romero to improve it. The developer who takes on this project could do some benchmarking of the @randall-romero routines against compecon-Matlab and BasisMatrices.jl to gauge performance and coordinate with lead quantecon devs and @randall-romero to speed up the slower routines. (a side benefit of this would be pitting BasisMatrices.jl against compecon-Matlab, which hasn't been throughly done).

[albop]

I completely agree with @sglyon and @cc7768 about development plans and with @davidrpugh in the sense that, to me, a clear notebook explaining the various flavours of interpolation and the various available libraries, would be the more useful thing to do.

There are probably several reasons why compecon users find it hard to transition to python: there is the quest for an exact api equivalent (strange since I never found funnode, funfitxy to be especially transparent names) and the fact that interpolation tasks can be described with different vocabulary sets. For instance, for cubic splines alone, one can refer to blending formula, linear bases, or prefiltering/filtering operations...

My initial position was that an efficient interpolation object object was all that was needed, like the one supplied by scipy and that what was happening inside should be an implementation detail. The library interpolation.py was initially intended to provide that for cubic splines and smolyak interpolation in any dimension.

But then, many lecture notes, and some algorithmic papers explain everything in term of basis matrices B and Phi, and people have learnt to use them when coding their algorithms. Also, in some cases, it leads to performance gains to reuse these matrices. This is why I started this summer to write objects to represent bases and products of linear bases and to manipulate them as efficiently as possible. Using objects, it should be possible to do:

basis = CartesianProduct(LineSplines(0,1,10),CubicSplines(0,1,10))   # define 2d interpolator
c = basis.B \ x       # fit coefficients using data x
y = basis.Phi(s)*c  # interpolate them on s

Essentially, this is the work @sglyon has already done in the BasisMatrices.jl library and what the linear_basis branch is supposed to bring. It is evolving slowly as I have tried several things (like code generation or optional typing), but I'm happy to hurry up if there is some interest, or to help anyone who would want to work on it.

[jstac]

Clearly everyone else has thought about this much more carefully than me. It sounds like QuantEcon should be trying to help development of the linear_basis branch and possibly the compecon-python library by @randall-romero. Compared to the naive ports of Matlab code, these tasks will require a more sophisticated coder. Right now the people I know that fall into that set are all very busy. But it's possible we'll hire a strong coder for our predoc position, and they might be able to help push this forward. We'll be making that hire pretty soon.