Resolvent matrix

Question

Suppose $A$ is a triangular matrix. What is the most efficient known algorithm to compute the polynomial (in $x$) matrix $(xI-A)^{-1}$?

Of course, $(xI-A)^{-1}= N(x)/p_A(x)$, where $p_A$ is the characteristic polynomial of $A$, which is easy to compute once we know an eigendecomposition of $A$. But what about $N(x)$?

I am aware of the Leverrier-Fadeev algorithm, which requires $O(n^4)$ operations if $A$ is $n\times n$. Moreover, it makes use of power iteration, which can lead to numerical instability.

Answer 1

Since $A$ is triangular, you may try to first diagonalize if, $A=PDP^{-1}$. You already know what the eigenvalues of $A$ are. Then, $$(xI-A)^{-1} = (P(xI-D)P^{-1})^{-1} = P(xI-D)^{-1}P^{-1}$$ and $(xI-D)^{-1}$ is trivial to calculate. Does this help?

Answer 2

$(xI-A)^{-1}=Adj(xI-A)/p(x)$ where $p(x)$ is the (known) characteristic polynomial of $A$. On the other hand, $Adj(sI-A)=\Delta p(s,A)$ where $\Delta p(s,t)=\dfrac{p(s)-p(t)}{s-t}$ (a symmetric polynomial).

For example, let $n=3,p(x)=x^3-2x^2-3x+2$; then

$\Delta p(s,t)=\dfrac{s^3-t^3-2s^2+2t^2-3s+3t}{s-t}=s^2+st+t^2-2s-2t-3$.

Finally $Adj(xI-A)=(x^2-2x-3)I+(x-2)A+A^2$.

EDIT. The resolvent is essentially used to calculate $f(A)$ when $f$ is an holomorphic function.

Here $(xI-A)^{-1}=\sum_j a_j(x)/p(x) A^j$ where $a_j$ is a polynomial.

Then $f(A)=1/2i\pi\sum_j (\int_{\gamma} f(x)a_j(x)/p(x)dx) A^j$ where $\gamma$ is a ad hoc curve. Consequently, we must only calculate $n$ numerical integrals.

In particular, when $\lambda_1,\cdots,\lambda_n$, the roots of $p$, are simple:

$f(A)=\sum_j(\sum_k f(\lambda_k)a_j(\lambda_k)/p'(\lambda_k)) A^j$.

Anyway, the complexity of the calculation is that of the computation of the $(A^j)_j$, that is here $\sim n^4/6$ (of course, if the calculation of the Jordan form of $A$ is stable, then we can do the job with complexity about $30 n³$).