Constructing equivalent algebraic expressions for matrix equations

Question

I have an expression involving matrices, of the form:

$$f(k)=x^T A_k^{-1}A x$$

where $x$ is a $1\times N$ vector, $A_k = A + k I$ and $A$ is an $N\times N$ matrix ($A_k$ is invertible for all $k$) and $k>0$. It is known that $f(k)$ is real and monotonic increasing, but nothing more. I need to further analyze the behaviour of $f(k)$ and the simplest way would be to plugin my matrices and plot it. However, this becomes computationally intensive for large $N$ as it involves calculating the inverses repeatedly.

One thought I had was to create a one-to-one map from $f(k)$ to an equivalent algebraic expression $g(k)$, which is easy to evaluate and investigate analytically. Indeed, if the expression had involved traces of inverses, one could've used the Stieltjes transform for some insight, but that doesn't seem likely here.

My question is: Are there general approaches/references I can look at to learn how to tackle such problems?

Answer 1

Another approach would be an improved version of what René Pannekoek suggests in a comment: compute the Schur form (not the Jordan form, which is hopelessly numerically unstable) of $A=QTQ^H$; get for free the Schur form of $A+kI=Q(T+kI)Q^H$; set $y=Q^H x$; rewrite as $f(x)=y^H(T+kI)^{-1}Ty$. This lowers your computational cost for evaluating the function in $K$ different points from $O(n^3K)$ to $O(n^3+n^2K)$ (precompute the Schur form, and do for each point one matrix product with $Q$, one with $T$, one triangular system to solve with $T+kI$).

I did not check the details, but I assume that Dima Pasechnik's method is $O(n^4+Kn)$, since you need $O(n)$ matmuls to precompute the Padé approximation. If you combine the two ideas (first reduce to Schur form, then compute a Padé approximation of $g(y)=y^H(T+kI)^{-1}Ty$), I think that you can get $O(n^3+nK)$.

I am afraid that the Padé part could be more numerically unstable though. Since $f(k)$ has poles in the eigenvalues of $A$, I guess that the denominator of the order-$n$ Padé approximant is the characteristic polynomial of $A$, and generally working with characteristic polynomials isn't a good idea numerically.

Answer 2

Just use the usual Cramer rule for the inverse of $A+kI$. This shows that $f(k)$ is a rational function of $k$ with degrees of numerator and denominator bounded by $N$, which you can find, e.g. numerically, by Pade approximation.

More sophisticated techniques, used in combinatorics, for tackling this kind of stuff, are called "transfer matrices". See e.g. Vol.I of R. Stanley's book, where in the 2nd edition, available online, it starts on p.573.

EDIT: while numerically this might not be optimal, at least this tells you upfront what kind of function you are dealing with. And if you're only interested in the asymptotic behaviour of $f(k)$ for large $k$, then with a bit of luck you might not even need to compute the full rational form expression for $f(k)$.