Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Bromwhich contour: Hyperbola; comments on inversion at t = 0 and fixed contours #23

Open
wants to merge 15 commits into
base: master
Choose a base branch
from

Conversation

heltonmc
Copy link
Member

@heltonmc heltonmc commented May 7, 2021

Hyperbola

This pull requests implements the numerical inversion of the Laplace transform by using a hyperbola for the Bromwhich contour. It's important to say that the inversion is very ill posed so different contours will lead to different convergence and accuracy.

Major Changes:

  • A new file bromwhich_contour.jl was added that utilizes a hyperbola for the Bromwhich contour. hyperbola uses a contour that depends on the time point which yields most accurate results. hyper_fixed computes a single optimal contour for a given range of times needed that is much more efficient at the cost of accuracy.
  • Main file was updated to export hyperbola and hyper_fixed.
  • Test suite was added for known analytical transforms.

Comparison to Talbot contour

As an example, let's compare the convergence of the Talbot and hyperbola contours for t = 1.0 and 10.0 for the functions
F(s) = 1/s - 1/(s + 1), f(t) = 1 - exp(-t).
Screen Shot 2021-05-06 at 8 11 19 PM

Screen Shot 2021-05-06 at 8 07 31 PM

Take aways:

  • Hyperbola contours have a faster convergence determined by the decay slope (on log scale so falls exp(-x))
  • There is an optimal N (number of nodes on the contour evaluated) depending on precision used. For example, the increase after the minimum is due to round-off errors and the ill conditioned inverse transform.
  • The contour and the parameters defining that contour will affect the convergence rate. There can be separate contours and parameters that give the most accurate results depending on the problem.
  • Accuracy of the method is closely related to singularities of transform F(z)
  • Utilizing hyperbola contours instead of Talbot allow for singularities with unbounded imaginary part given close proximity to negative real axis
  • I'm unsure why the hyperbola contour presents a smooth decay that fits exceptionally well with theoretical convergence estimates while the Talbot decay is much noisier.

Inversion at t = 0

I have seen a couple issues regarding this specifically #12 and #5 which has led to some confusion with packages using this JuliaRheology/RHEOS.jl#32 and JuliaRheology/RHEOS.jl#110.

The integrand we are trying to calculate looks like exp(st)*F(s)ds. The exp(st) term is highly oscillatory on the Bromwhich line and F(s) typically exhibits a rather slow decay. So now we have an integral that is both highly oscillatory and decays slowly..... this is hard to compute. The ingenious idea by talbot was to circumvent this slow decay by deforming the contour to begin and end in the left plane, i.e., Re s -> -inf at each end. Now, exp(st) when evaluated at s -> -inf forces a rapid decay of the integrand so we can then approximate the integral by the midpoint or trapezoidal rule.

Now, if t=0 we see that the exponential term no longer forces rapid decaying of the integrand so we must evaluate larger values and more terms.

Alternatively:

From the optimization approach we are trying to determine the contour parameters that give fastest convergence. We choose these parameters to be proportional to N/t. In fact if you look at the source code you will see:

# talbot
r = (2 * bM) / (5 * t)

#hyperbola 
μ = 4.492075287 * N / t

So if t = 0 you get one of the parameters to be infinity which leads to undefined behavior. Thus, the original Talbot implementation was operating appropriately, it just can't be evaluated when t = 0. **Note: talbotarr works when t[1] = 0 because the contour is fixed by r = (2 * bM) / (5 * tmax), so we never end up dividing by zero. This approach though is less accurate though I have not extensively used fixed talbot. The original suggestion was to approximate the value at t=0 by evaluating at a small positive value close to 0. Let's look at the contour we are integrating in the real and positive imaginary plane as a function of time. (only need to consider half of contour for real valued functions).

Screen Shot 2021-05-06 at 9 13 11 PM

If you have really good vision you can see the blue curve when t = 1.0 is a tiny dot near 0. So we can see that the contour size depends on t which is important for the evaluation near 0 as well as the topic I'll discuss later considering a fixed path. For values close to zero we have to accommodate the larger contour with more function evaluations.

Solution

As originally shown by @jlapeyre, we can't directly evaluate at t = 0, but we can approximate at t = 0 by evaluating very close to 0. In #5, we had this original example:

