Package for matrix polynomials

[fph]

What package do you recommend to work with matrix polynomials, i.e., expressions of the form A(x) = A_0 + A_1 x + \dots + A_d x^d?

I will use mainly the coefficients to do linear algebra, so I am interested in a package that keeps the coefficients A_i stored in contiguous arrays: a Matrix{Polynomial} wouldn’t be convenient.

After a quick search I found these packages, but they all have some drawbacks:

• GitHub - neveritt/PolynomialMatrices.jl: Polynomial matrices in Julia : but it was last updated 3 years ago and doesn’t even precompile without errors on 1.9.3.
• GitHub - jagot/MatrixPolynomials.jl : despite the encouraging name, it deals with a completely different concept (polynomials with scalar coefficients evaluated in a matrix, as in Krylov subspaces and matrix functions).
• GitHub - andreasvarga/MatrixPencils.jl: Matrix pencil manipulations using Julia : while it has many useful decompositions implemented, it does not contain good Julian syntactic sugar to work with matrix polynomials: there isn’t even a convenience type for matrix polynomials to define A * B to compute A(x)B(x), for instance.

[stevengj]

@fph

What package do you recommend to work with matrix polynomials, i.e., expressions of the form A(x) = A_0 + A_1 x + \dots + A_d x^d?

What operations do you want to perform? The Polynomials.jl package can do simple computations with these already, since from its perspective matrices are just a particular kind of coefficient:

julia> using Polynomials

julia> p = Polynomial([rand(3,3) for i=1:4])
Polynomial([0.8608493451549502 0.40249646273782136 0.29736298082105506; 0.3033857369382127 0.360125250924879 0.1391618742877383; 0.6172709803058753 0.09508887526836374 0.43790723011832033] + [0.9596676501415019 0.43430610756543553 0.2846821519890206; 0.8739974216258194 0.6927972973867174 0.3170401151329463; 0.09750151351768044 0.7350982960890662 0.6266149502835732]*x + [0.0925452448501407 0.44676221806701466 0.585415953263665; 0.25485246297681297 0.3320962411878218 0.9943801889841738; 0.27532271123727503 0.996163909889251 0.940717353711509]*x^2 + [0.35043298755900976 0.7187106937852754 0.44758569352667943; 0.5436421830854086 0.19671391692418616 0.8245177117682397; 0.43835896242201156 0.013564423506638446 0.4054034609556135]*x^3)

julia> p(3)
3×3 Matrix{Float64}:
 14.0345  25.1315  18.505
 19.8974  10.7387  32.3017
 15.2234  11.6321  21.7301

Of course, if if you want to compute A(x) v where v is a vector, then it is inefficient to add all of the matrices first — better to multiply each term by v and accumulate them in-place. I don’t think Polynomials.jl has a specialized method for this but it would be easy to add.

But if you need more specialized linear-algebra functions (which?) they would be better off in another package (maybe built-in top of Polynomials).

[fph]

@stevengj

What operations do you want to perform? The Polynomials.jl package can do simple computations with these alre

Thanks, that should work! I had incorrectly assumed that Polynomials.jl wouldn’t work for matrix coefficients, as Google took me to an older version of the docs which had a constraint T <: Number.

I don’t need complicated operations, mostly it’s enough to have a a container type that supports sums / poly-poly products / pretty-printing, so Polynomials.jl should do the job.