Possible future internal improvements

[rdeits]

Continuing the discussion from over at #28 ...

I've spent some time working on two different possible alternative designs for multivariate polynomials:

1. TypedPolynomials, which stores the variable names in the type parameters of each Monomial
2. StaticPolynomials, which stores the variable names and exponents in the type parameters of each Monomial

They're both still extremely rough (only addition and multiplication are tested). Both are using some new type system features from Julia 0.6, so they're both 0.6-only.

A consequence of both of these designs is that a variable's name is its identity. That is, any two variables with the name :x are the same variable. Otherwise code like:

for i in 1:10
  @polyvar x
  x^2
end

would have to compile new methods for each loop iteration and would be type-unstable. That's a big change from MultivariatePolynomials, where every variable has a unique ID. It's similar, however, to the way MultiPoly.jl works, and I think I'm OK with it.

TypedPolynomials

In this design, we have, roughly:

immutable Variable{Name} end

immutable Monomial{Variables}
  exponents::NTuple
end

immutable Term{T, M <: Monomial}
  coefficient::T
  monomial::M
end

immutable Polynomial{T <: Term, V <: AbstractVector{T}}
  terms::V
end

StaticPolynomials

This design has:

immutable Variable{Name} end

immutable Power{Variable}
  exponent::Int
end

immutable Monomial{Powers} end

immutable Term{Mono, T}
    coefficient::T
end

immutable Polynomial{T <: Tuple}
  terms::T
end

Those two designs have led to some different design choices. In TypedPolynomials, I've kept the standard Vector of Terms in each polynomial. That's necessary because when adding two polynomials, we can't infer at compile time how many terms the result will have. In StaticPolynomials, however, we can infer the number of terms, but each term is a different type, so the terms are stored in a heterogeneous tuple of Terms.

Advantages of both packages:

• speed: constructing and evaluating polynomials is much faster. Also, promotion from variable, monomial, or term up to Polynomial is essentially free and does not allocate
• type-stable substitution: since the variables live in the type, Julia can infer when substitution should result in a polynomial or a plain value. I've also found a different syntax for substitution:

p = x^2 + y
subs(p, x=>1, y=>2)

The advantage of this syntax is that it doesn't require allocating a vector of variables and a vector of values. In fact, complete substitution of all variables in a polynomial now allocates no memory at all.

Notable differences:

• TypedPolynomials was easier to write, and I have more features working there
• StaticPolynomials is much faster for operations like constructing a polynomial, and almost never allocates any memory. On the other hand, it's also likely to be slower to compile, harder to read, and it's been harder to actually get working
• TypedPolynomials allocates memory when adding two polynomials, which makes long expressions like 3x + 2y + 1 + x + z^3 inefficient. Some macro magic similar to what JuMP uses could help here.
• In StaticPolynomials, the return type of x^n depends on the value of n. This is less of an issue in julia 0.6 where x^n is lowered to x^Val{n}, but it still means code like [x^i for i in 0:3] is inherently type-unstable.
  • The flip side of that is that StaticPolynomials could potentially infer that the type of differentiate(3x, x) is Int, while TypedPolynomials could not.

Performance

As an example, constructing the polynomial x^2 + y + x * x + 3 * x * y + x * y on my laptop with Julia 0.6 takes:

MultivariatePolynomials:

memory estimate:  11.47 KiB
  allocs estimate:  197
  --------------
  minimum time:     135.708 μs (0.00% GC)

TypedPolynomials

  memory estimate:  1.05 KiB
  allocs estimate:  12
  --------------
  minimum time:     478.665 ns (0.00% GC)

StaticPolynomials

memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     1.974 ns (0.00% GC)

And substitution takes:

MultivariatePolynomials:

julia> @benchmark subs($p, [1, 2], [$x, $y])
BenchmarkTools.Trial:
  memory estimate:  3.14 KiB
  allocs estimate:  56
  --------------
  minimum time:     71.974 μs (0.00% GC)
  median time:      73.575 μs (0.00% GC)
  mean time:        85.366 μs (0.00% GC)
  maximum time:     342.221 μs (0.00% GC)

TypedPolynomials

julia> @benchmark subs($p, $x=>1, $y=>2)
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     17.546 ns (0.00% GC)
  median time:      17.559 ns (0.00% GC)
  mean time:        18.482 ns (0.00% GC)
  maximum time:     69.566 ns (0.00% GC)

StaticPolynomials

