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 indual.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.

[lbenet]

But, to sum of what's needed, the idea is to define isless from the zeroth order coefficient, right?

[dlfivefifty]

I think if we decide x is infinitesimal we would have the following behaviour:

1. 1 + x + O(x^2) < 2 returns true
2. 1 + x + O(x^2) < 1 returns false or throws an error (since x may be positive or negative and this is inconclusive)
3. 1 - x^2 + O(x^3) < 1 returns true

(EDIT: I changed isless to < since isless requires a total ordering)

[KeithWM]

I know I said I was tired of the discussion, but given the traction it's getting, I'm going to chip in anyway.

@MikaelSlevinsky I did realize we were talking about comparing TaylorN object to numbers.

@dlfivefifty I think the three examples you post make sense. They seem to fit the notion that given an infinitessimal x, you observe the condition within (-x, +x) and conclude true i.f.f. the condition holds almost everywhere (meaning the only case where the condition doesn't hold, at x=0, has measure zero). But this implies that 1 + x + O(x^2) is neither greater than nor smaller than 1. Total ordering is not really my branch of mathematics, but I suppose/fear this implies we must say 1 + x + O(x^2) == 1 holds true.

Maybe in conclusion: first the zeroth order term is considered, if that is strictly greater or smaller, that is the conclusion. If they are equal:

• If the first nonzero term after that is odd, we conclude based on the zeroth order alone that the two are equal.
• If the first term corresponding to a nonzero power is even, the sign of that term is included:
  • isless(1 - x^2 + O(x^3), 1) returns true
  • isless(1 + x^2 + O(x^3), 1) returns false
  • isless(1, 1 - x^2 + O(x^3)) returns false
  • isless(1, 1 + x^2 + O(x^3)) returns true

Higher orders than the zeroth and first nonzero term are not considered.

Of course to TaylorN objects it becomes a more complicated matter than these 1D examples.

[lbenet]

I agree, it's better to think and implement < rather than in isless. I also agree that x, the independent variable of whatever Taylor-polynomial type, should (must?) be considered as infinitesimal (small?).

Now, in order to be concrete for the implementation (thanks @KeithWM for the algorithmic summary), I'd like to point out that, as it stands, it applies directly to Taylor1 objects, but it's a bit subtle for TaylorN. Uncomfortable/subtle examples (for me) are: 1+x-y < 1, or 1+x^2-y^2 < 1, 1+x^2-2y^2 < 1. In terms of what @KeithWM wrote, all those cases should return false because there are sets (of positive measure) which do not fulfill the condition. So, to get a consistent answer in those cases, all coefficients of the first non-zeroth order should have the same sign. Would that be consistent?

As a side remark, another way to get around this would be to use intervals in the evaluation of the polynomial part. We would then need to define some sort of default "infinitesimal domain" for the infinitesimal variables. In some context, the [-1,1] interval is often used, which doesn't sound too infinitesimal...

[dlfivefifty]

Is there such a thing as "infinitesimally-small interval arithmetic"?

[lbenet]

Not that I know...

[lbenet]

I just opened #315.

So far, only the comparisons < and > involving a Taylor1 and a number are there and seemingly work. Essentially, I use the constant term of the series (a0), its eps value, and the sign of the leading order correction (a[k]*eps^k), and consider (for < the max of a0+eps(a0) and a0-eps(a0). I think this works fine.

julia> s, c = sincos(Taylor1(5))
( 1.0 t - 0.16666666666666666 t³ + 0.008333333333333333 t⁵ + 𝒪(t⁶),  1.0 - 0.5 t² + 0.041666666666666664 t⁴ + 𝒪(t⁶))

julia> c < 1
true

julia> c > 1
false

julia> c == 1 # compares the whole series with 1
false

julia> s < 1
true

julia> s < 0 
false

julia> s > 0 
false

julia> s == 0 # ordering is partial !
false

Any ideas how to generalize this to N variables?

Is there such a thing as "infinitesimally-small interval arithmetic"?

I though about using the symmetric box eps(a0)*[-1,1]^n as a "small interval", and evaluate the constant term and the next leading order there. But, for cos(x-y) the comparison < 1 returns false...

Other ideas?

[KeithWM]

There should really be some tests for these (in)equalities. The tests not only ensure correct behaviour, but also make for a clear requirements definition.

Also, cases where the first non-zero term is a higher order than 1 or 2 should also be considered.

[lbenet]

There should really be some tests for these (in)equalities. The tests not only ensure correct behaviour, but also make for a clear requirements definition.

Sure, tests will come later... I'm dealing now with the "proof of concept" part now. 😄

Also, cases where the first non-zero term is a higher order than 1 or 2 should also be considered.

For the Taylor1 part, this is already included:

julia> t = Taylor1(8)
 1.0 t + 𝒪(t⁹)

julia> p < 1
false

julia> p > 1
false

julia> p < 2
true

julia> q = 1 - 0.5*t^4
 1.0 - 0.5 t⁴ + 𝒪(t⁹)

julia> q < 1
true

[lbenet]

I just pushed a commit to #315, allowing to use < and > with TaylorN, with the idea I outlined above; again, higher order terms are considered...

julia> using TaylorSeries

julia> x, y = set_variables("x", numvars=2, order=6);

julia> sn, cn = sincos(x+y);

julia> sn < 1
ERROR: MethodError: no method matching isless(::TaylorN{Float64}, ::Int64)

#=
This is because `IntervalArithmetic.jl` is not a direct dependency of `TaylorSeries`, but if it is 
loaded, the specific methods will be
=#

julia> using IntervalArithmetic

julia> sn < 1 # ok
true

julia> sn > 1 # ok
false

julia> cn > 1 # so far, so good!
false

julia> cn < 1  # grrrr.... it should be true
false

julia> p = 1 - 0.5*(x-y)^3
 1.0 - 0.5 x₁³ + 1.5 x₁² x₂ - 1.5 x₁ x₂² + 0.5 x₂³ + 𝒪(‖x‖⁷)

julia> pol < 1  # this is correct
false

julia> q = 1 - 0.5*(x+y)^4
 1.0 - 0.5 x₁⁴ - 2.0 x₁³ x₂ - 3.0 x₁² x₂² - 2.0 x₁ x₂³ - 0.5 x₂⁴ + 𝒪(‖x‖⁷)

julia> q < 1 # again, this should return true...
false

[PerezHz]

Just noted that even though #323 was merged to fix this issue, it remains open, is this intended?

[Li-shiyue]

There is author issue:

x = set_variables("x", order=100)[1] Base.AbstractFloat(x::TaylorN{T}) where {T} = x.coeffs[1].coeffs[1] _₂F₁(1.0,2.0,3.0,(x.-1)) 0.6137056388801092

However,I can only get the first coeff, How to fix it? Can anybody help me

[lbenet]

Just noted that even though #323 was merged to fix this issue, it remains open, is this intended?

I guess I left it open just to check if everything was stable. Once said this, i think this can be closed...

[lbenet]

There is author issue:

x = set_variables("x", order=100)[1] Base.AbstractFloat(x::TaylorN{T}) where {T} = x.coeffs[1].coeffs[1] _₂F₁(1.0,2.0,3.0,(x.-1)) 0.6137056388801092

However,I can only get the first coeff, How to fix it? Can anybody help me

I guess this comment is related to #335, isn't it? If so, we should better keep the discussion there...

[Li-shiyue]

There is author issue: x = set_variables("x", order=100)[1] Base.AbstractFloat(x::TaylorN{T}) where {T} = x.coeffs[1].coeffs[1] _₂F₁(1.0,2.0,3.0,(x.-1)) 0.6137056388801092 However,I can only get the first coeff, How to fix it? Can anybody help me

I guess this comment is related to #335, isn't it? If so, we should better keep the discussion there...

Yes, It really confuses me