dynare/doc/sylvester.tex

\documentclass[11pt,a4paper]{article}

\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{amsthm}

\usepackage{fullpage}

\begin{document}
\title{Solution of Specialized Sylvester Equation}
\author{Ondra Kamenik}
\date{2004}
\maketitle

\renewcommand{\vec}{\mathop{\hbox{vec}}}
\newcommand{\otimesi}{\mathop{\overset{i}{\otimes}}}
\newcommand{\iF}{\,^i\!F}
\newcommand{\imF}{\,^{i-1}\!F}
\newcommand{\solvi}{\mathbf{solv1}}
\newcommand{\solvii}{\mathbf{solv2}}
\newcommand{\solviip}{\mathbf{solv2p}}

\newtheorem{lemma}{Lemma}

Given the following matrix equation
$$AX+BX\left(\otimesi C\right)=D,$$ where $A$ is regular $n\times n$
matrix, $X$ is $n\times m^i$ matrix of unknowns, $B$ is singular
$n\times n$ matrix, $C$ is $m\times m$ regular matrix with
$|\beta(C)|<1$ (i.e. modulus of largest eigenvalue is less than one),
$i$ is an order of Kronecker product, and finally $D$ is $n\times m^i$
matrix.

First we multiply the equation from the left by $A^{-1}$ to obtain:
$$X+A^{-1}BX\left(\otimesi C\right)=A^{-1}D$$
Then we find real Schur decomposition $K=UA^{-1}BU^T$, and
$F=VCV^T$. The equation can be written as
$$UX\left(\otimesi V^T\right) +
KUX\left(\otimesi V^T\right)\left(\otimesi F\right) =
UA^{-1}D\left(\otimesi V^T\right)$$
This can be rewritten as
$$Y + KY\left(\otimesi F\right)=\widehat{D},$$
and vectorized
$$\left(I+\otimesi F^T\otimes K\right)\vec(Y)=\vec(\widehat{D})$$
Let $\iF$ denote $\otimesi F^T$ for the rest of the text.

\section{Preliminary results}

\begin{lemma}
\label{sylv:first-lemma}
For any $n\times n$ matrix $A$ and $\beta_1\beta_2>0$, if there is
exactly one solution of
$$\left(I_2\otimes I_n +
\begin{pmatrix}\alpha&\beta_1\cr-\beta_2&\alpha\end{pmatrix}
\otimes A\right)\begin{pmatrix} x_1\cr x_2\end{pmatrix} = 
\begin{pmatrix} d_1\cr d_2\end{pmatrix},$$
then it can be obtained as solution of
\begin{align*}
\left(I_n + 2\alpha A+(\alpha^2+\beta^2)A^2\right)x_1 & = 
\widehat{d_1}\\
\left(I_n + 2\alpha A+(\alpha^2+\beta^2)A^2\right)x_2 & = 
\widehat{d_2}
\end{align*}
where $\beta=\sqrt{\beta1\beta2}$, and 
$$
\begin{pmatrix} \widehat{d_1}\cr\widehat{d_2}\end{pmatrix} =
\left(I_2\otimes I_n +
\begin{pmatrix}\alpha&-\beta_1\cr\beta_2&\alpha\end{pmatrix}
\otimes A\right)\begin{pmatrix} d_1\cr d_2\end{pmatrix}$$
\end{lemma}

\begin{proof} Since 
$$
\begin{pmatrix} \alpha&\beta_1\cr-\beta_2&\alpha\end{pmatrix}
\begin{pmatrix} \alpha&-\beta_1\cr\beta_2&\alpha\end{pmatrix}
=
\begin{pmatrix} \alpha&-\beta_1\cr\beta_2&\alpha\end{pmatrix}
\begin{pmatrix} \alpha&\beta_1\cr-\beta_2&\alpha\end{pmatrix}
=
\begin{pmatrix} \alpha^2+\beta^2&0\cr 0&\alpha^2+\beta^2\end{pmatrix},
$$
it is easy to see that if the equation is multiplied by
$$I_2\otimes I_n +
\begin{pmatrix}\alpha&-\beta_1\cr\beta_2&\alpha\end{pmatrix}
\otimes A$$
we obtain the result. We only need to prove that the matrix is
regular. But this is clear because matrix
$$\begin{pmatrix} \alpha&-\beta_1\cr\beta_2&\alpha\end{pmatrix}$$
collapses an eigenvalue of $A$ to $-1$ iff the matrix
$$\begin{pmatrix} \alpha&\beta_1\cr-\beta_2&\alpha\end{pmatrix}$$
does.
\end{proof}