julia> Jbar(s) = (1/s^2)*(1 + s^(0.5))/(1 + 1/s+ s^(-0.5));
julia> hyperbola(s-> Jbar(s), 1e-30) # using hyperbola presented here
1.1283791670955786e-15

I'm not for sure what the analytical f(t) looks like here, but it appears to go to 0 as t -> 0 and towards 1 as t -> inf. Can @moustachio-belvedere comment on how this function should look? But let's quickly look at a case when we know the solution as shown in #12. Again this is expected behavior as discussed previously, and I believe the work around is being implemented correctly and further resolves #5 and resolves #12.

julia> hyperbola(s-> 1/s - 1/(s + 1), 1e-100) # pick t close to zero
0.0 # correct analytical value

Evaluating for many time points

So I wanted to try and discuss different issues originally brought up in #5 but seem to be open in issues JuliaRheology/RHEOS.jl#32 and JuliaRheology/RHEOS.jl#110.

First, talbot and hyperbola require N evaluations of the function F(s) for each time point. The parameters of the contour depend on t, so for each new value of t, F(s) needs to be recomputed on a different contour. This is quite slow when F(s) is computationally expensive or many time points are needed. As I've discussed and shown in the previous figure, the contour size and shape depends on t. In fact, it is optimized for each individual t point for best convergence.

