OPs for a=0, b=0, c=-1 case and generic lowering methods #18
Conversation
|
For now, here's how we can use this to compute finite dimensional julia> N = 20
20
julia> t = 1.1
1.1
julia> α = zeros(N)'
1×20 LinearAlgebra.Adjoint{Float64,Array{Float64,1}}:
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
julia> α[1] = initialα(t)
0.4430825224938979
julia> αcoefficients!(α,t,2:N)
julia> α
1×20 LinearAlgebra.Adjoint{Float64,Array{Float64,1}}:
0.443083 0.521542 0.555073 0.573826 0.585848 … 0.622648 0.623722 0.624682 0.625545 0.626325
julia> B = BandedMatrices._BandedMatrix(Vcat((-1).^(1:N)' .* α,(-1).^(0:N-1)'), N, 0,1)
20×20 BandedMatrix{Float64} with bandwidths (0, 1) with data vcat(1×20 Array{Float64,2}, 1×20 Array{Int64,2}):
1.0 0.521542 ⋅ ⋅ … ⋅ ⋅ ⋅
⋅ -1.0 -0.555073 ⋅ ⋅ ⋅ ⋅
⋅ ⋅ 1.0 0.573826 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ -1.0 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 0.624682 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ -1.0 -0.625545 ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 1.0 0.626325
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0edit: swapped order of bands to be correct |
Codecov Report
@@ Coverage Diff @@
## master #18 +/- ##
==========================================
+ Coverage 90.15% 92.71% +2.56%
==========================================
Files 3 4 +1
Lines 264 508 +244
==========================================
+ Hits 238 471 +233
- Misses 26 37 +11
Continue to review full report at Codecov.
|
|
This is Uvarov's appproach? http://www.mathnet.ru/links/c8f4a5dd76f9cbfb879b7eaee6440f76/zvmmf7099.pdf |
|
The idea is based on Price https://epubs.siam.org/doi/abs/10.1137/0716073, which technically works for more general cases. But as I'm only doing 1/(t-x) right now, I think it is probably equivalent to what Uvarov is doing? |
|
What happens when t is large? Is there any instability with the second kind Legendre functions by forward recurrence used to define the connection coefficients? |
|
Cool! Busy revamping ....in the process I decided I wanted to be able to |
That sounds great, looking forward to it. |
In principle of course one of those terms in the initial condition goes to 1/0 as I also just pushed an update that simplifies the initial condition, it was needlessly complicated. |
|
No I think the failure is catastrophic (through no fault of your own!), and it's even bad for t not-so-far off the support of the orthogonality measure. julia> # inital value n=0 for α_{n,n-1}(t) coefficients
initialα(t) = t-2/(log1p(t)-log(t-1))
initialα (generic function with 1 method)
julia> # takes a previously computed vector of α_{n,n-1}(t) that has been increased in size and fills in the missing data guided by indices in inds
function αcoefficients!(α,t,inds)
@inbounds for n in inds
α[n] = (t*(2*n-1)-(n-1)/α[n-1])/n
end
end
αcoefficients! (generic function with 1 method)
julia> function αcoefficients(t, n)
α = zeros(eltype(float(t)), n)
α[1] = initialα(t)
αcoefficients!(α, t, 2:n)
α
end
αcoefficients (generic function with 1 method)
julia> using PyPlot
julia> n = 1000
1000
julia> t = 2.0 # Far but not really
2.0
julia> setprecision(256) # Default BigFloat precision
256
julia> semilogy(1:n, abs.(Float64.(αcoefficients(BigFloat(t), n)) - αcoefficients(t, n)), ".r")
1-element Array{PyCall.PyObject,1}:
PyObject <matplotlib.lines.Line2D object at 0x7ff6d0a9fb50>
julia> setprecision(65536) # Increase precision just in case
65536
julia> semilogy(1:n, abs.(Float64.(αcoefficients(BigFloat(t), n)) - αcoefficients(t, n)), ".k")
1-element Array{PyCall.PyObject,1}:
PyObject <matplotlib.lines.Line2D object at 0x7ff6d0aa4880>
julia> xlabel("\$n\$")
PyObject Text(0.5, 47.04444444444444, '$n$')
julia> ylabel("Absolute error")
PyObject Text(51.84335222032334, 0.5, 'Absolute error')
julia> grid()
I think the error is so large (and grows so rapidly) that standard BigFloat precision does not suffice because it eventually converges to the same wrong result as the 64-bit floating-point. P.S. it just so happens that I'm deeply interested in this at the moment. |
|
Fascinating! Consistent convergence to the wrong result is one of those things one reads about but I can't say I've encountered them directly in the wild yet before now. Let me just go through and check all of this on my end as well. |
|
In Gautschi's paper "Repeated modifications of orthogonal polynomials by linear divisors", he states that for large t, a continued fraction can be used for computing the modifications to the recurrence coefficients. I haven't looked into it at all, just suggesting that this might be relevant and could save you time. Cheers. |
|
@marcusdavidwebb Yes, I assume the continued fraction form for the n-th term results by just plugging in the previous term's recurrence, at least that certainly does give a continued fraction here. I'll have to have another look at Gautschi's paper to confirm, though, thanks! @MikaelSlevinsky I think I see where the error is coming from too: The recurrence terms basically all look like |
|
@snowball13 and I just used higher-and-higher |
|
Okay well if that's alright then yes that works as is and it's definitely still fast-ish, even quite high orders in tens of thousands are taking me just a few seconds on first computation with some ridiculous setprecision values. If we want to be cheeky we could even check if the values are getting closer to the wrong result I still want to check a few things to see if I can get it more stable, though. |
|
Since it was also going to be useful for more serious tests, I've now made a direct computation of the coefficients available too. Comparing the two shows that we're buying stability for computing cost julia> # goal parameters
julia> t = BigFloat("2.");
julia> n = 3000;
julia> # recurrence version
julia> setprecision(65536);
julia> αr = zeros(eltype(float(t)), n);
julia> αr[1] = initialα(t);
julia> @time αcoefficients!(αr, t, BigInt.(2:n));
0.782400 seconds (83.98 k allocations: 191.259 MiB, 2.65% gc time)
julia> # explicit version
julia> setprecision(256);
julia> αd = zeros(eltype(float(t)), n);
julia> @time αdirect!(αd,t,BigInt.(1:n));
4.463976 seconds (99.61 M allocations: 5.183 GiB, 14.96% gc time)
julia> errs = abs.(αr-αd);
julia> maximum(errs)
6.099089145333413279831345120482707841585231842798410926507117790845627952530002e-76
julia> scatter(errs,yscale=:log10,legend=false,markersize=2)
And of course we can check that at standard BigFloat precision the recurrence shows the wrong convergence to 2*t we talked about while the explicit version doesn't julia> setprecision(256);
julia> αr = zeros(eltype(float(t)), n);
julia> αr[1] = initialα(t);
julia> αcoefficients!(αr, t, BigInt.(2:n));
julia> αr[end]
3.731428791077769921599369234517361314155892288682987165316973210714883071500137
julia> αd[end]
0.267904542249388512670253160328905937104041057150894248970197013940127741415331And actually depending on how precise we want the coefficients to be, we can save a lot of computing time by using less precision in the direct version while still getting the correct results julia> # lower prec direct
julia> setprecision(65536);
julia> αr = zeros(eltype(float(t)), n);
julia> αr[1] = initialα(t);
julia> αcoefficients!(αr, t, BigInt.(2:n));
julia> setprecision(64);
julia> αd = zeros(eltype(float(t)), n);
julia> @time αdirect!(αd,t,BigInt.(1:n));
1.875301 seconds (56.64 M allocations: 2.102 GiB, 19.43% gc time)
julia> errs = abs.(αr-αd);
julia> maximum(errs)
3.32807740918091828711e-18 |
|
Use the hypergeometric functions as a final condition and recurse backwards? |
|
Don't spoil what I'm working on. 😆 I just need to clean it up a bit, add some tests etc. julia> using SpecialFunctions, HypergeometricFunctions
julia> setprecision(256)
256
julia> t = BigFloat("2.");
julia> n = 10000;
julia> α = zeros(eltype(float(t)), n);
julia> @time backαcoeff!(α,t,BigInt.(2:n))
0.043297 seconds (334.61 k allocations: 14.433 MiB, 47.63% gc time)
julia> abs(α[1]-initialα(t))
6.47712641632083346903976389710029967833700027332721103876777760262644385237182e-78
julia> t = BigFloat("100.");
julia> n = 10000;
julia> α = zeros(eltype(float(t)), n);
julia> @time backαcoeff!(α,t,BigInt.(2:n))
0.010631 seconds (223.57 k allocations: 8.502 MiB)
julia> abs(α[1]-initialα(t))
1.740751338909616831176459088151580383395646093770184472125752404895033027568138e-74 |
TSGut
changed the title
|
Yes, I think it should be fine to merge it as is and add those things after. I've changed the name of the PR, you can merge it at your convenience. |
|
Oh, and did I misunderstand, I thought you already had a way to lower If this isn't the case, I think the approach I used to lower |
|
No if I knew how to do a and b then c would follow.. I don’t have any of them |
|
Okay, I'll do a and b as well then, in a separate PR. |
|
Since it went smoother than expected and is close to done now, I guess I'll just add the As it stands the lowering approach seems to have a speed bottleneck in the computation of entries of the higher Jacobi matrix which from what I can tell is computed via Lanczos. I'll see what I can do about that but the basic functionality is there. This also fixes most of the things you mentioned above, e.g. Camel case typing and some other redundancies. I'll go over it one more time to make sure, though. |
TSGut
changed the title
|
This does all of the lowering things we want now, including iterative lowering, mixed lowering and raising and the special case (0,0,-1) is treated separately with very fast explicit expressions. Comparing with results obtained via raising instead of Lanczos there also doesn't seem to be an issue with stability. So I think this is getting ready to merge @dlfivefifty, maybe you can have another look. |
There seems to be a bug(?) in both 2F1 and 2F1general2 that is causing loss of full Float64 accuracy for certain high inputs unless the inputs are of BigFloat type. Presumably it's happening because of some intermediate allocations in HypergeometricFunctions.jl. Either way, this provides a workaround for that so that this method remains accurate to machine precision in high ranges of a, b and c.
The plan is to generate a=0, b=0, c=-1 case directly from Legendre(). This can later also be expanded to any negative integer c but that will be more complicated.