julia> @benchmark $p($x=>1, $y=>2)
BenchmarkTools.Trial:
  memory estimate:  80 bytes
  allocs estimate:  4
  --------------
  minimum time:     84.476 ns (0.00% GC)
  median time:      87.768 ns (0.00% GC)
  mean time:        109.203 ns (7.50% GC)
  maximum time:     4.456 μs (92.01% GC)

Summary

I'm currently leaning towards TypedPolynomials. It's easier to understand and develop, and I think its performance is good enough for what we're doing. But I'm interested in getting your opinion on both of these designs. I have a demo of TypedPolynomials here: https://github.com/rdeits/TypedPolynomials.jl/blob/master/demo.ipynb which shows off some of the basic features.

[blegat]

I am impressed, good work ! This might have a big impact in the performance of PolyJuMP and SumOfSquares :)
Why do you add monomial in the fields of TypedPolynomials:

immutable Term{T, M <: Monomial}
  coefficient::T
  monomial::M
end

couldn't you just do:

immutable Term{T, M <: Monomial}
  coefficient::T
end

I haven't yet checked the details of the implementation but it seems to me that we could keep several implementations. They would share common functions thanks to abstract types like PolynomialLike.

We would could have

@polyvar x :Static

that creates a variable x of type StaticVariable{:x}.
Then the type of x^2 will be StaticMonomial{(Power{StaticVariable{:x}()}(2),)}, and so on.
That is, you could do everything with :Static and then you could just change :Static to :Typed in the @polyvar line and everything will change accordingly!
If the user does

@polyvar x

then it would default to :Typed.
What do you think ? It would be insanely great to be able to switch between implementation so easily :-P

[rdeits]

Why do you add monomial in the fields of TypedPolynomials

