Solving Higher Order Matrix Polynomials

Hi, guys:

I wonder if there is some numerical method to solve a general matrix polynomial in the form: I+A1*X+A2*X^2+...+Aq*X^q=0. where X is supposed to be a matrix with the same dimension as As. A1 to Aq are k-by-k square matrices.

Any hint or reference is highly appreciated.

[Walter Roberson]

To cross-check, q is a positive integer?

[Jing]

Yes, q is a positive integer. q=1 is a simple question, but cases for q>=2 are much more complicated. If k=1, which is just a scalar case, we can always covert it to a matrix 1st order equation and the solution is simple. However, to stack matrix up to reduce the order does not work for the multivariate case.

Answer 1

[Teja Muppirala]

If k and q are not too large, one idea is to try to solve it as an optimization problem.

"Which elements of X will yield the smallest norm of the residual"

(you will have to save this as a function):

function solvepoly

k = 4;
qmax = 3;

for q = 1:qmax
    A(:,:,q) = randn(k); %Make some random A matrices
end

[xf,fval] = fminunc(@doCost,zeros(k))

      function COST = doCost(x)
          COST = eye(k);
          for q = 1:qmax
              COST = COST + A(:,:,q)*x^q;
          end
          COST = sum(COST(:).^2);
      end
end

[Jing]

Teja:

Thanks for the input. This searching mechanism can give a good approximation for the solution. However, I wonder if there are more efficient ways, since it is pretty easy to calculate the eigenvalues of this polynomial. The main problem is how to use these eigenvalues efficiently. In addition, There will be multiple solutions for sure, how can way make sure the searching mechanism to converge.