Power series do not distinguish between different precisions!

[Sebastian Oehms]

Dear all,

obviously Sage does not distinguish between precisions of power series:

sage: R.<t> = PowerSeriesRing(ZZ)
sage: O(t).is_zero()
True
sage: O(t) == O(t**2)
True

Similarly for p-adics:

sage: O(3).is_zero()
True
sage: O(3) == O(3^2)
True
sage:

This seems to be explicitly intended:

   def __nonzero__(self):
        """
        Return True if this power series is not equal to 0.

        EXAMPLES::

            sage: R.<q> = ZZ[[ ]]; R
            Power Series Ring in q over Integer Ring
            sage: f = 1 + 3*q + O(q^10)
            sage: f.is_zero()
            False
            sage: (0 + O(q^2)).is_zero()
            True
            sage: R(0).is_zero()
            True
            sage: (0 + O(q^1000)).is_zero()
            True
        """
        return not not self.polynomial()

Why?

[David Roe]

We only store power series and p-adics to finite precision, so it's impossible to tell whether the mathematical objects underlying two elements in Sage are actually equal.  In practice, the goal of a p-adic or power series computation is to determine some result approximately, since that's all that's possible.  In this context, the notion of equality implemented in Sage is the right choice.  For example, suppose you want to do some long computation and determine whether you get 1 at the end.  At the end of your computation, you find the answer 1 + O(t^33).  I think it's more useful to say that this is equal to 1 than that it's not.

I'm not sure what kind of equality you would imagine, but if we defined equality by something like "equal absolute precision and equal value modulo that absolute precision" then you no longer have additive and multiplicative inverses.  For example, if x = 1 + O(3^18) then there is no value of y so that x+y == 0.  That's pretty devastating for thinking about algebraic operations.

If you want to be more careful with the precision of your equality testing you can use the `is_equal_to` method.

sage: x = 1 + O(3^14)

sage: y = 1 + O(3^16)

sage: x.is_equal_to(y, absprec=10)

True

sage: x.is_equal_to(y, absprec=15)
PrecisionError (most recent call last):

...

PrecisionError: Elements not known to enough precision

Finally, I'll note that there a bunch of precision models in Sage for p-adics (they all apply in a power series context as well but haven't been implemented there because it's a lot of work).  The tutorial is still missing some of the newer ones (floating point precision and the two lattice precisions), but you can pick out the details from the reference manual if you're willing to wade through a bunch of text.

David

--
You received this message because you are subscribed to the Google Groups "sage-devel" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sage-devel/10981ee0-efea-4885-850c-7e4dcd0e4d32%40googlegroups.com.

[Nils Bruin]

On Wednesday, January 15, 2020 at 7:09:29 AM UTC-8, David Roe wrote:

I'm not sure what kind of equality you would imagine, but if we defined equality by something like "equal absolute precision and equal value modulo that absolute precision" then you no longer have additive and multiplicative inverses.  For example, if x = 1 + O(3^18) then there is no value of y so that x+y == 0.  That's pretty devastating for thinking about algebraic operations.

There is also the more computer-science type equality: will these two numbers lead to the same results if substituted in otherwise identical computations. This can be relevant for testing algorithmic aspects of what you're doing. That means asking: do the elements have the same valuation, digits, and precision (i.e., are they represented by the same p-adic disc). This requires some kind of normalization of 0-centered discs to be well-defined.

I think at the moment you can get a tuple capturing this information by something like

   x.__reduce__()[1][1:][1:]

but this is obviously not very robust if pickle formats change. I would have expected canonicalized construction parameters to be available efficiently somehow, but I didn't find an interface routine for it easily. (I don't have an immediate need for them and I could work around it if I did, but perhaps it's worth including -- and perhaps it's already there and I didn't find it).

[Sebastian Oehms]

For example, suppose you want to do some long computation and determine whether you get 1 at the end.  At the end of your computation, you find the answer 1 + O(t^33).  I think it's more useful to say that this is equal to 1 than that it's not.

Sorry, that I was not clear enought about which precision I was referring to. I didn't mean the default precision of the ring, as you seem to do, but the individual precision of elements which can be set using the add_bigoh method (or global O function). In my examples the default precision has been:

sage: R
Power Series Ring in t over Integer Ring
sage: R.default_prec()
20

But the individual precisions of the elements where below that.

I'm not sure what kind of equality you would imagine, but if we defined equality by something like "equal absolute precision and equal value modulo that absolute precision" then you no longer have additive and multiplicative inverses.  For example, if x = 1 + O(3^18) then there is no value of y so that x+y == 0.  That's pretty devastating for thinking about algebraic operations.

My expectation is that two power series are considered to be equal if there polynomials agree, but also their individual precisions as long as one of them is below the default precision of the ring. Accordingly the `is_zero` method should not return True if the individual precision is below the default precision.

If you want to be more careful with the precision of your equality testing you can use the `is_equal_to` method.

But this doesn't help for the following irritating answer (which was the original issue pointed out to me by Bill Hart):

sage: R.<t>=PowerSeriesRing(ZZ)
sage: S.<u>=PowerSeriesRing(R)
sage: 0 + O(t) + O(u)
O(u^1)

BTW: There is an inconsistency in the representation string which doesn't look very nice:

sage: O(t)
O(t^1)
sage: 1+O(t)
1 + O(t)

Best,

Sebastian

To unsubscribe from this group and stop receiving emails from it, send an email to sage-...@googlegroups.com.

[Nils Bruin]

 
On Wednesday, January 15, 2020 at 12:45:10 PM UTC-8, Sebastian Oehms wrote:

My expectation is that two power series are considered to be equal if there polynomials agree, but also their individual precisions as long as one of them is below the default precision of the ring. Accordingly the `is_zero` method should not return True if the individual precision is below the default precision.

The default precision is just that: a default. It's used to get a finitely presented result in cases where a computation could normally lead to a non-finite one. The setting of the default doesn't affect the mathematical properties of the elements. In particular, whether elements are equal or not is determined by their direct properties, not by the default set in the parent ring. Therefore, your expectation is not compatible with the design (as you discovered). But perhaps by seeing how it was designed and why, you'll warm to the design decisions made :-).

The arithmetic model used is basically one of "interval" or "ball" arithmetic: In terms of actual power series rings, an element like 1+t+O(t^2) stands for the set of all power series with initial part 1+t. In that model, "equality" gets translated to "have non-empty intersection".

This model has the advantage that (sqrt(1+t)^2 -1)/t == 1 returns true, as one would expect mathematically.

(insisting equality up to the default precision would have to lead to false, because the arithmetic on the LHS leads to precision loss)

 

But this doesn't help for the following irritating answer (which was the original issue pointed out to me by Bill Hart):

sage: R.<t>=PowerSeriesRing(ZZ)
sage: S.<u>=PowerSeriesRing(R)
sage: 0 + O(t) + O(u)
O(u^1)

Indeed, but I think the actual bug there is: Power series rings over non-exact base rings don't work properly. This is a problem that you should expect: algebra algorithms normally designed to work with exact base rings will have problems with non-exact ones, and usually for the reason you see here: to work properly with power series over exact rings, you eliminate (exact) zeros, but the inexact zeros that you encounter in power series rings should probably be kept for precision tracking reasons (also the ones up to default precision!). But keeping them around indiscriminately would lead to impractical results very quickly. Special care is needed. You might try working with ZZ[['t','u']] instead.

(and indeed, power series rings with p-adic coefficients need a whole set of extra care too -- and usually convergence properties on the coefficients anyway to make meaningful arithmetic possible)

[Sebastian Oehms]

... This model has the advantage that (sqrt(1+t)^2 -1)/t == 1 returns true, as one would expect mathematically.  (insisting equality up to the default precision would have to lead to false, because the arithmetic on the LHS leads to precision loss)

Your example is convincing and I'm not going to mistrust the ball design. The truth is that we can't say whether O(t) is zero or not. But since the answer is a boolean we have to make a decision (giving an "exact" answer for an "inexact" question). Independend on how this decision is made cases of irritations will be left. For example, if a user defines R.<t>=PowerSeriesRing(ZZ, default_prec=100) and types in p + O(21) with some polynomial p of degree 20 than he would not expect `p + O(21) == p` to be True since that looks as if we have lost something (similar as in Bill Harts example).

Indeed, but I think the actual bug there is: Power series rings over non-exact base rings don't work properly. This is a problem that you should expect: algebra algorithms normally designed to work with exact base rings will have problems with non-exact ones, and usually for the reason you see here: to work properly with power series over exact rings, you eliminate (exact) zeros, but the inexact zeros that you encounter in power series rings should probably be kept for precision tracking reasons (also the ones up to default precision!). But keeping them around indiscriminately would lead to impractical results very quickly. Special care is needed. You might try working with ZZ[['t','u']] instead.

So, what should we do than? Print a warning if somebody tries to define them (and suggest to use multivariate power series in cases that would lead to an alternative)?

[Marc Mezzarobba]

[re-posting a reply from a week ago that apparently did not go through
because gmane was moving]

Nils Bruin wrote:
> This model has the advantage that (sqrt(1+t)^2 -1)/t == 1 returns
> true, as one would expect mathematically.

Do you mean it has the advantage that cos(sin(tan(t^2)) -
tan(sin(t^2))) == 1 returns True, as one would expect mathematically?
;-)

Joking aside, if I hadn't be hit by that issue before, I would also
expect to be able to trust equality more than that. And I would
interpret the presence of an explicit O() term as a strong indication
that inexact series won't be considered equal to anything.

Perhaps more importantly, I find the fact that series and p-adics (but
not intervals and balls) are doing that problematic for writing generic
code. Suppose that a is a symbolic expression, or an element of any
other parent where inequality is not decidable. Would you expect
a.is_zero() to return True whenever Sage is unable to prove that a is
nonzero? If not, what can code written for generic coefficient rings do
to work with both expressions and series?

Regarding power series in particular, the structure where
cos(sin(tan(t^2)) - tan(sin(t^2))) == 1 exists and makes perfect sense,
but it's the ring of polynomials mod t^20, not the ring of power series
with precision 20.

--
Marc

[Nils Bruin]

 On Wednesday, January 22, 2020 at 5:07:18 AM UTC-8, Marc Mezzarobba wrote:

 

Perhaps more importantly, I find the fact that series and p-adics (but
not intervals and balls) are doing that problematic for writing generic
code. Suppose that a is a symbolic expression, or an element of any
other parent where inequality is not decidable. Would you expect
a.is_zero() to return True whenever Sage is unable to prove that a is
nonzero? If not, what can code written for generic coefficient rings do
to work with both expressions and series?

It has the big advantage that if you compute a root of a polynomial, then evaluating the polynomial at that root will give you a result that is equal to zero. That's extremely useful in for instance writing routines that test local solvability. It also means that if you assume you've started out with enough precision to distinguish any non-equal elements that occur in your computation, then equality gives you a reliable result.

Under the assumption you are starting out with sufficient excess precision to start, I think it's actually the model that makes generic algorithms still work: even gaussian elimination works properly "in the limit"; without some numerical stability tricks it likely just needs ridiculous amounts of initial precision.

I don't think that for any non-trivial applications generic code is going to perform properly for both exact and non-exact base rings, but actually the equality choice that has been made, allows the heuristic approach: start with very high precision or double the precision a couple of times and see if results stabilize. This can be extremely useful in getting some experimental verification of conjectures.

 

Regarding power series in particular, the structure where
cos(sin(tan(t^2)) - tan(sin(t^2))) == 1 exists and makes perfect sense,
but it's the ring of polynomials mod t^20, not the ring of power series
with precision 20.

--
Marc

[Sebastian Oehms]

On Wednesday, January 22, 2020 at 2:07:18 PM UTC+1, Marc Mezzarobba wrote:

[re-posting a reply from a week ago that apparently did not go through
because gmane was moving]

Nils Bruin wrote:
> This model has the advantage that (sqrt(1+t)^2 -1)/t == 1 returns
> true, as one would expect mathematically.

Do you mean it has the advantage that cos(sin(tan(t^2)) -
tan(sin(t^2))) == 1 returns True, as one would expect mathematically?
;-)

This is a good example, too. Here, even comparison using the predefined default precision wouldn't help. As I said before: It's impossible to avoid such wrong answers. The question is how we can make that transparent to the user and that question is independent on what realization of comparison has more advantages.

Let me say just this without thinking about the work and annoyance for users practice that would cause: The cleanest solution would be to print a warning the first time a user applies `==` or `!=` to elements of non-exact rings (once per session or once per instance of such a ring). In addition a second pair of comparison operators (using infix operator) should be supplied giving according results but without warning. These new operators should be used in the doctests.

But comparison is only one part of the irritations. The original issue (printing `O(u) + O(t)` as `O(u)`) has another reason, which can be explained by a look at the behavior in Nemo. Nemo does comparison in the same way as we do:

julia> R, t = PowerSeriesRing(ZZ, 100, "t")
(Univariate power series ring in t over Integer Ring, t+O(t^101))
julia> O(t) == 0
true

But there seems to be a principle difference: Sage allows the coexistence of exact (prec = infinity) and non-exact (prec = integer) elements, while Nemo seems to do not. Thus, in Nemo you have:

julia> R(0)
0+O(t^100)

while in Sage:

sage: R.<t> = PowerSeriesRing(ZZ, default_prec=100)
sage: R.zero()
0
sage: R.zero().prec()
+Infinity

Accordingly:

julia> S, u = PowerSeriesRing(R, 20, "u")
(Univariate power series ring in u over Univariate power series ring in t over Integer Ring, u+O(u^21))

julia> O(u) + O(t)
0+O(t^100)+O(u^1)

but:

sage: S.<u>  = PowerSeriesRing(R)
sage: O(u) + O(t)
O(u^1)

On the one hand, this coexistence is reasonable, since mathematically polynomials are special power series. But since there are principle differences between exact and non-exact rings in computer algebra, this is another source of irritations.
 

[Marc Mezzarobba]

Nils Bruin wrote:
> I don't think that for any non-trivial applications generic code is
> going to perform properly for both exact and non-exact base rings, but
> actually the equality choice that has been made, allows the heuristic
> approach: start with very high precision or double the precision a
> couple of times and see if results stabilize. This can be extremely
> useful in getting some experimental verification of conjectures.

I agree that it can be convenient, but I still don't see how you make
that work in a system where some other objects have a more rigid notion
of equality. The point being that you can't really choose a convention
independently for, say, power series and symbolic expressions; things
like generic polynomial and matrices need equality to behave roughly the
same for both.

(The real answer is probably that you need both to distinguish "is
certainly nonzero" for "may be nonzero" everywhere, *and* to carefully
choose between several "equality" predicates, e.g., "is equal", "is
approximately equal", "is formally/syntactically equal". But that's not
really an option for Sage.)

--
Marc