I think you maybe mean StaticPolynomials (aside: these names are not great. I'm open to more clear names)? For TypedPolynomials, the monomial instance holds the exponents, so we need to keep the instance as a member of the Term. For StaticPolynomials you're totally right: in fact, what you suggested is actually what's in the code; I just put the wrong thing in the issue above.

it seems to me that we could keep several implementations

Wow, that's true. I hadn't thought of that, but I think you're right. These two implementations can probably share a lot of code. OK, I'll see if it's possible.

[joehuchette]

This looks very impressive, @rdeits. Two quick questions/comments:

• When you say that StaticPolynomials is "harder to read," you mean the implementation, correct? AFAICT the user-facing syntax is identical.
• It'd be interesting to test out the compile time for StaticPolynomials on longer expression. While it's intuitive that the compilation time would be (much) longer, I've been pleasantly surprised by these types of things in the past.

[rdeits]

When you say that StaticPolynomials is "harder to read," you mean the implementation, correct? AFAICT the user-facing syntax is identical.

Yes, that's right. I'm still working the implementation out, though, so it might turn out to be cleaner than my current version.

It'd be interesting to test out the compile time for StaticPolynomials on longer expression. While it's intuitive that the compilation time would be (much) longer, I've been pleasantly surprised by these types of things in the past.

Absolutely. I'm doing some reworking of TypedPolynomials now to try to generalize as much as possible, then I'll try to wire up StaticPolynomials through the same general functions so that we can get a fair test.

[blegat]

Yes, indeed, I meant StaticPolynomials.
I have thought little bit more about it.
I kind of prefer TypedPolynomials over StaticPolynomials even though the latter can detect it when the sum of two terms is a term (for this reason, at some point I thought about doing something like StaticPolynomials but never actually did anything).

I also like the fact that with TypedPolynomials, when doing computations, the type rapidly converges. Indeed, it rapidly goes to something like Polynomial{Term{Monomial{(x, y, z)}, Float64}, Vector{...}}. With StaticPolynomials, a variable might constantly change its type. For instance, with the expression sum([x*y, x^2, y^2]) in StaticPolynomials you won't be able to promote the types of the x*y, x^2, y^2 to a common concrete type. sum won't be able to infer its return type. And in the best case it will do something like p = 0; for m in array; p += m; done in which p will constantly change its type. However, in TypedPolyniomals, the type of the array will be Monomial{(x, y)} and it will correctly infer that the type of the sum will be Polynomial{Term{Monomial{(x, y)}, Int}, Vector{...}}, will initialize p with zero(Polynomial{...}) and the sum will be type stable.

One other concern I have with StaticPolynomials is that it cannot see that terms cancel like in x + x^2 - x and it will have to store zero terms which is something we can avoid with TypedPolynomials. Moreover, we cannot even write simplify(::StaticPolynomials) that will remove zero terms since it will be type unstable. If someone wants to do that we might add a conversion TypedPolynomial(::StaticPolynomial) which will be type stable.

In fact it seems to me that TypedPolynomials is a strict improvement over the current implementation.
When I think about it, it seems that TypedPolynomals was what I was trying to implement the whole time ! This whole thing with MonomialVector is the fact that I want to promote the monomials to the same TypedMonomial type since what I do in MonomialVector is putting the variables in common. From x + y, I extract the vector of variable [x, y] and the MonomialVector contains the exponents [[1, 0], [0, 1]].
This is the reason a Polynomial is not a vector of Term because each Term has a field vars with [x, y] and I do not want to maintain such copy.

In short, I think we need to choose between

• replacing the current implementation with TypedPolynomials which is the perfect implementation the universe have been looking for since its creation billions of years ago.
• Keeping both TypedPolynomials and StaticPolynomials and being able to choose which one to use when doing @polyvar.

About the naming, if we remove the current implementation we might rename Typed* into * ans Static* in to Static* or just S*

[rdeits]

@blegat thanks for all the comments 😄 . I agree with all of your concerns about StaticPolynomials. One more weirdness of StaticPolynomials is that if you want to create a concretely typed vector containing [1, x, x^2], then the type of that vector has to be:

Vector{Tuple{Term{Int, Monomial{()}}, Term{Int, Monomial{(x)}}, Term{Int, Monomial{(x^2}}}

which is crazy. Essentially, StaticPolynomials only ever makes sense if you always store tuples, not vectors.

However, I think the good news is that it's not much more work to get StaticPolynomials functioning alongside TypedPolynomials. Providing both implementations would be neat. I think my plan is:

• try to get StaticPolynomials up to basic functionality (it's close)
• meanwhile, implement as much of the Typed* and Static* functionality in an abstract way that works for both types
  • I've done a lot of this already, and it's had the additional benefit of clarifying the Typed* code too
• try to get Typed* all the way up to full compatibility with Multivariate* and get the Multivariate* tests all passing
  • this will require implementing non-commutative variables, which I've kept in the back of my mind but haven't implemented yet
• try to swap out Typed* as the implementation inside Multivariate*

And yes, I think Polynomial and SPolynomial/StaticPolynomial are good name choices.

[blegat]

I like the plan, thanks for the hard work ! Let me know if you need any help.
Not that I am not strongly attached to my abstract types TermContainer, ... so you can totally change the abstract structure if needed.
Also the fact that you can now do {T, S{T}} with the types parametrization with Julia v0.6 may allow to change a few things inside MultivariatePolynomials :-P

[blegat]

I was thinking about something. Since TypedPolynomials and StaticPolynomials use @generated, this prevents compilation. This might be a reason for keeping an alternative DynamicPolynomials that might look like the current implementation but that we try to make as close as TypedPolynomials as possible but without using @generated functions. What do you think ? Would you say that this limitation with compilation and @generated functions is only temporary ?

[rdeits]

I'm not sure I know what you mean. v0.6 brought in some tighter restrictions on what can be done inside the body of a @generated function, but I think precompilation should work regardless.

[blegat]

I was talking about static compilation, see here.

[rdeits]

Ah, understood. Yes, then that would be a reasonable thing to do. I'll keep trying to move as much behavior as possible to the abstract types.

[rdeits]

Good news! I've managed to eliminate all of the @generated functions in TypedPolynomials: JuliaAlgebra/TypedPolynomials.jl#3 Performance is slightly worse (maybe 1.5x slower than before), but I suspect that can be fixed.

Of course, StaticPolynomials is a whole other story. Doing that without @generated seems unlikely.

[rleegates]

First of all, thank you for your great work, the development of this package will probably carry our (admittedly very MultiPoly-dependent) code into the future! I am starting to migrate our code to MultivariatePolynomials in the course of migrating to julia v0.6, however, I would like to circumvent the performance gap (ideally gain performance) which appears (#3) to exist without the improvements made in TypedPolynomials. Is there a branch in MultivariatePolynomials which integrates TypedPolynomials at this time?

[rdeits]

Sorry, not yet. It's definitely still on my radar, but real life has a habit of intruding.

The primary missing features are basic calculus (pretty easy) and non-commutative variables (possibly not so easy).

[blegat]

Finally \o/ !