How to increase precision of solution using Roots.jl and ArbNumerics.jl

[roi.holtzman]

@JeffreySarnoff

Before inverting the function, what is it doing – to what end.

I am inverting it, because I need to solve some differential equations.

@JeffreySarnoff

How does x go from a variable (an argument to f(x,) that constrains bounds) to a function ( to solve for x(f,) )?

I have some value of f, say f_0, and I need to find the x that solves f(x) = f_0. So basically, I am trying to implement the root-finding function to give me this inverse. I am not sure this answers your question…

[roi.holtzman]

@joeldahne

The Newton iterations are converging quadratically but you are only increasing the tolerance linearly, so most of the time it does the same number of iterations. Indeed the Mathematica result is also not accurate to 100 digits.

Thanks! you are indeed correct as I wrote

@roi.holtzman

Indeed the Mathematica result was not accurate to 100 digits (I copied the wrong thing). This is the Mathematica result up to 100 digits:

mat_sol = 5970995261535437979056290157285245402529726426867049646646376914313944654874563187306730863627235781

I still do not understand why I am getting that the deviation from the mathematical result is so big compared to the tolerance I am using

setprecision(ArbReal, bits=2000)

mat_sol = ArbReal(0.5970995261535437979056290157285245402529726426867049646646376914313944654874563187306730863627235781)

f(x) = elliptic_e(x) - 13//10
fp(x) = (elliptic_e(x) - elliptic_k(x)) / (2*x)

f_arb(x) = f(ArbReal(x))
fp_arb(x) = fp(ArbReal(x))
f_df_for_problem = (f_arb,fp_arb)