\begin{lemma}
\label{sylv:second-lemma}
For any $n\times n$ matrix $A$ and $\delta_1\delta_2>0$, if there is
exactly one solution of
$$\left(I_2\otimes I_n +
2\alpha\begin{pmatrix}\gamma&\delta_1\cr-\delta_2&\gamma\end{pmatrix}\otimes A
+(\alpha^2+\beta^2)
\begin{pmatrix}\gamma&\delta_1\cr-\delta_2&\gamma\end{pmatrix}^2\otimes A^2\right)
\begin{pmatrix} x_1\cr x_2\end{pmatrix}=
\begin{pmatrix} d_1\cr d_2\end{pmatrix}
$$
it can be obtained as
\begin{align*}
\left(I_n+2a_1A+(a_1^2+b_1^2)A^2\right)\left(I_n+2a_2A+(a_2^2+b_2^2)A^2\right)
x_1 & =\widehat{d_1}\\
\left(I_n+2a_1A+(a_1^2+b_1^2)A^2\right)\left(I_n+2a_2A+(a_2^2+b_2^2)A^2\right)
x_2 & =\widehat{d_2}
\end{align*}
where
$$
\begin{pmatrix} \widehat{d_1}\cr\widehat{d_2}\end{pmatrix} =
\left(I_2\otimes I_n +
2\alpha\begin{pmatrix}\gamma&-\delta_1\cr\delta_2&\gamma\end{pmatrix}\otimes A
+(\alpha^2+\beta^2)
\begin{pmatrix}\gamma&-\delta_1\cr\delta_2&\gamma\end{pmatrix}^2\otimes A^2\right)
\begin{pmatrix} d_1\cr d_2\end{pmatrix}
$$
and
\begin{align*}
a_1 & = \alpha\gamma - \beta\delta\\
b_1 & = \alpha\delta + \gamma\beta\\
a_2 & = \alpha\gamma + \beta\delta\\
b_2 & = \alpha\delta - \gamma\beta\\
\delta & = \sqrt{\delta_1\delta_2}
\end{align*}

\end{lemma}

\begin{proof}
The matrix can be written as
$$\left(I_2\otimes I_n+(\alpha+\mathrm i\beta)
\begin{pmatrix}\gamma&\delta_1\cr-\delta_2&\gamma\end{pmatrix}\otimes A\right)
\left(I_2\otimes I_n +(\alpha-\mathrm i\beta)
\begin{pmatrix}\gamma&\delta_1\cr-\delta_2&\gamma\end{pmatrix}\otimes A\right).
$$
Note that the both matrices are regular since their product is
regular. For the same reason as in the previous proof, the following
matrix is also regular
$$\left(I_2\otimes I_n+(\alpha+\mathrm i\beta)
\begin{pmatrix}\gamma&-\delta_1\cr\delta_2&\gamma\end{pmatrix}\otimes A\right)
\left(I_2\otimes I_n +(\alpha-\mathrm i\beta)
\begin{pmatrix}\gamma&-\delta_1\cr\delta_2&\gamma\end{pmatrix}\otimes A\right),
$$
and we may multiply the equation by this matrix obtaining
$\widehat{d_1}$ and $\widehat{d_2}$. Note that the four matrices
commute, that is why we can write the whole product as
\begin{align*}
\left(I_2\otimes I_n + (\alpha+\mathrm i\beta)
\begin{pmatrix}\gamma&\delta_1\cr-\delta_2&\gamma\end{pmatrix}\otimes A\right)\cdot
\left(I_2\otimes I_n + (\alpha+\mathrm i\beta)
\begin{pmatrix}\gamma&-\delta_1\cr\delta_2&\gamma\end{pmatrix}\otimes A\right)&\cdot\\
\left(I_2\otimes I_n + (\alpha-\mathrm i\beta)
\begin{pmatrix}\gamma&\delta_1\cr-\delta_2&\gamma\end{pmatrix}\otimes A\right)\cdot
\left(I_2\otimes I_n + (\alpha-\mathrm i\beta)
\begin{pmatrix}\gamma&-\delta_1\cr\delta_2&\gamma\end{pmatrix}\otimes A\right)&=\\
\left(I_2\otimes I_n + 2(\alpha + \mathrm i\beta)
\begin{pmatrix}\gamma&0\cr 0&\gamma\end{pmatrix}\otimes A +
(\alpha + \mathrm i\beta)^2
\begin{pmatrix}\gamma^2+\delta^2&0\cr 0&\gamma^2+\delta^2\end{pmatrix}\otimes A^2
\right)&\cdot\\
\left(I_2\otimes I_n + 2(\alpha - \mathrm i\beta)
\begin{pmatrix}\gamma&0\cr 0&\gamma\end{pmatrix}\otimes A +
(\alpha - \mathrm i\beta)^2
\begin{pmatrix}\gamma^2+\delta^2&0\cr 0&\gamma^2+\delta^2\end{pmatrix}\otimes A^2
\right)&
\end{align*}

