ENH: better build B function

[sglyon]

[cc7768]

This might be some food for thought. The recursive building of the chebyshev polynomials probably isn't our limiting factor, but the nth chebyshev polynomial of the first kind evaluated at a point x is equal to cos(n*arccos(x)).

Like I said, the recursion might not be a limiting factor but cos function is a little faster (and it might be easier to implement with the map idea from earlier today):

In [24]: %timeit np.cos(6*np.arccos(.75))
100000 loops, best of 3: 3.97 µs per loop

In [25]: %timeit chebyshev.chebval(.75, [0, 0, 0, 0, 0, 0, 1])
100000 loops, best of 3: 19.2 µs per loop

In [26]: %timeit math.cos(6*math.acos(.75))
1000000 loops, best of 3: 278 ns per loop

[albop]

Are you saying it's faster to repeatedly compute chebychev polynomials with the direct formulas instead of using the recursive ones ? That would be very strange.

[cc7768]

It depends on how we are setting up the evaluation at each point. This is a part that I understand, but probably not as thoroughly as you do since this is my first time putting some code together for this method.

One of the things we are thinking about right now is the fact that only the chebyshev polynomials of degree > 1 need to be evaluated. We are trying to create B by using the same set of indices as we used to evaluate the grid, but we add one more step and create a set of indices (poly_inds) that tell us which degree chebyshev polynomials we are evaluating.

It was just a thought and a couple of ideas are still bouncing around. Going to try to implement one tonight if I finish a satisfactory amount of my current homework.

[sglyon]

One thought I had that isn't totally fleshed out yet, but that I wanted to document is to define a class for the smolyak polynomial. The class would be very basic and would define __getitem__ in python and getindex in julia to return slices of the interpolation matrix that are computed on the fly.

The idea still needs some work, but there is some potential there.

[cc7768]

Several things:

• Most of our indexes are 1. I finally think I have a way to take advantage of this. Once we have d and mu then we know how many Chebyshev polynomials we will need, I will refer to this as num_cheb. The first d*(num_cheb-1) columns are simply the grid evaluated at each of the Chebyshev polynomials from 2 to num_cheb (To make indexing easier, I think putting the column of ones at the end will be easiest). So you were right Pablo, I think that using the recursive formulation is going to be optimal. This step will take care of all of the polynomial indexes where all but one value is a 1.
• The other terms that we have to take care of are the cross terms. This is a little more difficult to explain without drawing matrices, but I think basically what we will want to do is multiply the num_points x d subsections that correspond to each of the specific Chebyshev polynomials by the next polynomials corresponding subsection offset by a step. For example. If we wanted to get the block of B that corresponded to all of the points being evaluated at the 2nd polynomial and the third polynomial then we could do: cheb_block_2[:, :] * np.roll(cheb_block_3[:, :], 1). I'm pretty sure this works for at least the case of 2 term cross products. Furthermore, I think that doing a similar thing for the higher term cross products will work, just product each of the matrices rolled back by 1, 2, ..., n-1. Sorry if this is convoluted. I will try and get an example up soon.
  ^ Not sure this one works, but something similar to it will.