Skip to content
Snippets Groups Projects
Verified Commit c15602c6 authored by Sébastien Villemot's avatar Sébastien Villemot
Browse files

Dynare++: convert doc of Sylvester module to LaTeX 

Previously it was plain TeX.

This allows us to remove the test for pdfetex.
parent 268d540d
No related branches found
No related tags found
No related merge requests found
Pipeline #700 failed
Stage: build
Stage: test
Show whitespace changes
Inline Side-by-side
configure.ac
+ 0
− 3
View file @ c15602c6
......@@ -87,9 +87,6 @@ AM_CONDITIONAL([HAVE_MATIO], [test "x$has_matio" = "xyes"])
AC_CHECK_PROG([MAKEINFO], [makeinfo], [makeinfo])
AC_CHECK_PROG([PDFETEX], [pdfetex], [pdfetex])
AM_CONDITIONAL([HAVE_PDFETEX], [test "x$PDFETEX" != "x"])
AC_CHECK_PROG([PDFLATEX], [pdflatex], [pdflatex])
AM_CONDITIONAL([HAVE_PDFLATEX], [test "x$PDFLATEX" != "x"])
......
dynare++/sylv/Makefile.am
+ 5
− 3
View file @ c15602c6
......@@ -2,11 +2,13 @@ SUBDIRS = cc testing
EXTRA_DIST = sylvester.tex change_log.html matlab
if HAVE_PDFETEX
if HAVE_PDFLATEX
pdf-local: sylvester.pdf
sylvester.pdf: sylvester.tex
$(PDFETEX) sylvester
$(PDFLATEX) sylvester
$(PDFLATEX) sylvester
endif
CLEANFILES = *.pdf *.log
CLEANFILES = *.pdf *.log *.aux
dynare++/sylv/sylvester.tex
+ 165
− 170
View file @ c15602c6
\input amstex
\documentstyle{amsppt}
\def\vec{\mathop{\hbox{vec}}}
\def\otimesi{\mathop{\overset{\ssize i}\to{\otimes}}}
\def\iF{\,^i\!F}
\def\imF{\,^{i-1}\!F}
\def\solvi{\bold{solv1}}
\def\solvii{\bold{solv2}}
\def\solviip{\bold{solv2p}}
\topmatter
\title Solution of Specialized Sylvester Equation\endtitle
\author Ondra Kamenik\endauthor
\email ondrej.kamenik@cnb.cz\endemail
\endtopmatter
\document
\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
......@@ -36,173 +43,168 @@ 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.
\proclaim{Lemma 1}
\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 +
\pmatrix\alpha&\beta_1\cr-\beta_2&\alpha\endpmatrix
\otimes A\right)\pmatrix x_1\cr x_2\endpmatrix =
\pmatrix d_1\cr d_2\endpmatrix,$$
\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
$$\align
\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}
\endalign$$
\end{align*}
where $\beta=\sqrt{\beta1\beta2}$, and
$$
\pmatrix \widehat{d_1}\cr\widehat{d_2}\endpmatrix =
\begin{pmatrix} \widehat{d_1}\cr\widehat{d_2}\end{pmatrix} =
\left(I_2\otimes I_n +
\pmatrix\alpha&-\beta_1\cr\beta_2&\alpha\endpmatrix
\otimes A\right)\pmatrix d_1\cr d_2\endpmatrix$$
\endproclaim
\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}
\demo{Proof} Since
\begin{proof} Since
$$
\pmatrix \alpha&\beta_1\cr-\beta_2&\alpha\endpmatrix
\pmatrix \alpha&-\beta_1\cr\beta_2&\alpha\endpmatrix
\begin{pmatrix} \alpha&\beta_1\cr-\beta_2&\alpha\end{pmatrix}
\begin{pmatrix} \alpha&-\beta_1\cr\beta_2&\alpha\end{pmatrix}
=
\pmatrix \alpha&-\beta_1\cr\beta_2&\alpha\endpmatrix
\pmatrix \alpha&\beta_1\cr-\beta_2&\alpha\endpmatrix
\begin{pmatrix} \alpha&-\beta_1\cr\beta_2&\alpha\end{pmatrix}
\begin{pmatrix} \alpha&\beta_1\cr-\beta_2&\alpha\end{pmatrix}
=
\pmatrix \alpha^2+\beta^2&0\cr 0&\alpha^2+\beta^2\endpmatrix,
\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 +
\pmatrix\alpha&-\beta_1\cr\beta_2&\alpha\endpmatrix
\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
$$\pmatrix \alpha&-\beta_1\cr\beta_2&\alpha\endpmatrix$$
$$\begin{pmatrix} \alpha&-\beta_1\cr\beta_2&\alpha\end{pmatrix}$$
collapses an eigenvalue of $A$ to $-1$ iff the matrix
$$\pmatrix \alpha&\beta_1\cr-\beta_2&\alpha\endpmatrix$$
does.\qed
\enddemo
$$\begin{pmatrix} \alpha&\beta_1\cr-\beta_2&\alpha\end{pmatrix}$$
does.
\end{proof}
\proclaim{Lemma 2}
\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\pmatrix\gamma&\delta_1\cr-\delta_2&\gamma\endpmatrix\otimes A
2\alpha\begin{pmatrix}\gamma&\delta_1\cr-\delta_2&\gamma\end{pmatrix}\otimes A
+(\alpha^2+\beta^2)
\pmatrix\gamma&\delta_1\cr-\delta_2&\gamma\endpmatrix^2\otimes A^2\right)
\pmatrix x_1\cr x_2\endpmatrix=
\pmatrix d_1\cr d_2\endpmatrix
\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
$$
\align
\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}
\endalign$$
\end{align*}
where
$$
\pmatrix \widehat{d_1}\cr\widehat{d_2}\endpmatrix =
\begin{pmatrix} \widehat{d_1}\cr\widehat{d_2}\end{pmatrix} =
\left(I_2\otimes I_n +
2\alpha\pmatrix\gamma&-\delta_1\cr\delta_2&\gamma\endpmatrix\otimes A
2\alpha\begin{pmatrix}\gamma&-\delta_1\cr\delta_2&\gamma\end{pmatrix}\otimes A
+(\alpha^2+\beta^2)
\pmatrix\gamma&-\delta_1\cr\delta_2&\gamma\endpmatrix^2\otimes A^2\right)
\pmatrix d_1\cr d_2\endpmatrix
\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
$$
\align
\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}
\endalign$$
\endproclaim
\end{align*}
\demo{Proof}
\end{lemma}
\begin{proof}
The matrix can be written as
$$\left(I_2\otimes I_n+(\alpha+\roman i\beta)
\pmatrix\gamma&\delta_1\cr-\delta_2&\gamma\endpmatrix\otimes A\right)
\left(I_2\otimes I_n +(\alpha-\roman i\beta)
\pmatrix\gamma&\delta_1\cr-\delta_2&\gamma\endpmatrix\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)
\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+\roman i\beta)
\pmatrix\gamma&-\delta_1\cr\delta_2&\gamma\endpmatrix\otimes A\right)
\left(I_2\otimes I_n +(\alpha-\roman i\beta)
\pmatrix\gamma&-\delta_1\cr\delta_2&\gamma\endpmatrix\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)
\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
$$
\align
\left(I_2\otimes I_n + (\alpha+\roman i\beta)
\pmatrix\gamma&\delta_1\cr-\delta_2&\gamma\endpmatrix\otimes A\right)\cdot
\left(I_2\otimes I_n + (\alpha+\roman i\beta)
\pmatrix\gamma&-\delta_1\cr\delta_2&\gamma\endpmatrix\otimes A\right)&\cdot\\
\left(I_2\otimes I_n + (\alpha-\roman i\beta)
\pmatrix\gamma&\delta_1\cr-\delta_2&\gamma\endpmatrix\otimes A\right)\cdot
\left(I_2\otimes I_n + (\alpha-\roman i\beta)
\pmatrix\gamma&-\delta_1\cr\delta_2&\gamma\endpmatrix\otimes A\right)&=\\
\left(I_2\otimes I_n + 2(\alpha + \roman i\beta)
\pmatrix\gamma&0\cr 0&\gamma\endpmatrix\otimes A +
(\alpha + \roman i\beta)^2
\pmatrix\gamma^2+\delta^2&0\cr 0&\gamma^2+\delta^2\endpmatrix\otimes A^2
\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 - \roman i\beta)
\pmatrix\gamma&0\cr 0&\gamma\endpmatrix\otimes A +
(\alpha - \roman i\beta)^2
\pmatrix\gamma^2+\delta^2&0\cr 0&\gamma^2+\delta^2\endpmatrix\otimes A^2
\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)&
\endalign
$$
\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:
$$
\align
(I_n+(\alpha+\roman i\beta)(\gamma+\roman i\delta)A)\cdot
(I_n+(\alpha+\roman i\beta)(\gamma-\roman i\delta)A)&\cdot\\
(I_n+(\alpha-\roman i\beta)(\gamma+\roman i\delta)A)\cdot
(I_n+(\alpha-\roman i\beta)(\gamma-\roman i\delta)A)&
\endalign
$$
\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
$$
\align
(I_n+(\alpha+\roman i\beta)(\gamma+\roman i\delta)A)\cdot
(I_n+(\alpha-\roman i\beta)(\gamma-\roman i\delta)A)&\cdot\\
(I_n+(\alpha+\roman i\beta)(\gamma-\roman i\delta)A)\cdot
(I_n+(\alpha-\roman i\beta)(\gamma+\roman i\delta)A)&=\\
\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)&
\endalign
$$
\end{align*}
Now it suffices to compare $a_1=\alpha\gamma-\beta\delta$ and verify
that
$$
\align
\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
\endalign
$$
\end{align*}
For $b_2$ it is done in the same way.
\qed
\enddemo
\end{proof}
\head The Algorithm\endhead
\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)$.
\subhead Synopsis\endsubhead
\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
......@@ -215,7 +217,7 @@ 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.
\subhead Function $\solvi$\endsubhead
\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.$$
......@@ -226,7 +228,7 @@ 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:
\roster
\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
......@@ -234,46 +236,45 @@ 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=\pmatrix\alpha&\beta_1\cr -\beta_2&\alpha\endpmatrix$,
\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(\pmatrix\alpha&-\beta_1\cr \beta_2&\alpha\endpmatrix
$$y=\left(\begin{pmatrix}\alpha&-\beta_1\cr \beta_2&\alpha\end{pmatrix}
\otimes I_{m^{i-1}}\otimes I_n\right)
\pmatrix d_{\bar j} - x_{\bar j}\cr
\begin{pmatrix} d_{\bar j} - x_{\bar j}\cr
d_{\bar j+1} - x_{\bar j+1}
\endpmatrix$$
\end{pmatrix}$$
and eliminate guys below $F_j$. This is, for each $k=\bar j+2,\ldots, n$
we put
$$
\align
\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
\endalign
$$
\endroster
\end{align*}
\end{enumerate}
\subhead Function $\solvii$\endsubhead
\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 +
\pmatrix\alpha&\beta_1\cr-\beta_2&\alpha\endpmatrix
\begin{pmatrix}\alpha&\beta_1\cr-\beta_2&\alpha\end{pmatrix}
\otimes\iF\otimes K\right)x=d
$$
According to {\bf Lemma 1} the function returns
According to Lemma \ref{sylv:first-lemma} the function returns
$$
x=\pmatrix\solviip(\alpha,\beta_1\beta_2,\widehat{d_1},i)\cr
\solviip(\alpha,\beta_1\beta_2,\widehat{d_2},i)\endpmatrix,
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.
\subhead Function $\solviip$\endsubhead
\subsection{Function $\solviip$}
The function $x=\solviip(\alpha,\beta^2,d,i)$ solves equation
$$
......@@ -291,7 +292,7 @@ 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:
\roster
\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 +
......@@ -311,16 +312,14 @@ $$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
$$
\align
\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
\endalign
$$
\end{align*}
\item if $F_j=\pmatrix\gamma&\delta_1\cr -\delta_2&\gamma\endpmatrix$,
\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 +
......@@ -329,26 +328,24 @@ $$
takes the form
$$
I_{m^{i-1}}\otimes I_n +2\alpha
\pmatrix\gamma&\delta_1\cr -\delta_2&\gamma\endpmatrix\imF\otimes K
\begin{pmatrix}\gamma&\delta_1\cr -\delta_2&\gamma\end{pmatrix}\imF\otimes K
+(\alpha^2+\beta^2)
\pmatrix\gamma&\delta_1\cr -\delta_2&\gamma\endpmatrix^2\imF^2\otimes K^2
\begin{pmatrix}\gamma&\delta_1\cr -\delta_2&\gamma\end{pmatrix}^2\imF^2\otimes K^2
$$
According to {\bf Lemma 2}, we need to calculate
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
$$
\align
\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)
\endalign
$$
\end{align*}
Now we need to eliminate guys below $F_j$. Since $\Vert F^2_j\Vert <
\Vert F_j\Vert$, we precalculate
$$
\align
\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
......@@ -360,37 +357,36 @@ K\right)x_j\\
\left(
F_j^2
\otimes I_{m^{i-1}n}\right)y_2\right)\\
&=\left(\pmatrix\gamma&-\delta_1\cr\delta_2&\gamma\endpmatrix
&=\left(\begin{pmatrix}\gamma&-\delta_1\cr\delta_2&\gamma\end{pmatrix}
\otimes I_{m^{i-1}n}\right)(d_j-x_j)
-\left(\pmatrix\gamma&\delta_1\cr-\delta_2&\gamma\endpmatrix
-\left(\begin{pmatrix}\gamma&\delta_1\cr-\delta_2&\gamma\end{pmatrix}
\otimes I_{m^{i-1}n}\right)y_2
\endalign
$$
\end{align*}
Then we go through all $k=\bar j+2,\ldots, m$. For clearer formulas, let
$\bold f_k$ denote pair of $F^T$ elements in $k$-th line below $F_j$,
this is $\bold f_k=(f_{\bar jk}\quad f_{\bar j+1k})$. And let $\bold g_k$
denote the same for $F^{2T}$, this is $\bold g_k=(g_{\bar jk}\quad
$\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
$$
\align
d_k &= d_k - \left(2\alpha\bold f_k\otimes
\begin{align*}
d_k &= d_k - \left(2\alpha\mathbf f_k\otimes
\imF\otimes K +
(\alpha^2+\beta^2)\bold g_k\otimes
(\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(\bold f_k\otimes
\left(\mathbf f_k\otimes
I_{m^{i-1}n}\right)y_1
- \frac{1}{\gamma^2+\delta^2}
\left(\bold g_k\otimes
\left(\mathbf g_k\otimes
I_{m^{i-1}n}\right)y_2
\endalign
$$
\end{align*}
\endroster
\end{enumerate}
\head Final Notes\endhead
\section{Final Notes}
\subhead Numerical Issues of $A^{-1}B$\endsubhead
\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
......@@ -448,7 +444,7 @@ 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(\matrix 0&\alpha_i\cr-\alpha_i&0\endmatrix\right),$$
$$\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
......@@ -462,7 +458,7 @@ 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:
\roster
\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
......@@ -471,9 +467,9 @@ $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$)?
\endroster
\end{enumerate}
\subhead Making Zeros in $F$\endsubhead
\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
......@@ -482,10 +478,10 @@ 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=\pmatrix I &X\cr 0&I\endpmatrix,
\qquad F=\pmatrix A&C\cr 0&B\endpmatrix$$
$$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}=\pmatrix A& C + XB - AX\cr 0&B \endpmatrix$$
$$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
......@@ -503,7 +499,7 @@ 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.
\subhead Exploiting constant rows in $F$\endsubhead
\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
......@@ -529,13 +525,12 @@ 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:
$$\pmatrix P&0\cr 0&I\endpmatrix\pmatrix A&R\cr 0&B\endpmatrix
\pmatrix P^{-1}&0\cr 0&I\endpmatrix=
\pmatrix PAP^{-1}&PR\cr 0&B\endpmatrix$$
$$\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}$.
\enddocument
\end{document}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment