CAS Best Practices

Currently in Symbolics.jl (0.10.0) any numeric type (the default symbol) automatically transforms x - x to 0. If I understand right, @oscarbenjamin and @carette are both saying it should stay as x - x until explicitly evaluated with simplify()/evaluate()/etc?

1 Like

I was warning that x-x is tricky. It’s actually fairly safe to transform it to 0*x. For one, that’s type-safe, in the sense that x is still around, so you know the result’s type properly. (Polymorphic 0 is a major pain.)

If your numeric types include infinities, then this is required. If your numeric types only contain finite values, and will never contain anything other than classical numbers (vectors, matrices and series all have non-trivial claims to be kinds of numbers too) like the reals, complex, etc, that might be safe. But it depends very strongly on your meaning of ‘number’.

5 Likes

It’s zero(x) or zero(T) in Julia

6 Likes

Okay I see why the Nemo and Flint timings are so similar: Nemo just uses Flint for the calculation. This is what I meant by speed ultimately being determined by the lower-level libraries.

Yes, that is what I would do if I was making a new CAS. Not just x - x and also x/x. Also y + x would not automatically canonicalise to x + y. Functions like cos(0) should stay as cos(0). This is about both performance and usability. Many of the things that people want to use a CAS for are not necessarily about computing anything: they just want to be able to manipulate expressions in a controlled way.

In SymPy you can do that with cos(0, evaluate=False) and Add(y, x, evaluate=False). Using evaluate=False is a kludge though that tries to repair the damage caused by automatic evaluation retrospectively. It can never fully work because evaluation needs to default to False from the foundations. Any kind of automatic evaluation once introduced is very hard to remove because future code comes to depend on it.

How in Symbolics.jl do I create the expression cos(0) or sin(pi) or even sqrt(2)?

11 Likes

This doesn’t work if x is an array — zero(x) * x will throw if x is not a square matrix. In the context of rings it is fine, of course, though it still won’t be equivalent to x-x for square matrices that may contain non-finite values.

I guess you could do something like zeroscalar(x) * x where

zeroscalar(x) = zero(x)
zeroscalar(::Type{<:AbstractArray{T}}) where {T} = zeroscalar(T)
zeroscalar(x::AbstractArray) = zeroscalar(eltype(x))
8 Likes

15 posts were split to a new topic: Symbolic sqrt(2) and sin(pi)

julia> using Symbolics
julia> Term(sqrt,[2])
sqrt(2)

I don’t get why you’d work harder than that. You can build any expression to be lazily like that, or use @register on others.

I thought the reason for it was because the poster wanted “a zero in the type of x”:

In which case trying to get 0 * x as they posted is an X-Y problem: just writing it as zero(x) is more direct and means what they actually meant using Base Julia functionality.

4 Likes

I think you’re right, I was mistaken. Term(-, [x,x]) produces x - x which seems like what’s desired. Why doesn’t x - x produce the same result? For @syms x::Any, x-x errors, and for @syms x::Num, x-x returns 0::Float64.

0 * x is not the same as zero(x) if x is not finite (NaN or ±Inf). That’s why x - x can be simplified to 0*x (for an appropriate 0) but not to zero(x) in general, and simplifying x-x is the context here.

3 Likes

As the originator of this thread, it is my hope and request that we ( that means you Julians :slight_smile: ) return to being about “CAS Best Practices”. That there are different considerations held by CAS experts is helpful information. Let it be, so we move forward. Thank you all for participating.

5 Likes

If you want it exact, don’t ask for the Float64 result?

julia> ex = Term(sqrt,[big(2)])
sqrt(2)

julia> eval(toexpr(ex))
1.414213562373095048801688724209698078569671875376948073176679737990732478462102

It just follows Julia’s semantics.

I see. Yeah there’s a more interesting question to be answered here, which is what kind of number is the default symbol? The default is @variables x::Real, specifically not ::AbstractFloat. Since these rules are all type-based, I think simplifying x-x -> 0 on ::Real is the semantics you’d expect, while x-x -> 0 * x should be what happens on ::AbstractFloat. This keeps two algebras distinct: real numbers obey commutative and associative laws, they are a field, etc. while floating point numbers have a different set of rules.

The reason for keeping this difference is that the set of transformations that keep exact floating point semantics is already in LLVM, so I don’t think there’s too much to gain there. But it would be a good opt-in: “make simplifications only do what’s allowed on floating point numbers” instead of “make simplifications only do what’s allowed on real numbers”. The latter seems to be what people normally want on models, or at least out of a CAS, though there are times when exact floating point semantics need to be followed (in which case we have @register for someone to register a function and keep some calculation intact). What I want to do is add a call to Herbie, so that in many cases we aren’t just applying rules that change floating point behavior but actually improve it, so even further down the path of “keep floating point separate”.

4 Likes

This makes sense on its own, but I do not know how to square it with AbstractFloat <: Real. This is a question related to (or exactly about?) co/contra/invariance, and how the subtyping relation shows up, if at all, in figuring out what rules apply when. Is the following a bug? Or is it a necessary consequence of the rule defined on Real?

julia> y = SymbolicUtils.Sym{AbstractFloat}(:y)
y

julia> y - y
0
1 Like

I would say yes, it’s a missing specialization.

You have subtypes which override defaults of their supertype all of the time. AbstractArray for example has many array types which change underlying behaviors of broadcast, indexing, etc.

But subtypes should not override the defaults (or “fallbacks”) of their supertype “willy nilly”. There should be some kind of consistency with the fallback. e.g.:

julia> typeof(1:4)
UnitRange{Int64}

julia> typeof(1:4) <: AbstractArray
true

julia> (1:4) .+ 5
6:9

julia> collect(1:4) .+ 5
4-element Vector{Int64}:
 6
 7
 8
 9

julia> collect((1:4) .+ 5) == collect(1:4) .+ 5  # consistency
true

If broadcasted addition on UnitRange{Int64} were defined in such a way that the “# consistency” property did not hold, that would be bad for generic programming!

Does that badness apply in this CAS setting? What consistency should exist across rules that are defined on types that have a subtype relationship among them? I am not sure. That’s what the question is.

1 Like

That’s a good question. I would assume that for number types you’d want to have those simplifications by default, and choose on a subset of AbstractFloats to disable them. I could see someone going the other way and having @variables x default to SymbolicReal <: Real and define those rules on SymbolicReal. It’s worth a discussion on the SymbolicUtils.jl repo. Open an issue.

2 Likes

There is no problem whatsoever with defining the square root function (sqrt or \sqrt{}) to be the principal square root. This is the standard approach, as described in the Wikipedia article on square roots and as implemented in Julia, Python, Mathematica, Maple, Maxima, etc. Insisting that we define \sqrt{9} to be \pm3 is a pathway to confusion.

I would think that if you want to show that \sqrt{x^2} and x are unequal as expressions, you’d want to point out the example of x = -3. That way the argument does not depend on a non-standard interpretation of the square root symbol.

2 Likes

From a conventional programming language point-of-view, yes,

From a logico-symbolic point-of-view where
y = sqrt(x) means that x == y^2 , y = sqrt(9) means that 9 == y^2
asking Maple to solve(y^2 = 9, y) gives (3, -3)

A good deal of the back-and-forth is frisson between Numerics and Symbolics (both are math).

We are working toward a best (or a very good) case implementation approach that allows us to glean important benefits that each perspective brings, sometimes one better than the other and vice versa.

My point is exactly that the square root function is, in all common usage, distinct from the solve function. If we insist on conflating the two, how do we interpret an expression like \sqrt{2} + \sqrt{3}? It’s got to be a four-element set, I suppose. Here we’re already verging on chaos, wishing we had a concise way to say “the positive square root” when we want to.

I’m comfortable granting the point that there are many interesting loci of interplay between the demands of symbolic and numeric computation, some of which may require creative thinking about how common concepts need to be treated differently. But I remain skeptical about the specific proposal of abandoning the square root function (which is implied by an assertion like \sqrt{9} = \pm 3).

4 Likes

We are not going to abandon any of these print(names(Base.Math)). I understand how that impression may have been conveyed, inadvertently. What is being considered is how to respect the numerical math functions without abandoning the symbolic/theoretical math functions that may carry the same names in ordinary discourse (hence the tangling of notions). It is my guess that we will find a way to distinguish the desired context (symbolic, numeric, some of each) for any occurrence of the named function – there are multiple avenues for this. The simplest way to distinguish the Julia numerical square_root from the more constructively mathematical symbolic square_root would be use the name from Base.Math as the numerical function and to use its titlecase spelling as the symbolic function sqrt(2) !== Sqrt(2). That leaves out a mixed version (maybe Sqrt is mixed, and purely symbolic is a context, as with let, and do blocks).

see next post by @DNF for one of the more familiar approaches within Julia

Doesn’t that run a bit counter to the general Julia philosophy? Could one instead leverage multiple dispatch with sqrt(2) vs sqrt(Symbolic(2)) (just making up names here)?

Disclaimer: I know nothing about CASs in general, and hardly anything about Symbolics.jl, so this is a question rather than a suggestion.