compute the limit of a rational function

Question

Suppose I have a rational function defined by ($s$ complex) $$ f(s)=w^T s(sI-Q)^{-1} v $$ for nonzero column vectors $w,v$ and a (large) square matrix $Q$. Further assume that $Q$ is singular and that it is known that $f(0)=\lim_{s\rightarrow 0} f(s)$ exists finite. Clearly, for any $s$ that is not an eigenvalue of $Q$, one can compute $f(s)$ by $$ f(s)=s \cdot w^T y \text{ where $y$ is the unique solution of } (sI-Q)y=v $$ which is quite efficient and numerically stable.

Question: is there a similar way of computing $f(0)$, that is (ideally) by solving a single linear system?

[Context. In my problem, $f(s)$ arises as $sF(s)$, where $F(s)$ is the Laplace transform of a certain (exponential) function $g(t)$. The result is needed as a method to compute $f(0)=\lim_{t\rightarrow \infty} g(t)$, which is the final value theorem for Laplace transform.]

Answer 1

Amazingly, linear algebra is remarkably unstable in this respect. Here is an example that appeared recently in my research. (It is not exactly like your function, but it is very close: you can easily make an appropriate change of variables).

Let $\Theta$ be the direct sum of two copies of $$ \begin{bmatrix}0&b+1\\b&c\end{bmatrix}, $$ and let $\kappa:=[2b+1,c]^t\oplus[0,1]^t$ (for some integers $b\ne0,-1$ and $c\ne0$). Consider $H(s):=(1-s)(\Theta-s^{-1}\Theta^*)$. Then the function $f(s):=\langle H(s)^{-1}\kappa,\kappa\rangle$ is identically $0$ whenever defined (and, hence, all limits are $0$), whereas $\kappa$ is not in the image of $H(s)$ for $s=(1+1/b)^{\pm1}$, so these values (that should also be $0$) cannot be "computed directly".

Answer 2

Consider the Jordan canonical form $Q = U J U^{-1}$ for $Q$. Thus $w^T Q w = \sum_B w_B^T B v_B$ for Jordan blocks $B$, where $w_B$ and $v_B$ are the entries of $U^T w$ and $U^{-1} v$ corresponding to the block. For a Jordan block $B$ of size $m$ corresponding to the eigenvalue $0$, $s (sI - B)^{-1}$ is upper triangular with $1$ on the main diagonal and $1/s^j$ on the $j$'th diagonal above that. If it is known (for particular $w$ and $v$) that a finite limit exists, all the terms in $1/s^j$ in $w_B^T s (sI-B)^{-1} v_B$ must cancel, leaving the same result that you would have for $B = 0$. On the other hand, for a Jordan block corresponding to a nonzero eigenvalue, the limit is $0$. Thus if $P$ is the projection on the generalized eigenspace for eigenvalue $0$, the limit, if it exists, is $w^T P v$.

$P$ can be obtained as a polynomial $g(A)$, where if the largest Jordan block for eigenvalue $0$ has size $m_0$, $g(0) = 1$ and $g^{(j)}(0) = 0$ for $j < m_0$, while for every other eigenvalue $\lambda$ with largest Jordan block of size $m_\lambda$, $g(\lambda) = \ldots = g^{(m_\lambda-1)}(\lambda) = 0$.

Of course, this can't be numerically stable: eigenvalues that are exactly $0$ and those that are close to $0$ will produce very different results.