prob = ZeroProblem(f_df_for_problem, ArbFloat.((2//10, 7//10)))

tol_arr = ArbReal.([1e-10, 1e-20, 1e-40, 1e-80, 1e-160, 1e-320])
sol_diff = zeros(ArbReal{1024}, length(tol_arr))
for (i, tol) in enumerate(tol_arr)
    sol_diff[i] = solve(prob, Roots.LithBoonkkampIJzermanBracket(), atol=tol, rtol=tol, xatol=tol, xrtol=tol)
end

sol_diff .- mat_sol
# 6-element Vector{ArbReal{1024}}:
#  1.2638671375145267156858622267668263146486194631103066439239457343927831910566533896936293116052288778029257199931559187656131795963980324568386040031356265232081227664948940326357833243023160430254889004060107464311764950603434757473331346283872256753547221429720703375142021224453030934691368662734391262496e-13
#  1.7545462322884165336212496421883900334107875403756946307200351056640138408884790322533422106586522815394012458664323675976958448507164834376529961055721542090171459202123696446699269229040068395952584809704740000135811776502905386162369235021465761244329039989233135719311727522067456789461205007606883446998e-17
#  1.7545462322884165336079884568956236214664637691431394465487456318730673086362723578117626680882423863314123512122397648188809999514616284468318022816285430076900650010820250236575985366472186173954302239428583719428163347192246596616159009013210278691731167307196736105564879687806472283979634912581903772354e-17
#  1.7545462322884165336079884568956236214664637691431394465487456318730673086362723578117626680882423863314123512122397648188809999514616284468318022816285430076900650010820250236575985366472186173954302239428583719428163347192246596616159009013210278691731167307196736105564879687806472283979634912581903772354e-17
#  1.7545462322884165336079884568956236214664637691431394465487456318730673086362723578116701247318208247194189009150957932857999968194224787136290211688651052701079349918116855192091866405296436840426522586446300241603087130298569879618160424746435212787624616426553039569708383033944653456941061581000481656896e-17
#  1.7545462322884165336079884568956236214664637691431394465487456318730673086362723578116701247318208247194189009150957932857999968194224787136290211688651052701079349918116855192091866405296436840426522586446300241603087130298569879618160424746435212787624616426552989938655802225368867561742567542124464703805e-17

mat_sol = 5970995261535437979056290157285245402529726426867049646646376914313944654874563187306730863627235781

@joeldahne

Notice how the error jumps from 1e-18 to 1e-52 between one step. Then it doesn’t improve further due to already satisfying the required tolerance.

I fixed my example with the actual Mathematica proper result, but I do not get the same results as you. I do not understand why the deviation in my case is so large compared to yours. Do you?

@roi.holtzman

setprecision(ArbReal, bits=2000)

mat_sol = ArbReal(0.5970995261535437979056290157285245402529726426867049646646376914313944654874563187306730863627235781)

f(x) = elliptic_e(x) - 13//10
fp(x) = (elliptic_e(x) - elliptic_k(x)) / (2*x)

f_arb(x) = f(ArbReal(x))
fp_arb(x) = fp(ArbReal(x))
f_df_for_problem = (f_arb,fp_arb)

prob = ZeroProblem(f_df_for_problem, ArbFloat.((2//10, 7//10)))

tol_arr = ArbReal.([1e-10, 1e-20, 1e-40, 1e-80, 1e-160, 1e-320])
sol_diff = zeros(ArbReal{1024}, length(tol_arr))
for (i, tol) in enumerate(tol_arr)
    sol_diff[i] = solve(prob, Roots.LithBoonkkampIJzermanBracket(), atol=tol, rtol=tol, xatol=tol, xrtol=tol)
end
sol_diff .- mat_sol
# 6-element Vector{ArbReal{1024}}:
#  1.2638671375145267156858622267668263146486194631103066439239457343927831910566533896936293116052288778029257199931559187656131795963980324568386040031356265232081227664948940326357833243023160430254889004060107464311764950603434757473331346283872256753547221429720703375142021224453030934691368662734391262496e-13
#  1.7545462322884165336212496421883900334107875403756946307200351056640138408884790322533422106586522815394012458664323675976958448507164834376529961055721542090171459202123696446699269229040068395952584809704740000135811776502905386162369235021465761244329039989233135719311727522067456789461205007606883446998e-17
#  1.7545462322884165336079884568956236214664637691431394465487456318730673086362723578117626680882423863314123512122397648188809999514616284468318022816285430076900650010820250236575985366472186173954302239428583719428163347192246596616159009013210278691731167307196736105564879687806472283979634912581903772354e-17
#  1.7545462322884165336079884568956236214664637691431394465487456318730673086362723578117626680882423863314123512122397648188809999514616284468318022816285430076900650010820250236575985366472186173954302239428583719428163347192246596616159009013210278691731167307196736105564879687806472283979634912581903772354e-17
#  1.7545462322884165336079884568956236214664637691431394465487456318730673086362723578116701247318208247194189009150957932857999968194224787136290211688651052701079349918116855192091866405296436840426522586446300241603087130298569879618160424746435212787624616426553039569708383033944653456941061581000481656896e-17
#  1.7545462322884165336079884568956236214664637691431394465487456318730673086362723578116701247318208247194189009150957932857999968194224787136290211688651052701079349918116855192091866405296436840426522586446300241603087130298569879618160424746435212787624616426552989938655802225368867561742567542124464703805e-17

I thank you also for this last comment

@joeldahne

Regarding why the version with ArbReal returns ArbFloat I would guess it is because they use the type float(ArbReal) for the computations, which is ArbFloat. For most Julia types float(T) gives you a more capable type numerically so it is very commonly used is numerical code. For ArbReal it instead gives you a less capable type, which is a bit unfortunate.

I do not exactly understand what is more capable vs less capable. I guess it would be nice if the root-finding function would return an ArbReal with the error estimation from the root-finding algorithm.

In my other functions, I try to use ArbReal. Do you think since this is “less” capable, it is a bad choice?
If I were to keep everything ArbReal, would you just convert the value that is returned by solve to ArbReal?

[joeldahne]

I fixed my example with the actual Mathematica proper result, but I do not get the same results as you. I do not understand why the deviation in my case is so large compared to yours. Do you?

You are missing the citation marks around the value you give to ArbReal. If you replace it with

mat_sol = ArbReal("0.5970995261535437979056290157285245402529726426867049646646376914313944654874563187306730863627235781")

it works fine!

I do not exactly understand what is more capable vs less capable. I guess it would be nice if the root-finding function would return an ArbReal with the error estimation from the root-finding algorithm.

With more capable I mean that ArbReal is more capable than ArbFloat, it can do anything ArbFloat can do and at the same time it keeps track of the error. To say a bit more about this: The promotion system in Julia typically returns the more “general” of the input types, for example promote_type(Int32, Int64) = Int64 and promote_type(Float64, BigFloat) = BigFloat. In base Julia I believe it always holds that promote_type(T, float(T)) = float(T). But for ArbReal this is not the case, then float(ArbReal) = ArbFloat but promote_type(ArbReal, ArbFloat) = ArbReal.

In your case I believe you can just as well use ArbFloat. The only benefit with ArbReal is that you would notice if you lose too much precision in the numerical calculations, but it wont tell you anything about the error bound for the root anyway. If you want to get error bounds for the root you have to use some other algorithm. For example the interval Newton method, which is what the refine_root function I used in my previous post makes use of.

[dlakelan]

@joeldahne

You are missing the citation marks around the value you give to ArbReal.

This is the same error I mentioned early on. To understand what’s going on here, the Julia Parser reads the number as a Float64 thereby ignoring all the excess digits, then passes that to ArbReal … If you put the quote marks, the parser reads it as a string and passes the string to ArbReal which converts it to a full precision value. You simply can’t give a naked string of digits to the parser and get anything other than a native float type. The parser doesn’t know anything about Arb types.

[roi.holtzman]

Thanks so much for the answer! It is very clear.

@joeldahne

If you want to get error bounds for the root you have to use some other algorithm. For example the interval Newton method, which is what the refine_root function I used in my previous post makes use of.

It would be great to get the error. This refine_root method is part of your package ArbExtras right? Do you know if there is a similar analog in Roots?

And thank you both @joeldahne and @dlakelan for pointing out the issue with the string parsing!

[j_verzani]

Thanks for this informative response. Roots uses float for the input values as it seems like the right generic for different number types.

When I try to solve this with A42 which might be recommended as there is a known bracket, the convergence is quite quick for smooth enough functions, and the tolerance is like 2abs(r)*eps(one(r)) between iterates, I get an answer that is off from mat_sol, but yields an exact zero for f:

julia> setprecision(ArbFloat, bits=1000)
1000

julia> f(x) = elliptic_e(x) - 13 // 10
f (generic function with 1 method)

julia> r = find_zero(f, ArbReal.( (0.0, 1.0) ), A42())
0.5970995261535437979056290157285245402529726426867049646646376914313944654874563187306730863627235781167012473182082471941890091509579328579999681942247871362902116886510527010793499181168551920918664052964368404265225864463002416030871302985698796181604247464352127876246164265529899386558022253688676

julia> f(r)
0

julia> r - mat_sol
-2.496873060882422448783512389107976765165905935558302891285924378533353483003135167996405765830182705071755548296901317683611453726541415832228563754139160817061255022802266605200253564787212375383573447010061344197774631132438257432457875535296194985152788292500142950589819644532161239691396597893509e-78

I didn’t check with other zero-finding methods, and am not sure this happens with ArbReals which should be guaranteed intervals (as I understand it), but over floating point numbers some functions can have 0, 1 or more “zeros” near a mathematical zero.

[roi.holtzman]

@j_verzani

Roots uses float for the input values as it seems like the right generic for different number types.

Yes, that seems reasonable. Is there a way to get back from Roots the error estimate of the found solution?

[joeldahne]

@roi.holtzman

It would be great to get the error. This refine_root method is part of your package ArbExtras right? Do you know if there is a similar analog in Roots?

The refine_root is indeed part of my ArbExtras.jl package, which uses Arblib.jl instead of ArbNumerics.jl for the interface to Arb. The code is fairly short and self-contained, so you could take a look at it.

In general you need to use different algorithms to get rigorous error bounds and I think that would be outside the scope of Roots.jl. I’m only aware of one package in the general repository that does this, IntervalRootFinding.jl. But that package makes use of IntervalArithmetic.jl for the arithmetic and that doesn’t support the elliptic functions that you use. So if you want rigorous error bounds you will either have to use ArbExtras.jl, which is maybe not optimal since it is not in the general repository, or implement your own version of the interval Newton method, which is in principle a very simple method but might be hard if you are not used to working with rigorous numerics.

[roi.holtzman]

Thanks so much for the information! I will have a look at refine_root.