order argument

[matthieugomez]

Could you explain a little bit more the order argument in BasisStructure and in funeval? In the documentation, you write "The order field keeps track of what order of derivative or integral the arrays inside vals correspond to". Why is it a matrix? What's the best way to evaluate the zero order, first order, second order derivative of a set of Chebyshev coefficients on a given grid? I'm doing the following.

basis = Basis(Cheb, 50, 0.0, 1.0)
bs0 = BasisStructure(basis, Direct(), nodes(basis)[1], 0) 
bs1 = BasisStructure(basis, Direct(), nodes(basis)[1], 1) 
bs2 = BasisStructure(basis, Direct(), nodes(basis)[1], 2) 
funeval(coef, bs0, CompEcon._check_order(n, 0))
funeval(coef, bs1, CompEcon._check_order(n, 1))
funeval(coef, bs2, CompEcon._check_order(n, 2))

[sglyon]

The order argument is a matrix where each row specifies the order of the derivative to be evaluated along each dimension of the approximant.

So, suppose you have a 2d spline and you want to get the levels and first order mixed partials. The correct way to specify this would be order =[0 0; 1 0; 0 1]. The first row says to compute the levels in each dimension, second says to compute the first derivative along dimension 1 and levels in dimension 2, and the third row says to compute the levels in the first dimension and first derivative in the second dimension.

When there is more than one dimension this must be a matrix so that we can specify multiple orders for each dimension. When we have only 1 dimension we still request that this be a matrix to keep the API consistent.

Notice that in order =[0 0; 1 0; 0 1] there are 6 total basis matrices (length of order), but only 4 unique ones. This library will only compute the 4 unique ones and will simply apply them when needed. So, when it sees the first row it will compute basis matrices for both levels and then store them. Then on the second row it will compute a basis matrix for the first order derivative in dimension 1, but will not recompute the levels for dimension 2.

This is a design inherited from Miranda and Fackler's matlab code. It is a bit clunky to work with, but it is quite flexible and does gain you some efficiency if you have multiple dimensions and want multiple derivative orders.

With all that being said, the way to compute all 3 at the same time would be:

basis = Basis(Cheb, 50, 0.0, 1.0)
x = nodes(basis)[1]
y = sin(x)  # just an arbitrary function

