Excessive compilation times with v2.0

[andreasvarga]

Moving my project MatrixPencils to using Polynomials v0.2 I am facing the following problem: the execution times of the tests increased from about 16 minutes (typical run earlier before switching to v2.0) to around 50-60 Minutes when using Julia 1.2-1.5 and once again around 16 Minutes for Julia nightly version. I wonder if in general such excessive increases in compilation times are o be expected? Is any way to improve compilation times?

[jverzani]

Hmm, that is a huge regression.

Putting the symbol into the type forces a new compilation everytime there is a new symbol used, but I wouldn't have expected that to be this big an issue, as I would think there aren't many symbols used beyond :x, :z, :t,:u or so.

I'm guessing, but are some types not fully realized under the new setup? (By that, I'm guessing if there are any Polynomial{Float64,X} where {X} type things once you form the matrices.) That is something we can address through hook or crook. As seen in our last discussion, the promotion machinery isn't as automatic as I had thought it would be, but hopefully we can sort out the main issues with just a few new methods.) As well, if your structs don't also have the X in the type indication, then they will be slowup. Another slowup might be the type instability that comes about now by allowing constructs like p + q to ignore the symbol for constant terms, but I did test for regressions there (though using @btime and not @time so compilation isn't included.

Hopefully we can sort this out. The point of this was to make it easier and have more predictable behaviour.

[andreasvarga]

I followed your recommendation and used everywhere AbstractVecOrMat{<:Polynomial} and never used Polynomial{T,X} type. Could this be the cause of the slow down?

My last update of MatrixPencils needed 52m for 1.2, 50m for 1.3, 46m for 1.4, 39m for 1.5 and 16m for Julia nightly (I assume 1.6). It seems that improvemets in Julia can also improve this issue.

[jverzani]

I think that isn't an issue, as it only controls dispatch. I was thinking more of the structs (or types) you create and the fields therein. But in browsing the code, I only found 1 struct, which didn't have that problem. (I was expecting some more, so likely missed them).

Is there a smaller example I could look at (the whole test suite is a bit too broad) where there is a regression in execution time (presumably measured with @.` and not @.` so we can capture compile time extras)?

On Sun, Mar 14, 2021 at 10:38 AM Andreas Varga @.***> wrote:

I followed your recommendation and used everywhere AbstractVecOrMat{<:Polynomial} and never used Polynomial{T,X} type. Could this be the cause of the slow down?

My last update of MatrixPencils https://github.com/andreasvarga/MatrixPencils.jl needed 52m for 1.2, 50m for 1.3, 46m for 1.4, 39m for 1.5 and 16m for Julia nightly (I assume 1.6). It seems that improvemets in Julia can also improve this issue.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/JuliaMath/Polynomials.jl/issues/314#issuecomment-798918746, or unsubscribe https://github.com/notifications/unsubscribe-auth/AADG6TB7OMNFVS3YESEWB7TTDTC4TANCNFSM4ZE5UK3Q .

-- John Verzani Department of Mathematics College of Staten Island, CUNY @.***

[jverzani]

BTW, I didn't look that carefully, but did you have to do anything to get compatability with 2.0? If there aren't structs involved or extracting the symbol with p.var then there shouldn't have been anything special needed.

On Sun, Mar 14, 2021 at 12:16 PM John Verzani @.***> wrote:

I think that isn't an issue, as it only controls dispatch. I was thinking more of the structs (or types) you create and the fields therein. But in browsing the code, I only found 1 struct, which didn't have that problem. (I was expecting some more, so likely missed them).

Is there a smaller example I could look at (the whole test suite is a bit too broad) where there is a regression in execution time (presumably measured with @.` and not @.` so we can capture compile time extras)?

On Sun, Mar 14, 2021 at 10:38 AM Andreas Varga @.***> wrote:

I followed your recommendation and used everywhere AbstractVecOrMat{<:Polynomial} and never used Polynomial{T,X} type. Could this be the cause of the slow down?

My last update of MatrixPencils https://github.com/andreasvarga/MatrixPencils.jl needed 52m for 1.2, 50m for 1.3, 46m for 1.4, 39m for 1.5 and 16m for Julia nightly (I assume 1.6). It seems that improvemets in Julia can also improve this issue.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/JuliaMath/Polynomials.jl/issues/314#issuecomment-798918746, or unsubscribe https://github.com/notifications/unsubscribe-auth/AADG6TB7OMNFVS3YESEWB7TTDTC4TANCNFSM4ZE5UK3Q .

-- John Verzani Department of Mathematics College of Staten Island, CUNY @.***

-- John Verzani Department of Mathematics College of Staten Island, CUNY @.***

[andreasvarga]

λ = Polynomial([0,1],:λ) @time info, iz, nfp, ip = rmkstruct(λ, 1, fast = true)

john verzani @.***> schrieb am So., 14. März 2021, 17:22:

BTW, I didn't look that carefully, but did you have to do anything to get compatability with 2.0? If there aren't structs involved or extracting the symbol with p.var then there shouldn't have been anything special needed.

On Sun, Mar 14, 2021 at 12:16 PM John Verzani @.***> wrote:

I think that isn't an issue, as it only controls dispatch. I was thinking more of the structs (or types) you create and the fields therein. But in browsing the code, I only found 1 struct, which didn't have that problem. (I was expecting some more, so likely missed them).

Is there a smaller example I could look at (the whole test suite is a bit too broad) where there is a regression in execution time (presumably measured with @.` and not @.` so we can capture compile time extras)?

On Sun, Mar 14, 2021 at 10:38 AM Andreas Varga @.***> wrote:

I followed your recommendation and used everywhere AbstractVecOrMat{<:Polynomial} and never used Polynomial{T,X} type. Could this be the cause of the slow down?

My last update of MatrixPencils https://github.com/andreasvarga/MatrixPencils.jl needed 52m for 1.2, 50m for 1.3, 46m for 1.4, 39m for 1.5 and 16m for Julia nightly (I assume 1.6). It seems that improvemets in Julia can also improve this issue.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub < https://github.com/JuliaMath/Polynomials.jl/issues/314#issuecomment-798918746 , or unsubscribe < https://github.com/notifications/unsubscribe-auth/AADG6TB7OMNFVS3YESEWB7TTDTC4TANCNFSM4ZE5UK3Q

.

-- John Verzani Department of Mathematics College of Staten Island, CUNY @.***

-- John Verzani Department of Mathematics College of Staten Island, CUNY @.***

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaMath/Polynomials.jl/issues/314#issuecomment-798935528, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALJDHEAJK36VYCUXWVSDOR3TDTPEBANCNFSM4ZE5UK3Q .

[andreasvarga]

I am back to the old figures with a single change (replaced two lines of the commented code with the new ones):

# function poly2pm(PM::Polynomial; grade::Union{Int,Missing} = missing) 
#    d = degree(PM)
#    T = eltype(PM)
#    ismissing(grade) ? k = d+1 : k = max(d,grade)+1
#    k == 0 && (return zeros(T,1,1,1))
#    P = zeros(T,1,1,k)
#    d < 0 || (P[1,1,1:d+1] = coeffs(PM))
#    return P      
# end
function poly2pm(PM::Polynomial{T,X}; grade::Union{Int,Missing} = missing) where {T<:Number,X}
   d = degree(PM)
   ismissing(grade) ? k = d+1 : k = max(d,grade)+1
   k == 0 && (return zeros(T,1,1,1))
   P = zeros(T,1,1,k)
   d < 0 || (P[1,1,1:d+1] = coeffs(PM))
   return P      
end

Do you have any explanation for this ?

[jverzani]

I wouldn't have guessed this as the difference is basically related to dispatch, not run time. I assume it has to do with the T in the signature being known at compile time, rather than the look up at run time. Likely you would get the same with just Polynomial{T} ... where {T}. (I wouldn't hard code T<:Number, in a perfect world I'd like to remove that restriction...)

Anyways, I haven't had a time to look further here for the overall regression. Am I to believe this one change fixed all of that?

[andreasvarga]

Apparently yes. I am using now Polynomial{T} ... where {T}.

[jverzani]

Well that's good news. Now to tweak the promotion rules in #312 to get something reasonable.

[andreasvarga]

I will close this issue, with some remarks. I would be pleased to understand the implications of using types as Polynomial, <:Polynomial, Polynomial{T} and Polynomial{T,X} on the compilation and execution times (apparently there are no observable implications on execution times). It would also be important to know which form is best suited when implementing software involving polynomial matrices. For my purposes, the form Polynomial{T,X} (not used now in MatrixPencils) would allow to implicitly test that several matrices have the same variable X. However, I prefered to use the more generic form with <:Polynomial, which allowed to keep (most of) code unchanged when transiting to v2.0. Surprisingly, in one case, when only a polynomial was involved (not a matrix), I needed to use the form Polynomial{T} instead Polynomial to come back to reasonable compilation times, from about 50 minutes to 16 minutes (the usual value for this package from previous versions).

[jverzani]

Yes, my understanding is that you should have been able to keep p::Polynomial as the type signature of a function to control dispatch, the only real task type signatures do in a function. In that particular function that was modified, you want the element type, so using p::Polynomial{T} where {T} avoids the run-time lookup of the element type. Though this shouldn't be costly at all. I wouldn't have expected specifying the parameter to lead to fewer compilations, so the slow down is just surprising to me. My mental model is that the function gets compiled for every new type. So Polynomial{Float64, :x} and Polynomial{Float64, :y} arguments would require an additional compilation. (This was not the case when the symbol was a field and not in the type domain, so some compilation related slowdown is expected.)

However, were you to create a struct, you would want to have the full specification in the field so that you have a concrete type. For example, a struct for rational functions might look like:

struct RatFunction{T,X}
  num::Polynomial{T,X}
  den::Polynomial{T,X}
end

Just using Polynomial{T} would leave the fields parameterized as Polynomial{T,X} where {X}, hence abstract, not concrete.

Hope that helps. Others in the Julia community could no doubt explain this better.

[andreasvarga]

Thanks. I will do in this way in the second package DescriptorSystems which supports rational functions. And, if you decide at some time to implement a pure rational functions package, I would make the transition once again.

[jverzani]

Should I be looking to put the code from DescriptorSystems for rational functions into Polynomials? It might fit right in.

[andreasvarga]

In DescriptorSystems I would like to support a special rational function, which I call RationalTransferFunction. As it is now, this structure is

struct RationalTransferFunction{T} <: AbstractRationalTransferFunction
    num::Polynomial{T}        # numerator polynomial
    den::Polynomial{T}        # denominator polynomial
    Ts::Union{Real,Nothing}   # sampling time (0 - continuous-time, -1 or > 0 - discrete-time, nothing - none)
    function RationalTransferFunction{T}(num::Polynomial{T,X}, den::Polynomial{T,X}, Ts::Union{Real,Nothing}) where T <: Number where X
        length(num) > 1 && length(den) > 1 && Polynomials.indeterminate(num) != Polynomials.indeterminate(den) && 
              error("The numerator and denominator polynomials must have the same variable")
        if all(den == zero(den))
            error("Cannot create a rational function with zero denominator")
        elseif all(num == zero(num))
            # The numerator is zero, make the denominator 1
            den = one(den)
        end
        # Validate sampling time
        isnothing(Ts) || Ts >= 0 || Ts == -1 || 
             error("Ts must be either a positive number, 0 (continuous system), or -1 (unspecified)")
        new{T}(num, den, Ts)
    end
end

Probably, after having a pure type RationalFunction, this could have a simpler form

struct RationalTransferFunction{T,X} <: AbstractRationalTransferFunction
    r::RationalFunction{T,X}     # rational function
    Ts::Union{Real,Nothing}   # sampling time (0 - continuous-time, -1 or > 0 - discrete-time, nothing - none)
end

completed with suitable constructors. I guess, I still would have to support all kind of operations with rational matrices + sampling time (which probably would be simpler than using two polynomials). So, independently of this, I think it is a good idea to have a rational function structure (which would also fit with what you use for Pade approximations). I wonder how much can be reused from the already existing RationalFunctions package, which is apparently not further developed.

[andreasvarga]

I am about to finish the transition to Polynomials v2.0. The times are somewhat better for version 1.5 and 1.6 (20-25 minutes for a test run), but over 1 hour for vesion 1.2.

An issue which I encountered today is handling rational matrices with zero dimensions. I observed that building the first construction, e.g., [R R0] withR0 having zero columns can take a very, very long time to compile it. I tried to figure out why, the cause could be the following construction of such an R0:

julia> using DescriptorSystems

julia> R0 = rtf.(rand(3,0))
3×0 Array{RationalTransferFunction{Float64,X} where X,2}

I was not able to figure out which elementary constructor is behind this operation. My guess, is that the broadcast function provides automatically such a result. I tried this also for polynomials, where I obtained:

julia> using Polynomials

julia> P0 = Polynomial.(rand(3,0))
3×0 Array{Any,2}

which is also strange.

[jverzani]

I have to admit, I'm not sure what goes wrong with Polynomial.(rand(3,0)). It is different from both map(Polynomial, rand(3,0)) and [Polynomial(x) for x ∈ rand(3,0)], which keep the type.

I see this:

julia> Meta.@lower Polynomial.(rand(3,0))
:($(Expr(:thunk, CodeInfo(
    @ none within `top-level scope'
1 ─ %1 = rand(3, 0)
│   %2 = Base.broadcasted(Polynomial, %1)
│   %3 = Base.materialize(%2)
└──      return %3
))))

But just don't know what broadcasted defaults to for zero-length arguments.