Resolvent matrix 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.

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?
