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.
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
It comes from this line. 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.
The are numbers when you evaluate them. Keep them as variables
1 Like
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.
Yep. I think the idea is to use metaprogramming to remove accessing memory via c[k]
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).
1 Like
stevengj is spot on. This works:
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
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.
1 Like
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)
2 Likes
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:
f = Fun(...)
cdf = cumsum(f)
inv_cdf = Fun(x->roots(cdf-x)[1],0..1)
1 Like
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.
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
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).
Try ApproxFun.sample. Using bisection instead of roots Should fix the issues, look at ApproxFun.bisectioninv
1 Like
Tangentially, do you know of any literature that specifically deals with special functions and its computational considerations? This sounds really interesting!
1 Like
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, 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 … 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.
2 Likes
Numerical Computation of Special Functions, by Gil, Seguro and Temne
http://epubs.siam.org/doi/book/10.1137/1.9780898717822
The DLMF has brief sections on computation as well:
2 Likes
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.
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?
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
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
1 Like
Not necessarily for non-polynomial functions:
2 Likes