The product is a diagonal consiting of two $n\times n$ blocks, which are the
same. The block can be rewritten as product:
\begin{align*}
(I_n+(\alpha+\mathrm i\beta)(\gamma+\mathrm i\delta)A)\cdot
(I_n+(\alpha+\mathrm i\beta)(\gamma-\mathrm i\delta)A)&\cdot\\
(I_n+(\alpha-\mathrm i\beta)(\gamma+\mathrm i\delta)A)\cdot
(I_n+(\alpha-\mathrm i\beta)(\gamma-\mathrm i\delta)A)&
\end{align*}
and after reordering
\begin{align*}
(I_n+(\alpha+\mathrm i\beta)(\gamma+\mathrm i\delta)A)\cdot
(I_n+(\alpha-\mathrm i\beta)(\gamma-\mathrm i\delta)A)&\cdot\\
(I_n+(\alpha+\mathrm i\beta)(\gamma-\mathrm i\delta)A)\cdot
(I_n+(\alpha-\mathrm i\beta)(\gamma+\mathrm i\delta)A)&=\\
(I_n+2(\alpha\gamma-\beta\delta)A+
(\alpha^2+\beta^2)(\gamma^2+\delta^2)A^2)&\cdot\\
(I_n+2(\alpha\gamma+\beta\delta)A+
(\alpha^2+\beta^2)(\gamma^2+\delta^2)A^2)&
\end{align*}
Now it suffices to compare $a_1=\alpha\gamma-\beta\delta$ and verify
that 
\begin{align*}
b_1^2 & = (\alpha^2+\beta^2)(\gamma^2+\delta^2)-a_1^2 =\cr
      & = \alpha^2\gamma^2+\beta^2\gamma^2+\alpha^2\beta^2+\beta^2\delta^2-
          (\alpha\gamma)^2+2\alpha\beta\gamma\delta-(\beta\delta)^2=\cr
      & = (\beta\gamma)^2 + (\alpha\beta)^2 + 2\alpha\beta\gamma\delta=\cr
      & = (\beta\gamma +\alpha\beta)^2
\end{align*}
For $b_2$ it is done in the same way.
\end{proof}

\section{The Algorithm}

Below we define three functions of which
$\vec(Y)=\solvi(1,\vec(\widehat{D}),i)$ provides the solution $Y$. $X$
is then obtained as $X=U^TY\left(\otimesi V\right)$.

\subsection{Synopsis}

$F^T$ is $m\times m$ lower quasi-triangular matrix. Let $m_r$ be a
number of real eigenvalues, $m_c$ number of complex pairs. Thus
$m=m_r+2m_c$. Let $F_j$ denote
$j$-th diagonal block of $F^T$ ($1\times 1$ or $2\times 2$ matrix) for
$j=1,\ldots, m_r+m_c$. For a fixed $j$, let $\bar j$ denote index of the
first column of $F_j$ in $F^T$. Whenever we write something like
$(I_{m^i}\otimes I_n+r\iF\otimes K)x = d$, $x$ and $d$ denote column
vectors of appropriate dimensions, and $x_{\bar j}$ is $\bar j$-th
partition of $x$, and $x_j=(x_{\bar j}\quad x_{\bar j+1})^T$ if $j$-th
eigenvalue is complex, and $x_j=x_{\bar j}$ if $j$-th eigenvalue is real.

