evalpoly for matrix polynomials

[stevengj]

Closes #1161.

Includes a generic out-of-place fallback implementation as well as an in-place evalpoly! method; the latter is called by default only for X::StridedMatrix{<:Number} for now, to be conservative. (e.g. it's much trickier to allocate the in-place buffers correctly for matrices of matrices, which we don't support very well, and I'm not sure this would work well for sparse matrices.) More optimized cases can always be added later.

[mcabbott]

In addition to X::StridedMatrix{<:Number}, doesn't the in-place version need to widen based on the coefficients' types? Some failure cases:

julia> evalpoly([1 2; 3 4], (5, 6))  # error without this PR
2×2 Matrix{Int64}:
 11  12
 18  29

julia> evalpoly([1 2; 3 4], (5, 6.7))
ERROR: InexactError: Int64(6.7)
Stacktrace:
...
 [11] evalpoly!(Y::Matrix{Int64}, X::Matrix{Int64}, p::Tuple{Int64, Float64})
    @ Main ./REPL[4]:19

julia> evalpoly([1.0 2.0; 3.0 4.0], (2, 3+im, 4))
ERROR: InexactError: Float64(7.0 + 1.0im)

[stevengj]

Argh, this is a PITA to get right for heterogeneous tuples of coefficients, especially for a tuple of dimensionful coefficients. (e.g. suppose that X is in meters, and the coefficients p are in (s, s/m, s/m^2) — this is dimensionfully correct with an evalpoly result in seconds, but you can't simply promote the elements of p).

What is a good way to compute the correct element type of the result matrix?

For p a tuple of numbers, I'm tempted to do typeof(evalpoly(zero(eltype(X)), p)), but that will fail if eltype(X) is abstract.

One option would be to just punt on evalpoly! for now. Another option would be to only support homogeneous tuples.

[mcabbott]

Ah now I see some logic for p[begin] which my examples happened to avoid.

Only accepting p::NTuple (and concretely typed vectors) seems fine really. And perhaps even StridedMatrix{<:AbstractFloat} with real coefficients (to avoid my integer example), sim. complex with real-or-complex?

Easy to add methods to send your own fancy types to evalpoly! if you want something not covered by this.

[stevengj]

I tightened up the type signatures for dispatch to evalpoly! and it fixes your examples, along with abstractly typed arrays. (Also added tests.)

[mcabbott]

This is a little harder to fool, but:

julia> evalpoly([1 2; 3 4], Number[5.5, 6])
ERROR: InexactError: Int64(11.5)

julia> evalpoly([1.0 2.0; 3.0 4.0], Number[2, 3, 4+im])
ERROR: InexactError: Float64(4.0 + 1.0im)

[codecov]

Codecov Report

Attention: Patch coverage is 87.75510% with 6 lines in your changes missing coverage. Please review.

Project coverage is 92.02%. Comparing base (61e444d) to head (68b3f0f).

Files with missing lines	Patch %	Lines
src/evalpoly.jl	87.75%	6 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1163      +/-   ##
==========================================
- Coverage   92.03%   92.02%   -0.02%     
==========================================
  Files          34       35       +1     
  Lines       15501    15550      +49     
==========================================
+ Hits        14267    14310      +43     
- Misses       1234     1240       +6     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[LilithHafner]

I think it makes sense to dispatch to _evalpoly in the case of vectors with non-concrete eltype. We could do a reducion with promotion but I don't think it's worth trying particularly hard to optimize that case.

[ViralBShah]

Is this good to merge (assuming tests pass)?

[LilithHafner]

@mcabbott's comment is still not addressed. That said, those cases already error.

[ViralBShah]

Would it be better to merge this if it is at least addressing some of the issues?

[LilithHafner]

CI failure is in matmul so probably unrelated. If so, I think it's likely reasonable to merge this even with @mcabbott's unresolved comment. Do you agree, @mcabbott?

[mcabbott]

I have no strong feelings... my examples above invented to break this do seem unlikely to show up in the wild.