0
$\begingroup$

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.]

$\endgroup$
3
  • $\begingroup$ I'm not sure if it helps, but if we consider diagonal $Q$ (e.g. if we diagonalize), then $\lim_{s \to 0} s(sI-Q)^{-1}$ sends the zero diagonal entries to $1$ and the nonzero diagonal entries to $0$. $\endgroup$
    – user13113
    Mar 11, 2015 at 23:28
  • $\begingroup$ It may be simpler to write $sI = (sI-Q) + Q$ and use the alternative form $$ f(s) = w^T v + (w^T Q) (sI-Q)^{-1} v $$ $\endgroup$
    – user13113
    Mar 14, 2015 at 19:46
  • $\begingroup$ @Hurkyl. Hello, could you explain why this form is simpler? $\endgroup$
    – Michele
    Mar 15, 2015 at 22:38

2 Answers 2

1
$\begingroup$

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".

$\endgroup$
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.