search

simonbyrne / Remez.jl

Remez algorithm for computing minimax polynomial approximations
Other
38 stars 9 forks source link

SingluarException (not sure why) #6

Open oscardssmith opened 4 years ago

oscardssmith commented 4 years ago

When I run ratfn_minimax(x -> x == 0 ? one(x) : expm1(x)/x, (-1/512,1/512), 2, 0,(x,y) -> x)[3], I get the following error.

ERROR: LinearAlgebra.SingularException(4)
Stacktrace:
 [1] checknonsingular at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/factorization.jl:19 [inlined]
 [2] generic_lufact!(::Array{BigFloat,2}, ::Val{true}; check::Bool) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/lu.jl:178
 [3] #lu!#134 at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/lu.jl:130 [inlined]
 [4] #lu#136 at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/lu.jl:273 [inlined]
 [5] lu at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/lu.jl:272 [inlined] (repeats 2 times)
 [6] \(::Array{BigFloat,2}, ::Array{BigFloat,1}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/generic.jl:1116
 [7] ratfn_equal_deviation(::var"#82#84", ::Array{BigFloat,1}, ::Int64, ::Int64, ::BigFloat, ::var"#83#85") at /home/oscardssmith/.julia/packages/Remez/HLpT8/src/Remez.jl:407
 [8] ratfn_minimax(::var"#82#84", ::Tuple{Float64,Float64}, ::Int64, ::Int64, ::var"#83#85") at /home/oscardssmith/.julia/packages/Remez/HLpT8/src/Remez.jl:610
 [9] top-level scope at REPL[287]:1

That said, ratfn_minimax(x -> x == 0 ? one(x) : expm1(x)/x, (-1/512,1/512), 3, 0,(x,y) -> x)[3] works just fine. Is there a way to make this not happen or work arround it?

simonbyrne commented 4 years ago

Not sure either: unfortunately the code didn't cite any references. From what I understand, this is a known problem with the Remez algorithm. I think there are some potential modifications to work around it, but if you widen the interval so that is no longer symmetric, you get a valid answer:

ratfn_minimax(x -> x == 0 ? one(x) : expm1(x)/x, (-1/512,1/512+1e-8), 2, 0,(x,y) -> x)
oscardssmith commented 4 years ago

That mostly works, but if I only modify the interval slightly, I get very high errors (in the hundreds). Is there any way to prevent this from happening automatically?

simonbyrne commented 4 years ago

Probably, but I don't know.

simonbyrne commented 4 years ago

Might be possible to try something like this: https://people.maths.ox.ac.uk/trefethen/vander_revised.pdf