However, this has a drawback because for each new time point we must also compute a different contour and evaluate our function again. What talbotarr and hyper_fixed seek to do is a pick one contour shape that does not depend on t (only the ratio of t[end] and t[1] and evaluate the function only N times instead of N*length(t) times. So now this method is significantly faster proportional to however long it takes to compute F(s) and how many time points needed. For one of my applications it is over 5 orders of magnitude faster.

The question now is how accurate do you need the calculation to be? Instead of having an optimized contour for each time point we now have to compromise to an average best contour which will decrease accuracy. The decrease in accuracy will depend on the ratio between the start and end times. Let us try to show us an example of what I mean:

# set up example for start t = 1.0
julia> t = 1.0:10.0
1.0:1.0:10.0
julia> @test isapprox(hyperbola(s -> 1/s - 1/(s + 1), t), 1 .- exp.(-t))
Test Passed
julia> @test isapprox(hyper_fixed(s -> 1/s - 1/(s + 1), t), 1 .- exp.(-t)) 
Test Passed

# set up example for start value is 0.1
julia> t = 0.1:10.0
0.1:1.0:9.1

julia> @test isapprox(hyperbola(s -> 1/s - 1/(s + 1), t), 1 .- exp.(-t))
Test Passed
julia> @test isapprox(hyper_fixed(s -> 1/s - 1/(s + 1), t), 1 .- exp.(-t)) 
Test Failed

# though we can make this example work by increasing the number of node evaluations
julia> @test isapprox(hyper_fixed(s -> 1/s - 1/(s + 1), t, N = 36), 1 .- exp.(-t))
Test Passed

Ok so we can see that if the range of the start and end t values is too large that increases the size of our contour and will require more evaluations. To resolve JuliaRheology/RHEOS.jl#110, I am not for sure what the true values should be so I can't comment on the best way to invert this. @akabla or @moustachio-belvedere may have to provide expected values. Though using talbotarr as well as hyper_fixed will provide poor results if your start value is very close to zero and your end time value is very far. Going back to the previous example:

julia> t = 0.0001:10.0
julia> @test isapprox(hyper_fixed(s -> 1/s - 1/(s + 1), t, N = 96), 1 .- exp.(-t))
Test Passed

We can increase the number of evaluations so this test passes even with a really low start time. Honestly, this is pretty incredible that this test past when t[end]/t[1] is 100000. This approach will be orders of magnitude faster than using hyperbola or talbot. I'm hoping some of these explanations will help close JuliaRheology/RHEOS.jl#32 as well. One thing to mention that would be a hybrid approach is to utilize 2 contours. Use hyperbola to estimate a value close to zero and then use hyper_fixed for the contour so that the range of time values are smaller. This could actually allow for less computations if hyperbola converges very fast then by decreasing the t[end]/t[1] ratio you would need less evaluations there as well. Though this really depends on how close to zero you need to know f(t). If 0.1 is sufficient probably can just use the same contour but if we need 1e-100.. then might be best to consider two contours.

@jlapeyre There are also seems to be some confusion about some of the implementations of the approaches I've seen through the issues. I would have some time to work on documentation more extensively but wanted to know if you had a vision for this, maybe we can open an issue to discuss separately.

@heltonmc heltonmc marked this pull request as draft May 7, 2021 13:00
@heltonmc heltonmc marked this pull request as ready for review May 7, 2021 19:04
@moustachio-belvedere
Copy link

moustachio-belvedere commented May 16, 2021

Thank you, this seems like a very interesting and valuable PR and I appreciate the highly informative PR message. Thanks also for keeping us in the loop. JuliaRheology/RHEOS.jl#32 is certainly still a topic of interest for us, though work on RHEOS is now somewhat intermittent.

With regards to how the functions should look, there is an example in #12 , the time domain form is J(t) = 1 - exp(-t) which yields 0.0 at t=0.0. The Laplace transform of this function is Jbar(s) = 1/s - 1/(s + 1).

julia> J_i(s) = 1/s - 1/(s + 1);
julia> talbot(s -> J_i(s), 0.0)
NaN

I don't have to hand all the other functions we tested in time domain + Laplace domain but most of them rely on either summed exponentials (similar to example above) or some power law terms. So for another test case this might be a good one (hopefully I have both forms right):

julia> f(t) = t^0.5;            
julia> g(s) = gamma(1.5)/s^1.5;                  
julia> f(0.0)                   
0.0                             
julia> talbot(s -> g(s), 0.0)   
NaN       

Lastly, I would just note that on JuliaRheology/RHEOS.jl#110 I don't believe the title refers to any bad implementation of InverseLaplace.jl itself, rather that us RHEOS developers may have not integrated it into RHEOS.jl in the best way, and the interest in exploring the Weeks algorithm which may be faster in speed-critical part of RHEOS.

I hope the above is helpful, and thanks again!

@moustachio-belvedere
Copy link

Sorry I missed the important part:

I'm not for sure what the analytical f(t) looks like here, but it appears to go to 0 as t -> 0 and towards 1 as t -> inf

That sounds right to me. From physical principles I would expect it to certainly go to 0 as t->0, and if it didn't I would assume I had made a mistake in the analytical form. Unfortunately this was one of the moduli we were not able to analytically invert so I don't have more insight.

@heltonmc
Copy link
Member Author

heltonmc commented Jul 4, 2021

Thanks so much for the comments!! And very sorry for the long delay as I got busy writing proposals. These two examples you provided are excellent as they are so different, and maybe this can help with your discussion on JuliaRheology/RHEOS.jl#110. I'm going to try to discuss this in terms of 1. accuracy and 2. performance.

Accuracy

I attached the salient parts of the code I ran to produce the plots below.

## two function to look at
f1(t) = t^0.5
g(s) = gamma(1.5)/s^1.5

f1(t) = 1 - exp(-t)
g(s) = 1.0 / s - 1.0 / (s + 1.0)

## compute Weeks function weights, performance critical
w = Weeks(s -> g(s), 80)

## define our t array, the ratio of t_end/t_start affects accuracy of fixed hyperbolas 
t = range(0.05, 40.0, length = 60)

## compute absolute difference for hyperbola, talbot, fixed hyperbola, and weeks
## hyperbola, talbot, weeks take a float value, hyperbola fixed requires an array
h_err = hyperbola.(s -> g(s), t) .- f1.(t) # hyperbola 
t_err = talbot.(s -> g(s), t) .- f1.(t)    # talbot
w_err = w.(t) .- f1.(t)                    # weeks
hf_err46 = hyper_fixed(s -> g(s), t, N = 46) .- f1.(t)   # fixed hyperbola using 46 evaluations
hf_err76 = hyper_fixed(s -> g(s), t, N = 76) .- f1.(t)   # fixed hyperbola using 76 evaluations

Screen Shot 2021-07-04 at 10 53 29 AM

Screen Shot 2021-07-04 at 10 55 56 AM

For these two examples, the hyperbola contour converges the quickest and is the most accurate. There is a small caveat that the Weeks method produces surprisingly accurate results for J(t) = 1 - exp(-t) when t < 10.0. The Week's method often times produced relative errors of less than 1e-15 in that window. However, for the other example and for time windows outside of listed it produces by far the worst results. These results would not be acceptable I don't think. On the other hand, both the hyperbola, fixed contour hyperbola, and talbot method produced very accurate results for both examples and for all the time windows. It looks like there are other optional arguments for the Week's method to optimize better, but it appears that the hyperbola and talbot contours don't require as much optimization. I personally would only use Week's method if you have a good way to confirm that it works appropriately for your problem and time window.

Based on these examples and others as well. The hyperbola contour is able to converge quicker and more accurately than the talbot contour. Also interesting is the ability of the fixed hyperbola to match the accuracy of talbot given enough evaluations. In general the larger the ratio of t[end]/t[1] the more evaluations will be needed. I showed here 46 and 76 evaluations. Depending on your application, the 46 evaluations can give about 7 digits of accuracy. Again for a smaller time window, less evaluations would be needed for similar accuracy. It is very strange to note the consistency in the error in the fixed_hyper for both evaluations... I honestly, am not for sure why the error is so stable. It is quite interesting. The advantage of this approach will be highlighted in the performance section.

Performance

t = range(0.05, 40.0, length = 2000)
@btime hyperbola.(s -> g(s), $t)
  9.283 ms (1 allocation: 15.75 KiB)
@btime talbot.(s -> g(s), $t)
  10.472 ms (1 allocation: 15.75 KiB)
@btime w.($t)
  2.212 ms (3 allocations: 15.88 KiB)
@btime hyper_fixed(s -> g(s), $t, N = 46)
  1.866 ms (22 allocations: 18.59 KiB)
@btime hyper_fixed(s -> g(s), $t, N = 76)
  3.152 ms (22 allocations: 19.66 KiB)

In order of increasing computation time: hyper_fixed(N = 46) > weeks > hyper_fixed(N=76) >> hyperbola > talbot. The fixed hyperbola with 42 evaluation is the fastest followed closely by Week's method and then the fixed hyperbola with 76 evaluations. This also doesn't include the time to initialize by the Weeks' method which is computational demanding compared to the evaluation showed here. The performance will be affected by how many time points you need evaluated at and how expensive J(s) is to evaluate. For example, in the hyperbola and talbot method it requires ~16 evaluations for each time point so the computation time will be proportional to 16*length(t). Whereas, for the fixed method you need N evaluations and it is mostly independent by how many t values you need. I showed here 46 and 76 evaluations. So with that said, if you only have a few time points that you need evaluated at hyperbola is a really good method that is the most accurate. If you need to evaluate for many time points, the fixed_hyper method becomes a really good option which can get very close to the accuracy of hyperbola given enough evaluations.

There is this snag when t=0. This is something I don't have too much experience in because my time signals I work with are always 0 at t=0. I personally wouldn't recommend any of these methods for direct evaluation at t=0. You can approximate this as suggested by just inserting a value very close to zero. Again, I would be very careful because your absolute accuracy is limited by the machine precision. So for values that go to zero at t=0 you get below the machine precision which return pretty meaningless numbers. Now, you can use BigFloats but I have found performance to be not sufficient. In my applications, using BigFloats although they generate exceptionally accurate results can take 1000x longer to compute for complicated functions. My recommendation would probably be to have two evaluations... something like hyperbola(s -> g(s), 1e-5) for a very short time and then hyper_fixed(s -> g(s), t, N = 46) where t = range(0.05, 5.0, length = 100) is the rest of your time array. Again, it will depend on the range of you t values and the number of time points you need. For example, I just use hyper_fixed in my applications when I need a time vector between 0.03 and 8.0. I've found this to be accurate enough while being by far the fastest when I use multi-threading.

@akabla
Copy link

akabla commented Jul 28, 2021

Dear Michael @heltonmc ,
many thanks for your thorough analysis of the accuracy and performance. This looks very promising indeed. It also answers some questions I had about the accuracy of the Weeks method in our context. I would be delighted to see these methods integrated into the package.

@jlapeyre
Copy link
Collaborator

It looks like I somehow missed this PR. It's a year old! Thanks for doing this.

I do have a concern about the many floating point literals and conversions in this PR. I'm not certain that they matter. For the Talbot routine they certainly do when you use BigInt and BigFloat for stability.

For example 1/2 is a Float64.

And in the examples, it's better to use g(s) = 1/ s - 1/ (s + 1) than g(s) = 1.0 / s - 1.0 / (s + 1.0). Julia will convert integers to floats when necessary.

@jlapeyre
Copy link
Collaborator

BigFloats although they generate exceptionally accurate results can take 1000x longer to compute for complicated functions.

But, it makes sense that for some applications BigFloat is not possible.

@heltonmc
Copy link
Member Author

No problem! And no worries if this isn't wanted in this package, I should have filed an issue first. It provides similar features to talbot so the maintenance burden might not be what you want (though I'm happy to provide help there).