\subsection{Function $\solvi$}

The function $x=\solvi(r,d,i)$ solves equation
$$\left(I_{m^i}\otimes I_n+r\iF\otimes K\right)x = d.$$
The function proceedes as follows:

If $i=0$, the equation is solved directly, $K$ is upper
quasi-triangular matrix, so this is easy.

If $i>0$, then we go through diagonal blocks $F_j$ for
$j=1,\ldots, m_r+m_c$ and perform:
\begin{enumerate}
\item if $F_j=(f_{\bar j\bar j}) = (f)$, then we return
$x_j=\solvi(rf,d_{\bar j}, i-1)$. Then precalculate $y=d_j-x_j$, and
eliminate guys below $F_j$. This is, for each $k=\bar j+1,\ldots, m$, we
put
$$d_k = d_k-rf_{\bar jk}\left(\imF\otimes K\right)x_{\bar j}=
d_k - \frac{f_{\bar jk}}{f}y$$

\item if $F_j=\begin{pmatrix}\alpha&\beta_1\cr -\beta_2&\alpha\end{pmatrix}$,
we return $x_j=\solvii(r\alpha, r\beta_1, r\beta_2, d_j, i-1)$. Then
we precalculate 
$$y=\left(\begin{pmatrix}\alpha&-\beta_1\cr \beta_2&\alpha\end{pmatrix}
\otimes I_{m^{i-1}}\otimes I_n\right)
\begin{pmatrix} d_{\bar j} - x_{\bar j}\cr
         d_{\bar j+1} - x_{\bar j+1}
\end{pmatrix}$$
and eliminate guys below $F_j$. This is, for each $k=\bar j+2,\ldots, n$
we put
\begin{align*}
d_k &= d_k - r(f_{{\bar j}k}\quad f_{{\bar j}+1 k})
               \otimes\left(\imF\otimes K\right)x_j\\
    &= d_k - \frac{1}{\alpha^2+\beta_1\beta_2}
              \left((f_{{\bar j}k}\quad f_{{\bar j}+1 k})
               \otimes I_{m^{i-1}}\otimes I_n\right)y
\end{align*}

\end{enumerate}

\subsection{Function $\solvii$}

The function $x=\solvii(\alpha, \beta_1, \beta_2, d, i)$ solves
equation
$$
\left(I_2\otimes I_{m^i}\otimes I_n +
\begin{pmatrix}\alpha&\beta_1\cr-\beta_2&\alpha\end{pmatrix}
\otimes\iF\otimes K\right)x=d
$$

According to Lemma \ref{sylv:first-lemma} the function returns
$$
x=\begin{pmatrix}\solviip(\alpha,\beta_1\beta_2,\widehat{d_1},i)\cr
          \solviip(\alpha,\beta_1\beta_2,\widehat{d_2},i)\end{pmatrix},
$$
where $\widehat{d_1}$, and $\widehat{d_2}$ are partitions of
$\widehat{d}$ from the lemma.

\subsection{Function $\solviip$}

The function $x=\solviip(\alpha,\beta^2,d,i)$ solves equation
$$
\left(I_{m^i}\otimes I_n + 2\alpha\iF\otimes K +
(\alpha^2+\beta^2)\iF^2\otimes K^2\right)x = d
$$ 
The function proceedes as follows:

If $i=0$, the matrix $I_n+2\alpha K+(\alpha^2+\beta^2)K^2$ is
calculated and the solution is obtained directly.

Now note that diagonal blocks of $F^{2T}$ are of the form $F_j^2$,
since if the $F^T$ is block partitioned according to diagonal blocks,
then it is lower triangular.

If $i>0$, then we go through diagonal blocks $F_j$ for $j=1,\ldots, m_r+m_c$ 
and perform:
\begin{enumerate}
\item if $F_j=(f_{\bar j\bar j})=(f)$ then $j$-th diagonal block of 
$$
I_{m^i}\otimes I_n + 2\alpha\iF\otimes K +
(\alpha^2+\beta^2)\iF^2\otimes K^2
$$
takes the form
$$
I_{m^{i-1}}\otimes I_n +2\alpha f\imF\otimes K +
(\alpha^2+\beta^2)f^2\imF^2\otimes K^2
$$
and we can put $x_j = \solviip(f\alpha,f^2\beta^2,d_j,i-1)$.