# note the transpose. We want the matrix to have 1 column b/c the basis is 1d
bs = BasisStructure(basis, Direct(), x, [0 1 2]')
coefs = get_coefs(basis, bs, y)

out = funeval(coefs, bs, bs.order)

Now out will be (50, 1, 3) = (length(x), size(coefs, 2), size(bs.order, 1)). The reason for the middle dimension is that you might have fit the coefficients to multiple functions on the same basis (i.e. if y had been a matrix instead of a vector).

• out[:, 1, 1] will have the levels
• out[:, 1, 2] will have the first derivative
• out[:, 1, 3] will have the second derivative

You can check this with maxabs(slice(out, :, 1, :) - [sin(x) cos(x) -sin(x)]), which should be small.

[sglyon]

Oh I should add, you need to pull the latest master down because I fixed a lingering bug here as I wrote the previous message.

[matthieugomez]

Thanks a lot for your answer. I'll try to update the README once I practice it a little bit more. One more question, why does not the following give the higher order derivatives?

basisμ = Basis(Cheb, 50, 0.0, 1.0)
basisσ = Basis(Cheb, 30 ,0.0, 1.0)
basis = Basis(basisμ, basisσ)
S, (gridμ, gridσ) = nodes(basis)
bs = BasisStructure(basis, Direct(), S, [0 0 ; 1 0 ; 0 1 ; 2 0 ; 0 2]) 
#BasisStructure{CompEcon.Direct} of order [0 0]
bs.order
#1x2 Array{Int64,2}:
#0  0

[matthieugomez]

The derivatives are wrong in this example. Is this normal? (I don't have much experience with splines).

x = collect(linspace(0, 1, 10))
basis = Basis(Spline, x, 0, 3)
bs = BasisStructure(basis, Expanded(), x, [0 1 2]') 
funeval(get_coefs(basis, bs, map(exp, x)), bs, bs.order)

[sglyon]

This look like a bug to me...

One small note. When you use higher order spline (more than degree 1), the x you pass in is not actually the x used as the nodes of the spline. To preserve continuity at the points the user passed in, we actually pad the given x with a few extra points. To see what I mean check out nodes(basis)[1].

Let's open a separate issue for the derivatives.

[sglyon]

@matthieugomez I looked into this and we are doing exactly what compecon matlab are.

Here's the matlab:

xx = linspace(0, 1, 10);
basis = fundef({'spli', xx, 0, 3})
x = funnode(basis)
y = exp(x)
c = funfitxy(basis, x, y);
yy = funeval(c,basis, x, [0 1 2]')

And the julia:

xx = collect(linspace(0, 1, 10))
basis = Basis(Spline, x, 0, 3)
x = nodes(basis)[1]
bs = BasisStructure(basis, Expanded(), x, [0 1 2]') 
c = get_coefs(basis, bs, map(exp, x))
yy = funeval(c, bs, bs.order)

Then I compared the two (just copied and pasted output from matlab after doing format long):

julia> maxabs(c - [   1.000000000000000
          1.037037673405857
          1.115221701679381
          1.246281946464564
          1.392743726472456
          1.556417698162149
          1.739326466519631
          1.943730427560073
          2.172156069568134
          2.427424884464883
          2.617606373685348
          2.718281828459046])
1.3322676295501878e-15

julia> maxabs(yy[:, :, 1] - [   1.000000000000000
          1.037731454624735
          1.117519068741864
          1.248848869001682
          1.395612425086090
          1.559623497606781
          1.742908998633458
          1.947734041054676
          2.176629931716248
          2.432425454287207
          2.619446308912387
          2.718281828459046])
8.881784197001252e-16

julia> maxabs(yy[:, :, 2] - [   1.000017181958128
          1.037724394109161
          1.117513292379614
          1.248849111568837
          1.395610882639133
          1.559622330212289
          1.742907282290656
          1.947733213718261
          2.176625056071645
          2.432434719273511
          2.619464916970877
          2.718237278889816])
2.7977620220553945e-14

julia> maxabs(yy[:, :, 3] - [   0.998409595220153
          1.037779860935622
          1.116520392366610
          1.247524353039381
          1.394187526225934
          1.558018530090856
          1.741110607319717
          1.945756158377179
          2.174297003983753
          2.430276933649850
          2.619538404179252
          2.714169139443811])
9.094947017729282e-13

---

So, I think the problem is in the theory, not the implementation here. What compecon does when it takes derivatives of a spline approximated function is differentiate the spline itself, rather than approximating the derivative of the underlying function we want to interpolate over.

[matthieugomez]

Thanks for looking into it. To be clear, this is the output I get

x = collect(linspace(0, 1, 10))
basis = Basis(Spline, x, 0, 3)
bs = BasisStructure(basis, Expanded(), x, [0 1 2]') 
funeval(get_coefs(basis, bs, map(exp, x)), bs, bs.order)
10x1x3 Array{Float64,3}:
[:, :, 1] =
 1.0    
 1.11752
 1.24885
 1.39561
 1.55962
 1.74291
 1.94773
 2.17663
 2.43243
 2.71828

[:, :, 2] =
 -6588.24  
  1766.87  
  -472.53  
   130.758 
   -42.1098
    47.0585
  -135.645 
   507.233 
 -1880.2   
  7028.2   

[:, :, 3] =
      2.0543e5 
 -55038.0      
  14728.8      
  -3869.59     
    757.977    
    847.054    
  -4135.72     
  15707.5      
 -58681.4      
      2.19033e5

i.e. the derivative looks like garbage. It seems like we're getting different results, right?

[cc7768]

Just one thing to mention -- Be careful about using x and nodes(basis)[1]. They aren't necessarily the same. Here you seem to be applying exp to x instead of nodes(basis)[1].

julia> nodes(basis)[1]
12-element Array{Float64,1}:
 0.0     
 0.037037
 0.111111
 0.222222
 0.333333
 0.444444
 0.555556
 0.666667
 0.777778
 0.888889
 0.962963
 1.0   
julia> x
10-element Array{Float64,1}:
 0.0     
 0.111111
 0.222222
 0.333333
 0.444444
 0.555556
 0.666667
 0.777778
 0.888889
 1.0     

[sglyon]

Yep, after taking hint from @cc7768 we get:

julia> x = collect(linspace(0, 1, 10));

julia> basis = Basis(Spline, x, 0, 3);

julia> xx = nodes(basis)[1]; ## NOTE: not the same as x

julia> bs = BasisStructure(basis, Expanded(), xx, [0 1 2]');  # NOTE: using xx

julia> funeval(get_coefs(basis, bs, map(exp, xx)), bs, bs.order)
12x1x3 Array{Float64,3}:
[:, :, 1] =
 1.0
 1.03773
 1.11752
 1.24885
 1.39561
 1.55962
 1.74291
 1.94773
 2.17663
 2.43243
 2.61945
 2.71828

[:, :, 2] =
 1.00002
 1.03772
 1.11751
 1.24885
 1.39561
 1.55962
 1.74291
 1.94773
 2.17663
 2.43243
 2.61946
 2.71824

[:, :, 3] =
 0.99841
 1.03778
 1.11652
 1.24752
 1.39419
 1.55802
 1.74111
 1.94576
 2.1743
 2.43028
 2.61954
 2.71417