# Evaluate integral on many points (Cubature.jl ?)

**URL:** <https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723>\
**Category:** Numerics\
**Created:** [January 27, 2017, 6:48pm UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723 "2017-01-27T18:48:42Z")\
**Posts on this page:** 20\
**Page:** 2

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [January 30, 2017, 8:17pm UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723/21 "2017-01-30T20:17:33Z")

</div>

> [@dpsanders](#):
>
> Then you could just move plan\_evaluate inside the constructor as well?

I suppose there could be an `evaluateplan` field, though I’m hesitant to add fields to `Fun` at the moment. This would also be a bad idea if the plan involved compilation, as then every manipulation would be extremely slow.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [January 30, 2017, 8:25pm UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723/22 "2017-01-30T20:25:19Z")

</div>

I don’t see why it should by exponentially large. Isn’t the goal to manually unroll the for loop, evaluating the coefficients explicitly? This should just result in `n` lines

---

<div class="post-metadata">

**Author:** ![piever](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/piever/32/1815_2.png) [@piever](https://discourse.julialang.org/u/piever)\
**Post date:** [January 31, 2017, 1:23am UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723/23 "2017-01-31T01:23:03Z")

</div>

It comes from [this line](https://github.com/JuliaApproximation/ApproxFun.jl/blob/development/src/LinearAlgebra/clenshaw.jl#L30). The length of `bk1` and `bk2` grows like the Fibonacci sequence, i.e. exponentially. You can try for yourself and compare `ApproxFun.@clenshaw` and `Math.@horner` with `macroexpand`, the difference is quite striking. On the other hand I think that there should be some recursive formulas to convert from Chebyshev series to power series in a stable way, I’ll give it a try tomorrow and post the results.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [January 31, 2017, 1:54am UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723/24 "2017-01-31T01:54:23Z")

</div>

The are numbers when you evaluate them. Keep them as variables

---

<div class="post-metadata">

**Author:** ![piever](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/piever/32/1815_2.png) [@piever](https://discourse.julialang.org/u/piever)\
**Post date:** [January 31, 2017, 2:25am UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723/25 "2017-01-31T02:25:25Z")

</div>

I see what you mean about the conversion to power series being unstable. What I meant about the macro is that, as it is implemented, the macro leads to expressions of exponentially increasing length. Of course with the normal function implementation this doesn’t happen because those expressions get evaluated at every step. It could still be that a more careful use of metaprogramming could lead to some improvement in performance, I’ll try tomorrow more extensively and I’ll post the results if I get to something interesting.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [January 31, 2017, 2:29am UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723/26 "2017-01-31T02:29:19Z")

</div>

Yep. I think the idea is to use metaprogramming to remove accessing memory via `c[k]`

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 31, 2017, 3:05am UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723/27 "2017-01-31T03:05:29Z")

</div>

> [@dlfivefifty](#):
>
> Yep. I think the idea is to use metaprogramming to remove accessing memory via

You also remove the overhead of a loop by unrolling it.

It sounds like the exponential length increase comes because you are including the same expression more than once, which also hurts performance (unless the compiler can eliminate the common subexpressions). One thing you could do is to assign each expression to a temporary variable before using it. Alternatively, just overwrite the variables exactly as you would in a loop (but without the runtime loop, just an unrolled sequence of statements).

---

<div class="post-metadata">

**Author:** ![piever](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/piever/32/1815_2.png) [@piever](https://discourse.julialang.org/u/piever)\
**Post date:** [January 31, 2017, 3:38am UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723/28 "2017-01-31T03:38:51Z")

</div>

stevengj is spot on. This works:

```julia
clenshaw(x,c) = clenshaw_halved(2*x, c)
@generated function clenshaw_halved{N, R<:Real}(x::R, c::StaticArrays.SVector{N,R})
    a,b = :(zero(R)),:(zero(R))
    as = []
    for k = N:-1:2
        ak = Symbol("a",k)
        push!(as, :($ak = $a))
        a = :(muladd(x,$a,c[$k]-$b))
        b = :($ak)
    end
    Expr(:block,
    as...,
    :(muladd(x/2,$a,c[1]-$b)))
end

```

and seems to produce a substantial speed up (~ 30%). You need to give a static array as vector of coefficients so that the generated function can dispatch on the length of the array.

Btw, if you think it would be useful addition I can make a PR to ApproxFun.jl

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [February 16, 2017, 9:42am UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723/29 "2017-02-16T09:42:46Z")

</div>

A related question – does ApproxFun provide a convenient inversion function, so that for monotonic function F(x) = p, we get a function F^{-1} such that F^{-1}(p) = x?

The idea would be to use ApproxFun to approximate a probability density, cumsum to transform it to the CDF, and then find a quantile function for convenient calculation of intervals, the median, etc.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [February 17, 2017, 1:20am UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723/30 "2017-02-17T01:20:47Z")

</div>

you can treat any function as a PDF and call `sample(f,100)` to do precisely what you want (using bisection to invert the CDF)

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [February 17, 2017, 1:24am UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723/31 "2017-02-17T01:24:05Z")

</div>

For things besides sampling: there isn’t an inverse function yet, though if you know it’s smooth you can do it manually using roots:

```julia
f = Fun(...)
cdf = cumsum(f)
inv_cdf = Fun(x->roots(cdf-x)[1],0..1)

```

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [February 17, 2017, 7:45am UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723/32 "2017-02-17T07:45:39Z")

</div>

Okay, cool, thanks! Simply inverting manually was my first thought as well, but I was lost on implementation until I saw your response. I also ran into the road block that even though the test function below is continuous, roots found many false results.  
My test code used a folded normal for simplicity.

```nohighlight
evaled = Array{Float64}(0)

function pfnorm!{T <: Real}(x::T, evaled::Array{Float64,1})
  push!(evaled, x)
  exp(-x^2)
end

f = Fun(x -> pfnorm!(x, evaled), 0..Inf)
F = cumsum(f) / sum(f)
Q = Fun(x -> roots(F-x)[1], 0..1)

length(evaled), maximum(evaled)
# (248, 2.66e+4)
# Computing my density is costly, and increasingly so away from the MLE.
# Wanted to get an idea of the necessary range and # evaluations.
# Restricting the range to non-negligble prob mass seems like a good idea.
roots_test = roots(F - 0.3)
length(roots_test)
# 100
sum(roots_test .== 0.0) 
# 10
sum(roots_test .== Inf) 
# 10
sum(isapprox.(F(roots(F - 0.3)), 0.3))
# 1
root = roots(F - 0.3) |> y -> y[isapprox.(F(y), 0.3)]
F(root)
# 0.3

```

Simply running root reveals 10 0.0s, and 10 Infs – they’re the first and last 10 items, respectively. Thankfully the actual root is buried in the vector, so its just a matter of indexing it.

Trying to create an ApproxFun function using the anonymous rootfinder function above is slow (waited several minutes before killing it), although that function itself is relatively quick to evaluate. This suggests simply sticking with something along the lines of

```julia
function quantile{Tx <: Real, Tf <: ApproxFun.Fun}(F::Tf, x::Tx)
  for i ∈ roots(F-x)
    if isapprox(F(i), x)
      return(i)
    end
  end
  error("No root found.")
end

quantile.(F, [.025, .25, .5, .75, .975])
# 0.0222, 0.225, 0.477, 0.813, 1.58

```

Any more details on sample?  
`methods(sample)`  
Suggests it takes array inputs, rather than functions. I could make do by define a grid of points, and using f() to assign weights.

Or of course implement a metropolis algorithm. Then that invites comparison to Stan’s implementation of variational Bayesian inference as an alternative means of approximate sampling (my density function is an asymptotic approximation).

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [February 17, 2017, 7:54am UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723/33 "2017-02-17T07:54:51Z")

</div>

Try `ApproxFun.sample`. Using bisection instead of `roots` Should fix the issues, look at `ApproxFun.bisectioninv`

---

<div class="post-metadata">

**Author:** ![miguelraz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/miguelraz/32/631_2.png) [@miguelraz](https://discourse.julialang.org/u/miguelraz)\
**Post date:** [February 17, 2017, 1:24pm UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723/34 "2017-02-17T13:24:53Z")

</div>

Tangentially, do you know of any literature that specifically deals with special functions and its computational considerations? This sounds really interesting!

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 17, 2017, 6:57pm UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723/35 "2017-02-17T18:57:18Z")

</div>

> [@miguelraz](#):
>
> Tangentially, do you know of any literature that specifically deals with special functions and its computational considerations? This sounds really interesting!

Yes and no. I don’t know offhand of any good review or textbook that deals with computing special functions in _general_, though such a thing may well exist. (e.g. I haven’t read [this book by Zhang](https://www.amazon.com/Computation-Special-Functions-Shanjie-Zhang/dp/0471119636), but it might be helpful.) For every _individual_ well-known special function, however, there are lots and lots of papers on computing it efficiently.

For example, one function that I was looking at recently was the exponential-integral function, and I collected some notes on it as well as a Julia notebook implementing it efficiently (the preliminary version is 5x faster than the Fortran code in SciPy) here: [exponential integral (Ei, E₁, Eₙ...) function · Issue #19 · JuliaMath/SpecialFunctions.jl · GitHub](https://github.com/JuliaMath/SpecialFunctions.jl/issues/19) … however, I was mainly concerned about implementing it for complex arguments. For real arguments, there are additional techniques as I mentioned above, e.g. if you google “minimax rational approximation” you’ll find lots of links.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [February 18, 2017, 12:05am UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723/36 "2017-02-18T00:05:28Z")

</div>

Numerical Computation of Special Functions, by Gil, Seguro and Temne

[http://epubs.siam.org/doi/book/10.1137/1.9780898717822](http://epubs.siam.org/doi/book/10.1137/1.9780898717822)

The DLMF has brief sections on computation as well:

[http://dlmf.nist.gov](http://dlmf.nist.gov)

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [February 19, 2017, 2:56pm UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723/37 "2017-02-19T14:56:51Z")

</div>

Great – both those are perfect.

Now, before I start digging into what exactly is going on under the hood when I call Fun, perhaps someone can help me. My major concern is speed. A few details:

- So far Fun has logged over 4,700 cpu minutes. This is more than a hundred fold worse than my initial conservative estimates; function calls should average at under 10 cpu seconds.
- The purpose of what I’m doing is speed: trading accuracy for speed by using an asymptotic approximation instead of standard MCMC. I’ve already lost precision upstream, so maintaining machine precious downstream is preserving accuracy I don’t have in the first place.
- The number of evaluations is some sort of function of complexity (non-smoothness) of the target function. Clearly, number of evaluations can be scaled.

I’m also wondering if the specific behavior of my function may somehow cause issues. Perhaps the use of Optim can lead to slightly jagged (non-smooth) behavior when evaluating close points in succession, when viewed at a fine enough precision? That is, Optim lands somewhere in the neighborhood of the minimum before meeting the “|g(x)| \< 1.0e-08” criteria and halting the search. Where in the neighborhood it lands is erratic; noise that can not be modeled → problems for a function intending to be accurate to machine precision.  
[Also, the function stores the minimizer so that it is used as the initial value in future calls. After a couple indirect tests, I don’t think storing is likely to be the issue.]

Any suggestions?

EDIT:  
Basic example of the problem.

```julia
evaled = Array{Float64}(0)
function pfnorm!{T <: Real}(x::T, evaled::Array{Float64,1})
  push!(evaled, x)
  exp(-x^2)+1e-12randn()
end
f = Fun(x -> pfnorm!(x, evaled), 0..10)
evaled
# 32767 element array
4sum(f)^2
# 3.14159265

```

Adding noise on the order of 10^-12 increases the number of evaluations from 119 to 32,767 – more than 200 times the evaluations, because of irrelevant noise!  
Rounding off noise makes matters worse, increasing the number of evaluations to over 4 million.

Perhaps I should give up on ApproxFun and simply evaluate the function a hundred or so times along an interval, and then fit a high order polynomial via regression?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [February 19, 2017, 5:44pm UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723/38 "2017-02-19T17:44:30Z")

</div>

> [@Elrod](#):
>
> Perhaps I should give up on ApproxFun and simply evaluate the function a hundred or so times along an interval, and then fit a high order polynomial via regression?

Gaussian quadrature will almost certainly dominate that (in the sense of smaller error for a given number of points, or fewer points for the same error). Read up on it, then use the excellent

> **[GitHub - JuliaApproximation/FastGaussQuadrature.jl: Julia package for...](https://github.com/JuliaApproximation/FastGaussQuadrature.jl)**
>
> Julia package for Gaussian quadrature. Contribute to JuliaApproximation/FastGaussQuadrature.jl development by creating an account on GitHub.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [February 19, 2017, 5:53pm UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723/39 "2017-02-19T17:53:27Z")

</div>

Try plotting the coefficients, you may have a long tail of noise.

The easiest solution is fo explicitely specify the number of quadrature points: `Fun(f,0..1,n)`

The slightly harder solution is to expose support for an explicit tolerance, and make a pull request

The even better solution would be to copy the Chebfun heuristic convergence test

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [February 19, 2017, 5:57pm UTC](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723/40 "2017-02-19T17:57:24Z")

</div>

> [@Tamas\_Papp](#):
>
> Gaussian quadrature will almost certainly dominate that

Not necessarily for non-polynomial functions:

[http://eprints.maths.ox.ac.uk/1116/1/NA-06-07.pdf](http://eprints.maths.ox.ac.uk/1116/1/NA-06-07.pdf)

[Previous page](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723.md?page=1)

[Next page](https://discourse.julialang.org/t/evaluate-integral-on-many-points-cubature-jl/1723.md?page=3)