Then we need to eliminate guys below $F_j$. Note that $|f^2|<|f|$,
therefore we precalculate $y_2=(\alpha^2+\beta^2)f^2(\imF^2\otimes K^2)x_j$,
and then precalculate
$$y_1 = 2\alpha f(\imF\otimes K)x_j = d_j-x_j-y_2.$$
Let $g_{pq}$ denote element of $F^{2T}$ at position $(q,p)$. 
The elimination is done by going through $k=\bar j+1,\ldots, m$ and
putting
\begin{align*}
d_k &= d_k - \left(2\alpha f_{\bar j k}\imF\otimes K +
(\alpha^2+\beta^2)g_{\bar j k}\imF^2\otimes K^2\right)x_j\\
    &= d_k - \frac{f_{\bar j k}}{f}y_1 -
       \frac{g_{\bar j k}}{f^2}y_2
\end{align*}

\item if $F_j=\begin{pmatrix}\gamma&\delta_1\cr -\delta_2&\gamma\end{pmatrix}$,
then $j$-th diagonal block of
$$
I_{m^i}\otimes I_n + 2\alpha\iF\otimes K +
(\alpha^2+\beta^2)\iF^2\otimes K^2
$$
takes the form
$$
I_{m^{i-1}}\otimes I_n +2\alpha
\begin{pmatrix}\gamma&\delta_1\cr -\delta_2&\gamma\end{pmatrix}\imF\otimes K
+(\alpha^2+\beta^2)
\begin{pmatrix}\gamma&\delta_1\cr -\delta_2&\gamma\end{pmatrix}^2\imF^2\otimes K^2
$$
According to Lemma \ref{sylv:second-lemma}, we need to calculate
$\widehat{d}_{\bar j}$, and $\widehat{d}_{\bar j+1}$, and $a_1$,
$b_1$, $a_2$, $b_2$. Then we obtain
\begin{align*}
x_{\bar j}&=
 \solviip(a_1,b_1^2,\solviip(a_2,b_2^2,\widehat{d}_{\bar j},i-1),i-1)\\
x_{\bar j+1}&=
 \solviip(a_1,b_1^2,\solviip(a_2,b_2^2,\widehat{d}_{\bar j+1},i-1),i-1)
\end{align*}

Now we need to eliminate guys below $F_j$. Since $\Vert F^2_j\Vert <
\Vert F_j\Vert$, we precalculate 

\begin{align*}
y_2&=(\alpha^2+\beta^2)(\gamma^2+\delta^2)
\left(I_2\otimes\imF^2\otimes K^2\right)x_j\\
y_1&=2\alpha(\gamma^2+\delta^2)\left(I_2\otimes\imF\otimes
K\right)x_j\\
   &=(\gamma^2+\delta^2)
      \left(F_j^{-1}
         \otimes I_{m^{i-1}n}\right)
      \left(d_j-x_j-\frac{1}{(\gamma^2+\delta^2)}
        \left(
         F_j^2
         \otimes I_{m^{i-1}n}\right)y_2\right)\\
   &=\left(\begin{pmatrix}\gamma&-\delta_1\cr\delta_2&\gamma\end{pmatrix}
           \otimes I_{m^{i-1}n}\right)(d_j-x_j)
     -\left(\begin{pmatrix}\gamma&\delta_1\cr-\delta_2&\gamma\end{pmatrix}
           \otimes I_{m^{i-1}n}\right)y_2
\end{align*}

Then we go through all $k=\bar j+2,\ldots, m$. For clearer formulas, let
$\mathbf f_k$ denote pair of $F^T$ elements in $k$-th line below $F_j$,
this is $\mathbf f_k=(f_{\bar jk}\quad f_{\bar j+1k})$. And let $\mathbf g_k$
denote the same for $F^{2T}$, this is $\mathbf g_k=(g_{\bar jk}\quad
g_{\bar j+1k})$. For each $k$ we put

