Substituting matrices in polynomials

[mahrud]

This issue is inspired by this Stack Exchange question that I randomly ran into.

Would be great if the following worked:

S = QQ[t]
m = random(ZZ^2, ZZ^2)
f = det(m - t*m^0) -- characteristic polynomial
-- By Cayley-Hamilton theorem we should have:
0 == f[m] -- error: expected a matrix with 1 row
0 == sub(f, t => m) -- error: expected substitution values to be ring elements or numbers

The easiest workaround I can think of is still too complicated:

0 == sum(pairs reverse flatten entries last coefficients(f, Variables => t), (i, c) -> m^i * c)

But more generally, if I have a polynomial in some number of variables, it shouldn't be too difficult to substitute matrices (or even ring maps) for the variables or return an error if the matrices are not compatible in some way.

Related: #2863, #2792

[pzinn]

with univariate polynomials, I think this is a good idea and should be implemented. not sure about multivariate polynomials, since the matrices you feed may not commute and then it's ambiguous how to take products.

[mahrud]

I main had linear combinations of variables in mind, but you could use a multivariate free algebra instead.