Resolvent of a triangular 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.
 A: The Leverrier-Faddeev algorithm for a triangular matrix ? No kidding !
Here $N(x)=Adjoint(xI-A)$. It suffices to inverse a triangular matrix; cf. this algorithm, the complexity of which, is $\approx n^3/3$:
http://www.iaeng.org/publication/WCE2012/WCE2012_pp100-102.pdf
Yet, here, we multiply polynomials in $K[x]$ and not only elements in $K$.
EDIT 2: answer to Michele. 1. Of course, the complexity of the above cited method is $\approx n^4/3$ mult. in $K$ (using FFT for the product of polynomials).


*About the instability, let $A\in M_n(\mathbb{Z})$, where the $a_{i,j}$ have $k$ digits; then some coefficients of the entries of $N(x)$ have almost $kn$ digits (see $N(x)[1,n]$). You have the same problem when you calculate the gcd of $2$ polynomials over $\mathbb{Z}$. If you work with finite-precision arithmetic, then of course there is a real risk.

*Michele, about the complexity of L-F, you are right - I completely forgot the formula $N(x)=\sum_{i=1}^n x^{n-i}N_i$ where $N_i=AN_{i-1}+p_{i-1}I$  and $p_i=coeff(p_A(x),n-i)$. I tested this method and, indeed, it is fast and good if you work in extended precision. Otherwise you can proceed as follows (cf. point 4).

*Choose $n$ values $(x_i)_i$ s.t. the matrices $x_iI-A$ are well-conditioned and calculate the $((x_i-A)^{-1})_i$, that is the $(N(x_i))_i$ (a complexity in $n^4/3$ again). According to point 3., it remains to solve a linear system in the unknowns $(N_i)_i$; the essential task is to inverse a Vandermonde matrix, that can be done in $O(n^3)$ and even in $O(n^2)$ in a numerically stable way (cf. http://www.unix.eng.ua.edu/~japalmore/papers/vandermonde3.pdf ).