\begin{align*}
d_k &= d_k - \left(2\alpha\mathbf f_k\otimes
                   \imF\otimes K +
                   (\alpha^2+\beta^2)\mathbf g_k\otimes
                   \imF^2\otimes K^2\right)x_j\\
    &= d_k - \frac{1}{\gamma^2+\delta^2}
             \left(\mathbf f_k\otimes
                   I_{m^{i-1}n}\right)y_1
           - \frac{1}{\gamma^2+\delta^2}
             \left(\mathbf g_k\otimes
                   I_{m^{i-1}n}\right)y_2
\end{align*}

\end{enumerate}

\section{Final Notes}

\subsection{Numerical Issues of $A^{-1}B$}

We began the solution of the Sylvester equation with multiplication by
$A^{-1}$.  This can introduce numerical errors, and we need more
numerically stable supplement. Its aim is to make $A$ and $B$
commutative, this is we need to find a regular matrix $P$, such that
$(PA)(PB)=(PB)(PA)$. Recall that this is neccessary in solution of
$$
(I_2\otimes I_{m^i}\otimes (PA)+(D+C)\otimes\,\iF\otimes (PB))x=d,
$$
since this equation is
multiplied by $I_2\otimes I_{m^i}\otimes (PA)+(D-C)\otimes\,\iF\otimes PB$,
and the diagonal
result 
$$
I_2\otimes I_{m^i}\otimes (PAPA) + 2D\otimes\,\iF\otimes (PAPB) +
(D^2-C^2)\otimes\,\iF^2\otimes (PBPB)
$$
is obtained only if
$(PA)(PB)=(PB)(PA)$.

Finding regular solution of $(PA)(PB)=(PB)(PA)$ is equivalent to
finding regular solution of $APB-BPA=0$. Numerical error of the former
equation is $P$-times greater than the numerical error of the latter
equation. And the numerical error of the latter equation also grows
with the size of $P$. On the other hand, truncation error in $P$
multiplication decreases with growing the size of $P$. By intuition,
stability analysis will show that the best choice is some orthonormal
$P$.

Obviously, since $A$ is regular, the equation $APB-BPA=0$ has solution
of the form $P=\alpha A^{-1}$, where $\alpha\neq 0$. There is a vector
space of all solutions $P$ (including singular ones). In precise
arithmetics, its dimension is $\sum n^2_i$, where $n_i$ is number of
repetitions of $i$-th eigenvalue of $A^{-1}B$ which is similar to
$BA^{-1}$. In floating point arithmetics, without any further
knowledge about $A$, and $B$, we are only sure about dimension $n$
which is implied by similarity of $A^{-1}B$ and $BA^{-1}$. Now we try
to find the base of the vector space of solutions.

Let $L$ denote the following linear operator:
$$L(X)=(AXB-BXA)^T.$$

Let $\vec(X)$ denote a vector made by stacking all the
columns of $X$. Let $T_n$ denote $n^2\times n^2$ matrix representing
operator $\vec(X)\mapsto \vec(X^T)$. And
finally let $M$ denote $n^2\times n^2$ matrix represening the operator
$L$. It is not difficult to verify that:
$$M=T_n(B^T\otimes A - A^T\otimes B)$$
Now we show that $M$ is skew symmetric. Recall that $T_n(X\otimes
Y)=(Y\otimes X)T_n$, we have:
$$M^T=(B^T\otimes A - A^T\otimes B)^TT_n=(B\otimes A^T - A\otimes B^T)T_n=
T_n(A^T\otimes B - B^T\otimes A) = -M$$

We try to solve $M\vec(X) = T_n(0) = 0$. Since $M$ is
skew symmetric, there is real orthonormal matrix $Q$, such that
$M=Q\widehat{M}Q^T$, and $\widehat{M}$ is block diagonal matrix
consisting of $2\times 2$ blocks of the form
$$\left(\begin{matrix} 0&\alpha_i\cr-\alpha_i&0\end{matrix}\right),$$
and of additional zero, if $n^2$ is odd.

Now we solve equation $\widehat{M}y=0$, where $y=Q^T\vec(X)$. Now
there are $n$ zero rows in $\widehat{M}$ coming from similarity of
$A^{-1}B$ and $BA^{-1}$ (structural zeros). Note that the additional
zero for odd $n^2$ is already included in that number, since for odd
$n^2$ is $n^2-n$ even. Besides those, there are also zeros (esp. in
floating point arithmetics), coming from repetitive (or close)
eigenvalues of $A^{-1}B$. If we are able to select the rows with the
structural zeros, a solution is obtained by picking arbitrary numbers
for the same positions in $y$, and put $\vec(X)=Qy$.

