Should `<` be supported?

[dlfivefifty]

I'm under the impression that a TaylorN should be thought of like a dual number, but there are some inconsistencies:

julia> dual(1,2) < 5
true

julia> x = set_variables("x", order=100)[1]; x < 5
ERROR: MethodError: no method matching isless(::TaylorN{Float64}, ::Int64)
Closest candidates are:
  isless(::Missing, ::Any) at missing.jl:70
  isless(::InfiniteArrays.Infinity, ::Real) at /Users/sheehanolver/.julia/packages/InfiniteArrays/24ELy/src/Infinity.jl:47
  isless(::InfiniteArrays.OrientedInfinity{Bool}, ::Number) at /Users/sheehanolver/.julia/packages/InfiniteArrays/24ELy/src/Infinity.jl:146
  ...
Stacktrace:
 [1] <(::TaylorN{Float64}, ::Int64) at ./operators.jl:260
 [2] top-level scope at none:0

https://github.com/JuliaApproximation/HypergeometricFunctions.jl/issues/11

[lbenet]

Thanks for reporting. You have a point; would it be consistent to define isless caring only about the constant term? So in your example, it would return true

[dlfivefifty]

Yes. Adding that and a few other overrides we can compute Taylor series of Hypergeometric functions: https://github.com/JuliaApproximation/HypergeometricFunctions.jl/issues/11

[lbenet]

Excellent! Please go ahead and make the PR! If by some reason you can't, I'll try to do it, but I'm a bit busy now, so it will need to be delayed to next week.

What are the other few overrides that would be needed?

[dlfivefifty]

They are listed in that other issue:

julia> x = set_variables("x", order=100)[1];

julia> Base.isless(y::Number, x::TaylorN) = isless(y, x.coeffs[1].coeffs[1])

julia> Base.log1p(x::TaylorN) = log(1+x)

julia> Base.eps(::Type{TaylorN{T}}) where T = eps(T)

julia> Base.isless(x::TaylorN, y::Number) = isless(x.coeffs[1].coeffs[1], y)

julia> _₂F₁(1.0,2.0,3.0,(x-1))
 0.6137056388801087 + 0.22741127776023898 x + 0.09111691664014213 x² + 0.03815588885481866 x³ + 0.016444861063214 x⁴ + 0.007233833291493177 x⁵ + 0.003231138806309342 x⁶ + 0.0014605872604416172 x⁷ + 0.000666598108984577 x⁸ + 0.00030663682865957726 x⁹ + 0.00014198800183346362 x¹⁰ + 6.61175831375729e-5 x¹¹ + 3.0937284901315815e-5 x¹² + 1.4537030067464611e-5 x¹³ + 6.856078195624316e-6 x¹⁴ + 3.244138814946803e-6 x¹⁵ + 1.5395500556889938e-6 x¹⁶ + 7.325365567743256e-7 x¹⁷ + 3.4937741952650163e-7 x¹⁸ + 1.669918137650735e-7 x¹⁹ + 7.997399515112675e-8 x²⁰ + 3.8369378241164586e-8 x²¹ + 1.8439091439531812e-8 x²² + 8.874730907546905e-9 x²³ + 4.277408653697482e-9 x²⁴ + 2.064305243595276e-9 x²⁵ + 9.974712684076268e-10 x²⁶ + 4.825345651032184e-10 x²⁷ + 2.336803285122648e-10 x²⁸ + 1.1327664717625831e-10 x²⁹ + 5.495981038395394e-11 x³⁰ + 2.6688027662909334e-11 x³¹ + 1.2970572243188646e-11 x³² + 6.309438715257751e-12 x³³ + 3.071904695721014e-12 x³⁴ + 1.4967552080750751e-12 x³⁵ + 7.296290393609026e-13 x³⁶ + 3.557466214493406e-13 x³⁷ + 1.734781742062317e-13 x³⁸ + 8.4641167744794e-14 x³⁹ + 4.1351988664342396e-14 x⁴⁰ + 2.0245607188192277e-14 x⁴¹ + 9.933435583249127e-15 x⁴² + 4.87745577266054e-15 x⁴³ + 2.3898424070298422e-15 x⁴⁴ + 1.164730557369523e-15 x⁴⁵ + 5.638532846273414e-16 x⁴⁶ + 2.7201988011610697e-16 x⁴⁷ + 1.3204809121747408e-16 x⁴⁸ + 6.541814035827748e-17 x⁴⁹ + 3.3413699485305246e-17 x⁵⁰ + 1.750131222506521e-17 x⁵¹ + 9.158752144614354e-18 x⁵² + 4.592589955603457e-18 x⁵³ + 2.092072012797721e-18 x⁵⁴ + 8.029050207486038e-19 x⁵⁵ + 2.2317192379005413e-19 x⁵⁶ + 2.5459451555342405e-20 x⁵⁷ - 5.194425770307503e-22 x⁵⁸ + 2.508218649398269e-20 x⁵⁹ + 4.726277522731426e-20 x⁶⁰ + 5.099591855430971e-20 x⁶¹ + 4.0237927900139253e-20 x⁶² + 2.4012059810360116e-20 x⁶³ + 9.464022285185436e-21 x⁶⁴ - 8.560789937190647e-23 x⁶⁵ - 4.4530982758298455e-21 x⁶⁶ - 5.0935311392342046e-21 x⁶⁷ - 3.811703326348834e-21 x⁶⁸ - 2.0285756003983875e-21 x⁶⁹ - 5.543686046038847e-22 x⁷⁰ + 3.2889317600861456e-22 x⁷¹ + 6.678686505189806e-22 x⁷² + 6.466364352208879e-22 x⁷³ + 4.578004331394191e-22 x⁷⁴ + 2.410193593061127e-22 x⁷⁵ + 7.056801258420997e-23 x⁷⁶ - 3.118640628674264e-23 x⁷⁷ - 7.253314824525695e-23 x⁷⁸ - 7.383671851819276e-23 x⁷⁹ - 5.553094701417065e-23 x⁸⁰ - 3.2534231869591845e-23 x⁸¹ - 1.3172312889958448e-23 x⁸² - 4.6537002481632535e-25 x⁸³ + 5.852887338436875e-24 x⁸⁴ + 7.526522169523832e-24 x⁸⁵ + 6.5491773946846995e-24 x⁸⁶ + 4.527660619681191e-24 x⁸⁷ + 2.484347522530125e-24 x⁸⁸ + 9.128664351625636e-25 x⁸⁹ - 6.364263035599112e-26 x⁹⁰ - 5.26940020289893e-25 x⁹¹ - 6.36132983252991e-25 x⁹² - 5.492753722287055e-25 x⁹³ - 3.869126645659105e-25 x⁹⁴ - 2.2388952296311785e-25 x⁹⁵ - 9.622509299690597e-26 x⁹⁶ - 1.334906115815213e-26 x⁹⁷ + 3.0045944188318423e-26 x⁹⁸ + 4.508518096605067e-26 x⁹⁹ + 4.3209257579636724e-26 x¹⁰⁰ + 𝒪(‖x‖¹⁰¹)

