polyvalm2 - A faster matrix polynomial evaluator

polyvalm2 evaluates a polynomial with a matrix argument faster than the MATLAB function polyvalm.

polyvalm2 evaluates a polynomial with a square matrix argument faster than the MATLAB built-in functions polyvalm or mpower.

Y = polyvalm2(P,X), when P is a vector of length N+1 whose elements are the coefficients of a polynomial, is the value of the polynomial evaluated with matrix argument X. X must be a square matrix.

Y = P(1)*X^N + P(2)*X^(N-1) + ... + P(N)*X + P(N+1)*I

Class support for inputs P, X:
float: double, single

The polyvalm2 speed improvements come from the following:

1) The MATLAB built-in function polyvalm uses Horner's method. polyvalm2 uses a binary decomposition of the matrix powers to do the calculation more efficiently, reducing the total number of matrix multiplies used to calculate the answer.

2) polyvalm calculates the product of a scalar times a matrix as the product of diag(scalar*ones(etc))*matrix ... i.e. it does a matrix multiply. polyvalm2 will calculate this more efficiently as the simple product scalar*matrix.

3) polyvalm does all of the calculations shown above, even if the coefficient P(i) is zero. polyvalm2 does not do calculations for P(i) coefficients that are zero.

4) polyvalm converts sparse matrix inputs into full matrices to do the calculations, whereas polyvalm2 keeps the intermediate calculations and the answer sparse.

The trade-off is that polyvalm2 uses more memory for intermediate variables than polyvalm, so for very large matrices polyvalm2 can run out of memory. In these cases polyvalm2 will abandon the efficient calculation method and just call the built-in polyvalm. For sparse matrix inputs, however, polyvalm2 will typically be more memory efficient than the MATLAB polyvalm function.

Caution: Since polyvalm2 uses different calculations to form the matrix powers, the end result may not match polyvalm exactly. Also, for the case where only the leading coefficient is non-zero, polyvalm2 may not match mpower exactly. But the answer will be just as accurate. And if there are inf's or NaN's involved, then the end result will not, in general, match polyvalm or mpower results. This should not be a great drawback to using polyvalm2, however, since even the MATLAB built-in functions polyvalm and mpower will not match each other in these cases. (By reordering the calculations, the NaN's propagate differently)