The following questions need to be answered:
\begin{enumerate}
\item How to recognize the structural rows?
\item Is $A^{-1}$ generated by a $y$, which has non-zero elements only
on structural rows? Note that $A$ can have repetitive eigenvalues. The
positive answer to the question implies that in each $n$-partition of
$y$ there is exactly one structural row.
\item And very difficult one: How to pick $y$ so that $X$ would be
regular, or even close to orthonormal (pure orthonormality
overdeterminates $y$)?
\end{enumerate}

\subsection{Making Zeros in $F$}

It is clear that the numerical complexity of the proposed algorithm
strongly depends on a number of non-zero elements in the Schur factor
$F$. If we were able to find $P$, such that $PFP^{-1}$ has
substantially less zeros than $F$, then the computation would be
substantially faster. However, it seems that we have to pay price for
any additional zero in terms of less numerical stability of $PFP^{-1}$
multiplication. Consider $P$, and $F$ in form
$$P=\begin{pmatrix} I &X\cr 0&I\end{pmatrix},
\qquad F=\begin{pmatrix} A&C\cr 0&B\end{pmatrix}$$
we obtain
$$PFP^{-1}=\begin{pmatrix} A& C + XB - AX\cr 0&B \end{pmatrix}$$

Thus, we need to solve $C = AX - XB$. Its clear that numerical
stability of operator $Y\mapsto PYP^{-1}$ and its inverse $Y\mapsto
P^{-1}YP$ is worse with growing norm $\Vert X\Vert$. The norm can be
as large as $\Vert F\Vert/\delta$, where $\delta$ is a distance of
eigenspectra of $A$ and $B$. Also, a numerical error of the solution is
proportional to $\Vert C\Vert/\delta$.

Although, these difficulties cannot be overcome completely, we may
introduce an algorithm, which works on $F$ with ordered eigenvalues on
diagonal, and seeks such partitioning to maximize $\delta$ and
minimize $C$. If the partitioning is found, the algorithm finds $P$
and then is run for $A$ and $B$ blocks. It stops when further
partitioning is not possible without breaking some user given limit
for numerical errors. We have to keep in mind that the numerical
errors are accumulated in product of all $P$'s of every step.

\subsection{Exploiting constant rows in $F$}

If some of $F$'s rows consists of the same numbers, or a number of
distict values within a row is small, then this structure can be
easily exploited in the algorithm. Recall, that in both functions
$\solvi$, and $\solviip$, we eliminate guys below diagonal element (or
block) (of $F^T$), by multiplying solution of the diagonal and
cancelling it from right side. If the elements below the diagonal
block are the same, we save one vector multiplication. Note that in
$\solviip$ we still need to multiply by elements below diagonal of the
matrix $F^{T2}$, which obviously has not the property. However, the
heaviest elimination is done at the very top level, in the first call
to $\solvi$.

Another way of exploitation the property is to proceed all
calculations in complex numbers. In that case, only $\solvi$ is run.

How the structure can be introduced into the matrix? Following the
same notation as in previous section, we solve $C = AX - XB$ in order
to obtain zeros at place of $C$. If it is not possible, we may relax
the equation by solving $C - R = AX - XB$, where $R$ is suitable
matrix with constant rows. The matrix $R$ minimizes $\Vert C-R\Vert$
in order to minimize $\Vert X\Vert$ if $A$, and $B$ are given. Now, in
the next step we need to introduce zeros (or constant rows) to matrix
$A$, so we seek for regular matrix $P$, doing the
job. If found, the product looks like:
$$\begin{pmatrix} P&0\cr 0&I\end{pmatrix}\begin{pmatrix} A&R\cr 0&B\end{pmatrix}
\begin{pmatrix} P^{-1}&0\cr 0&I\end{pmatrix}=
\begin{pmatrix} PAP^{-1}&PR\cr 0&B\end{pmatrix}$$
Now note, that matrix $PR$ has also constant rows. Thus,
preconditioning of the matrix in upper left corner doesn't affect the
property. However, a preconditioning of the matrix in lower right
corner breaks the property, since we would obtain $RP^{-1}$.

\end{document}