I think this is pretty incredible. (Thanks @MikaelSlevinsky @jwscook)

[dlfivefifty]

@edwardcao3026 would you feel up to making a PR? I'm also quite busy.

[lbenet]

Indeed, that's pretty impressive! Thanks @MikaelSlevinsky and @jwscook !

[edwardcao3026]

Thanks for the solution! @dlfivefifty

[lbenet]

@edwardcao3026 I might be able to come with a PR; do you still want to submit something (and then I wait), or should I go ahead?

[edwardcao3026]

@lbenet please go ahead. I do not have much to say/write. Thanks.

[lbenet]

@dlfivefifty One question: implementing isless as you propose above, things like x >= 0 returns false, as well as x <= 0 and x == 0. The reason is that <= (using Base methods) evaluates (x < y) | (x == y). The second condition, in particular, is false because all terms of the series are checked.

Do you think this is consistent?

[dlfivefifty]

Hmm, if we followed dual numbers then we would want x == 0 to return true:

julia> x = dual(0.0,1);

julia> x ≤ 0
true

julia> x == 0
true

julia> x ≥ 0
true

Though I agree that this doesn't feel right here.

[dlfivefifty]

Maybe Cassette.jl is the answer here? So a user can choose whether the TaylorSeries act like numbers or like polynomials.

[KeithWM]

I know this is an old issue, but I would say that defining isless between two polynomials is a strange thing to do. If I have f(x) and g(x), what do I expect f(x) < g(x) to return? If anything it should return a function h(x) that equals true for those x for which f(x) < g(x), but that would still be a function.

[dlfivefifty]

But Taylor series are not polynomials. And the Issue makes clear we are discussing the setting where they are used for autodiff.

[KeithWM]

Umm, surely a Taylor series, a function of the form a_0 + a_1 x + a_2 x^2 + ..., is a polynomial. You might intend it for a particular setting, but I'd be really surprised to find 1 - x^2 + O(x^3) > x^2 + O(x^3) evaluate to true.

[dlfivefifty]

Your example the series goes on forever so is not a polynomial. In this package it's a_0 + a_1 x + … + a_n x^n + O(x^{n+1}) which is certainly also not a polynomial.

Again you are ignoring my statement about auto-diff, where this behaviour is essential.

[KeithWM]

Leaving the nomenclature aside then, as I don't really care if we call a sum of infinitely many monomials a polynomial or not.

I appreciate your enthusiasm to have my input on auto-diff, I took a look at where the example you provide actually fails. It appears that HyperGeometricFunctions is doing some case-switching based on the value of the number provided by means of comparison operators, but these fail for the Taylorseries. In DualNumbers there is an isless provided that compares only the 0th order term (the "realpart"). Some things that come to mind:

• In any case, your suggestion should not use the .coeffs fields directly, but Taylor1|N should implement something similar to the value method in DualNumbers
• It's a pity DualNumbers doesn't define an intermediary abstract type between itself and Number. If its methods were defined with regards to such a type, we could inherit from that, and only have to implement a few of the methods in dual.jl`.
• Julia's "astronomical composibility" is awesome, but does occasionaly raise the question who should be adapting to whom.
  • I think defining an unavoidable isless between a Taylor series and a number with behaviour that is not entirely unambiguous is a bad idea. For example, should Taylor1([1, 0, -1], 2) >= 1 evaluate to true?
  • I'm not sure if we can expect HyperGeometricFunctions to call a value function on its arguments before doing the comparison. This would allow us to define the value method for Taylor objects.
• I have my doubts about the use case of passing a high order (i.e. more than just a DualNumber) into a function that has such a case-switching in it. In such cases, you might be better advised to find the correct local approximation function (i.c. _₂F₁maclaurin) and use that directly.

This probably doesn't answer any questions, but I thought I'd share my thoughts on the matter.

[mcabbott]

into a function that has such a case-switching in it

ForwardDiff ran into problems especially with functions which test x==0, allowing this to be true often gave wrong answers. Tests like x>0 don't seem to cause problems in the wild, but may break total order. See https://github.com/JuliaDiff/ForwardDiff.jl/issues/480 and around.

It's a pity DualNumbers doesn't define an intermediary abstract type

ForwardDiff might be the obvious place to consider adding something like this. (DualNumbers isn't much used, I think.) Maybe worth mentioning https://github.com/SimonDanisch/AbstractNumbers.jl too.

[dlfivefifty]

Might point is if you called it ε instead of x you would have no issue with isless … that is when it’s infinitesimally small.

but probably there should be another type HyperDual that could wrap a Taylor1 to avoid this issue

[dlfivefifty]

DualNumbers isn't much used

the plan at one point was to have ForwardDiff.Dual to move to DualNumbers.jl. Just would need someone to take the initiative

[dlfivefifty]

(btw, the use of “Taylor series” in this package seems like a misnomer: that term implies convergence in a neighbourhood of 0. The actual data structure seems to actually represent an “asymptotic series”, ie there is no guarantee the series represents something convergent, rather describes behaviour as x -> 0)

[lbenet]

(btw, the use of “Taylor series” in this package seems like a misnomer: that term implies convergence in a neighbourhood of 0. The actual data structure seems to actually represent an “asymptotic series”, ie there is no guarantee the series represents something convergent, rather describes behaviour as x -> 0)

I am not sure I fully agree with you, though I admit that you have a point. True, the package produces the first $n$ (normalized) Taylor coefficients associated to a function $f(x)$ around a point $x_0$, so the series produced are indeed truncated, as most things are in the (finite memory) computer world are. Yet, this doesn't make the series asymptotic. Convergence involves much more than only the coefficients, since the series may have a finite radius of convergence, converging in a neighbourhood of $x_0$ (which is undetermined), and not beyond that. You are right that in the package we do not consider this point, but still the coefficients produced are those expected, actually by construction.

All in all, I think the package name is ok and descriptive enough, though it avoids important issues. Perhaps a better name would have been TruncatedTaylorPolynomials, but the package is already too old to make such a change... I think it is enough to add some warnings in the documentation for the too naive user 😄.

[KeithWM]

Might point is if you called it ε instead of x you would have no issue with isless … that is when it’s infinitesimally small.

but probably there should be another type HyperDual that could wrap a Taylor1 to avoid this issue

It's nice of you to think I would not have a problem with if it were to have a different name, but I think you underestimate me. Of course 1 + ε^2 <= 1 is still ridiculous.

[KeithWM]

(btw, the use of “Taylor series” in this package seems like a misnomer: that term implies convergence in a neighbourhood of 0. The actual data structure seems to actually represent an “asymptotic series”, ie there is no guarantee the series represents something convergent, rather describes behaviour as x -> 0)

A Taylor series is actually typically used when the point of expansion is not 0. If it were, the series could be called a Maclaurin series.

Whether or not the infinite series would converge, I think it's reasonable to call the first n terms a truncated Taylor series, much like people might truncate a Fourier series.

[dlfivefifty]

1 + ε^2 is greater than 1

[dlfivefifty]

A truncated Taylor series does not have the O(x^n) part. Neither does a truncated Fourier

[KeithWM]

1 + ε^2 is greater than 1

Not according to your definition above though.

But frankly I'm tired of this discussion. I have no power over this package, but do hope my input here is appreciated by those that do. I also hope they don't implement anything irreversible they're going to regret later.

[lbenet]

I also hope they don't implement anything irreversible they're going to regret later.

We try to... but we succeed in making our own mistakes 😄

[MikaelSlevinsky]

@KeithWM, I think you misread something: the proposed missing methods compare, through isless, a TaylorN and a Number. You jumped ahead and discussed comparing two TaylorN structs.

Maybe the side-discussion about x <= 0 is why this stalled?

[lbenet]

Maybe the side-discussion about x <= 0 is why this stalled?

From the discussion above, specially the time span involved, I think it coincided with my computer breaking down, getting another to run, and bad very bad memory.... Sorry.