With that said reviewing this code with a bit more julia experience, I did have a couple concerns I'll mark if you want to go ahead with this PR which I can fix.

"""
function hyperbola(f::Function, t::AbstractFloat; N = 16)
a = zero(Complex{eltype(t)})
N = convert(eltype(t), N)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Probably not needed to convert types here as it will be type stable without this if N is an Int. We should limit N to be an integer as line 57 will not be accurate if it's not.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. For a number of reasons, you don't want to convert N.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. For a number of reasons, you don't want to convert N.

src/bromwhich_contour.jl Show resolved Hide resolved

# parameterized hyperbola contour
# parameters are given by ref [2] which are given in increased precision from calculations given in [1]
function s(θ, N, t)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

s and ds should be combined in a single functions as they share several arithmetic operations. Can save at least one multiplication and division per loop doing this.

I don't think there's an option to return both sinh and cosh as there is a sincos(x) function.


# get the fixed integration components
function hyper_coef(N, t::AbstractArray; ϕ = 1.09)
A = acosh(((π - 2 * ϕ) * t[end] / t[1] + 4 * ϕ - π) / ((4 * ϕ - π) * sin(ϕ)))
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

O no!!! 😲 Considering all the recent blogs and attention on this t[1] seems like a potential problem here considering we allow t to be an AbstractArray. This also assumes t is ordered and t[1] is the start or first value in t.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, but it's really easy to do this. I think first(t) will give you what you want.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But, I'm pretty sure AbstractArray has to be indexable, and the indices fix an order.

dsk = ds((k + 1/2) * h, N, t)
a += f(sk) * exp(sk * t) * dsk
end
return imag(a) * h / pi
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I use pi here and π below which is a style inconsistency.

# compute the function values over the fixed contour nodes
function fixed_sk(f::Function, N, t::AbstractArray)
μ, h = hyper_coef(N, t)
a = zeros(Complex{eltype(h)}, Int(N))
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmmm I always go back and forth on the best way to handle this. But we are allocating quite a bit within this function that the user can't do anything about it. And then we call Threads which greatly increases the allocations as well. If the user is calling this within a loop they'll see a lot of allocations. Wondering if she provide some pre cache option here....

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You should probably used typeof rather than eltype for things that aren't containers. Although eltype works. Maybe because floats and ints are iterable.

It's common to provide fixed_sk! that takes a buffer as input and then fixed_sk just does this allocation and calls fixed_sk!.

I don't know the best way to handle threads; you can provide a keyword argument. You can also have a conditional that turns on threads only if the work is big enough.

@jlapeyre
Copy link
Collaborator

If I understand, these routines provide some significant improvement over other methods for some problems, right? Converge more quickly, or with better accuracy? If so, then I think we should go ahead with this. I am traveling for the next few days. But, ping me if you don't hear from me.

@heltonmc
Copy link
Member Author

@jlapeyre sorry for the delay. I think I've now fixed all these issues and ready to review.

All the functions will now preserve input types so you can use these functions with Float32, Float64, BigFloat just fine and will output correct type. I've also completely re-worked hyper_fixed and also added hyper_fixed! which gives the option of preallocating so the function has no allocations (and is now a much simpler function). This way is definitely more memory efficient than the previous as it doesn't have any internal calculations but I didn't notice any speed improvements in microbenchmarks.

But yes, these are improved versions of the original talbot contours that will converge faster so will provide more accuracy in less evaluations of F(s).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
4 participants