doc/tex/complex_matrices.tex

\chapter{Complex matrices}
\label{chap:complex}

\section{Introduction}
\label{sec:cmplx-intro}

Native support for complex matrices was added to gretl in version
2019d. Not all of hansl's matrix functions accept complex input, but
we have enabled a sizable subset of these functions which should
suffice for most econometric purposes.

Complex numbers are not used in most areas of econometrics, but there
are a few notable exceptions: among these, complex numbers allow for
an elegant treatment of univariate spectral analysis of time series,
and become indispensable if you consider multivariate spectral
analysis---see for example \cite{shumwaystoffer2017}. A more recent
example is the numerical solution of linear models with rational
expectations, which are widely used in modern macroeconomics, for
which the complex Schur factorization has become the tool of choice
\citep{klein2000}.

A first point to note is that complex values are treated as a special
case of the hansl \texttt{matrix} type; there's no \texttt{complex}
type as such. Complex scalars fall under the \texttt{matrix} type as
$1 \times 1$ matrices; the hansl \texttt{scalar} type is only for real
values (as is the \texttt{series} type). A $1 \times 1$ complex matrix
should do any work you might require of a complex scalar.

Before we proceed to the details of complex matrices in gretl, here's
a brief reminder of the revelant concepts and notation. Complex
numbers are pairs of the form $a + b\,i$ where $a$ and $b$ are real
numbers and $i$ is defined as the square root of $-1$: $a$ is the real
part and $b$ the imaginary part. One can specify a complex number
either via $a$ and $b$ or in ``polar'' form. The latter pertains to
the complex plane, which has the real component on the horizontal axis
and the imaginary component on the vertical. The polar representation
of a complex number is composed of the length $r$ of the ray from the
origin to the point in question and the angle $\theta$ subtended
between the positive real axis and this ray, measured
counter-clockwise in radians. In polar form the complex number
$z = a + b\,i$ can be written as
\[
  z = |z|\,(\cos \theta + i\,\sin \theta) = |z|\,e^{i\theta}
\]
where $|z| = r = \sqrt{a^2 + b^2}$ and $\theta = \tan^{-1}(b/a)$. The
quantity $|z|$ is known as the modulus of $z$, and $\theta$ as its
complex ``argument'' (or sometimes ``phase''). The notation $\bar{z}$
is used for the complex conjugate of $z$: if $z = a + b\,i$, then
$\bar{z} = a - b\,i$.


\section{Creating a complex matrix}
\label{sec:cmplx-create}

The unique explicit constructor for complex matrices is the
\cmd{complex()} function. This takes two arguments, giving the real
and imaginary parts respectively, and sticks them together, as in
\begin{code}
C = complex(A, B)
\end{code}
Four cases are supported, as follows.
\begin{itemize}
\item \texttt{A} and \texttt{B} are both $m \times n$ real matrices
  Then \texttt{C} is an $m \times n$ complex matrix such that
  $c_{kj} = a_{kj} + b_{kj}\,i$.
\item \texttt{A} and \texttt{B} are both scalars: \texttt{C} is a
  $1 \times 1$ complex matrix such that $c = a + b\,i$.
\item \texttt{A} is an $m \times n$ real matrix and \texttt{B} is a
  scalar: \texttt{C} is an $m \times n$ matrix such that
  $c_{kj} = a_{kj} + b\,i$.
\item \texttt{A} is a scalar and \texttt{B} is an $m \times n$ real
  matrix: \texttt{C} is an $m \times n$ matrix such that
  $c_{kj} = a + b_{kj}\,i$.
\end{itemize}

In addition, complex matrices may naturally arise as the result of
certain computations.

With both real and complex matrices in circulation, one may wish to
determine whether a particular matrix is complex. The function
\cmd{iscomplex()} can tell you. Passed an identifier, it returns 1
if it names a complex matrix, 0 if it names a real matrix, or
\texttt{NA} otherwise.

Note, however, that the \cmd{iscomplex()} function only tells you if a
certain matrix is endowed with an imaginary part, which may be
zero. The following code snippet should clarify the point:
\begin{code}
matrix z = complex(1,0)
scalar a = iscomplex(z)
scalar b = z[imag] == 0
printf "a = %g, b = %g\n", a, b
\end{code}
The code above gives
\begin{code}
a = 1, b = 1
\end{code}
The test \texttt{a} is non-zero (or ``true'') because the matrix
\texttt{z} is defined as complex, but \texttt{b}, which tests for an
all-zero imaginary part of \texttt{z}, is also true. In mathematical
terms, then, \texttt{z} is effectively a real matrix.

\section{Indexation}

Indexation of complex matrices works as with real matrices, on the
understanding that each element of a complex matrix is a complex
pair. So for example \texttt{C[i,j]} gets you the complex pair at row
\texttt{i}, column \texttt{j} of \texttt{C}, in the form of a
$1 \times 1$ complex matrix.

If you wish to access just the real or imaginary part of a given
element, or range of elements, you can use the functions \cmd{Re()}
or \cmd{Im()}, as in
\begin{code}
scalar rij = Re(C[i,j])
\end{code}
which gets you the real part of $c_{ij}$.

In addition the dummy selectors \cmd{real} and \cmd{imag} can be
used to assign to just the real or imaginary component of a complex
matrix. Here are two examples:
\begin{code}
# replace the real part of C with random normals
C[real] = mnormal(rows(C), cols(C))

# set the imaginary part of C to all zeros
C[imag] = 0
\end{code}
The replacement must be either a real matrix of the same dimensions as
the target, or a scalar.

Further, the \cmd{real} and \cmd{imag} selectors may be combined
with regular selectors to access specific portions of a complex matrix
for either reading or writing. Examples:
\begin{code}
# retrieve the real part of a submatrix of C
matrix R = C[1:2,1:2][real]

# set the imaginary part of C[3,3] to y
C[3,3][imag] = y
\end{code}

\section{Operators}
\label{sec:cmplx-ops}

Most of the operators available for working with real matrices are
also available for complex ones; this includes the ``dot-operators''
which work element-wise or by ``broadcasting'' vectors. Moreover,
mixed operands are accepted, as in \texttt{D = C + A} where \texttt{C}
is complex and \texttt{A} real; the result, \texttt{D}, will be
complex. In such cases the real operand is treated as a complex matrix
with an all-zero imaginary part.

The operators \textit{not} defined for complex values are:
\begin{itemize}
\item Those that include the inequality tests ``\verb+>+'' or
  ``\verb+<+'', since complex values as such cannot be compared as
  greater or lesser (though they can be compared as equal or not
  equal).
\item The (real) modulus operator (percent sign), as in \texttt{x \%
    y} which gives the remainder on division of \texttt{x} by
  \texttt{y}.
\end{itemize}

As for real matrices, the transposition operator ``\cmd{'}'' is
available in both unary form, as in \texttt{B = A'}, and binary form,
as in \texttt{C = A'B} (transpose-multiply). But note that for complex
\texttt{A} this means the conjugate transpose, $A^\mathrm{H}$. If you
need the non-conjugated transpose you can use \cmd{transp()}.

You may wish to note: although none of gretl's explicit regression
functions (or commands) accept complex input you can calculate
parameter estimates for a least-squares regression of complex $Y$
($T \times 1$) on complex $X$ ($T \times k$) via \verb|B = X \ Y|.

\section{Functions}
\label{sec:cmplx-funcs}

To give an idea of what works, and what doesn't, for complex
matrices, we'll walk through the hansl function-space using the
categories employed in gretl's online ``Function reference'' (under the
\textsf{Help} menu in the GUI program).

\subsection{Linear algebra}

The functions that accept complex arguments are: \cmd{cholesky},
\cmd{det}, \cmd{ldet}, \cmd{eigensym} (for Hermitian matrices),
\cmd{ffti}, \cmd{inv}, \cmd{ginv}, \cmd{hdprod}, \cmd{mexp},
\cmd{mlog}, \cmd{qrdecomp}, \cmd{rank}, \cmd{svd}, \cmd{tr}, and
\cmd{transp}. Note, however, that \cmd{mexp} and \cmd{mlog} require
that the input matrix be diagonalizable, and \cmd{cholesky} requires a
positive definite Hermitian matrix.

In addition the new functions \cmd{eigen} and \cmd{fft2} are
complex-supporting versions of \cmd{eigengen} and \cmd{fft},
respectively (see section~\ref{sec:cmplx-compat} for details). And
there are the complex-only functions \cmd{ctrans}, which gives the
conjugate transpose,\footnote{The \cmd{transp} function gives the
  straight (non-conjugated) transpose of a complex matrix.} and
\cmd{schur} for the Schur factorization.

\subsection{Matrix building}

Given what was said in section~\ref{sec:cmplx-create} above, several
of the functions in this category should be thought of as applying to
the real or imaginary part of a complex matrix (for example,
\cmd{ones} and \cmd{mnormal}), and are of course usable in that
way.  However, some of these functions can be applied to complex
matrices as such, namely, \cmd{diag}, \cmd{diagcat},
\cmd{lower}, \cmd{upper}, \cmd{vec}, \cmd{vech} and
\cmd{unvech}.

Please note: when \cmd{unvech} is applied to a suitable real
vector it produces a symmetric matrix, but when applied to a complex
vector it produces a Hermitian matrix.

The only functions \textit{not} available for complex matrices are
\cmd{cnameset} and \cmd{rnameset}. That is, you cannot name the
columns or rows of such matrices (although this restriction could
probably be lifted without great difficulty).

\subsection{Matrix shaping}

The functions that accept complex input are: \cmd{cols},
\cmd{rows}, \cmd{mreverse}, \cmd{mshape}, \cmd{selifc},
\cmd{selifr} and \cmd{trimr}.

The functions \cmd{msortby}, \cmd{sort} and \cmd{dsort} are
excluded for the reason mentioned in section~\ref{sec:cmplx-ops}.

\subsection{Statistical}

Supported for complex input: \cmd{meanc}, \cmd{meanr},
\cmd{sumc}, \cmd{sumr}, \cmd{prodc} and \cmd{prodr}. And
that's all.

\subsection{Mathematical}

In the matrix context, these are functions that are applied element by
element. For complex input the following are supported: \cmd{log},
\cmd{exp} and \cmd{sqrt}, plus all of the trigonometric
functions with the exception of \cmd{atan2}.

In addition there are the complex-only functions \cmd{cmod}
(complex modulus, also accessible via \cmd{abs}), \cmd{carg}
(complex argument), \cmd{conj} (complex conjugate), \cmd{Re}
(real part) and \cmd{Im} (imaginary part). Note that
$\mbox{carg}(z) = \mbox{atan2}(y,x)$ for $z=x +
y\,i$. Listing~\ref{cmplx-modes} illustrates usage of \cmd{cmod}
and \cmd{carg}.

\begin{script}[htbp]
  \scriptcaption{Variant representations of complex numbers. We picked 8
    points on the unit circle in the complex plane, so their modulus
    is constant and equal to 1. The \texttt{Polar} matrix below shows
    that the complex argument is expressed in radians; multiplying by
    180/$\pi$ gives degrees. The \texttt{chk} matrix verifies that
    we can retrieve the orginal representation of the complex values
    from the polar form in either of the two ways mentioned at the
    start of the chapter: $z = |z|\,(\cos \theta + i\,\sin \theta)$ or
    $z = |z|\,e^{i\theta}$.}
  \label{cmplx-modes}
\begin{scodebit}
# complex values in a + b*i form
scalar rp5 = sqrt(0.5)
matrix A = {1, rp5, 0, -rp5, -1, -rp5,  0,  rp5}'
matrix B = {0, rp5, 1,  rp5,  0, -rp5, -1, -rp5}'
matrix Z = complex(A, B)

# calculate modulus and argument
matrix zmod = cmod(Z)
matrix theta = carg(Z)
matrix Polar = zmod ~ theta ~ (theta * 180/$pi)
cnameset(Polar, "modulus radians degrees")
printf "%12.4f\n", Polar

# reconstitute the original Z matrix in two ways
matrix Z1 = zmod .* complex(cos(theta), sin(theta))
matrix Z2 = zmod .* exp(complex(0, theta))
matrix chk = Z ~ Z1 ~ Z2
print chk
\end{scodebit}


  Printing of \texttt{Polar} and \texttt{chk}
\begin{outbit}
     modulus     radians     degrees
      1.0000      0.0000      0.0000
      1.0000      0.7854     45.0000
      1.0000      1.5708     90.0000
      1.0000      2.3562    135.0000
      1.0000      3.1416    180.0000
      1.0000     -2.3562   -135.0000
      1.0000     -1.5708    -90.0000
      1.0000     -0.7854    -45.0000

 1.00000 + 0.00000i   1.00000 + 0.00000i   1.00000 + 0.00000i
 0.70711 + 0.70711i   0.70711 + 0.70711i   0.70711 + 0.70711i
 0.00000 + 1.00000i   0.00000 + 1.00000i   0.00000 + 1.00000i
-0.70711 + 0.70711i  -0.70711 + 0.70711i  -0.70711 + 0.70711i
-1.00000 + 0.00000i  -1.00000 + 0.00000i  -1.00000 + 0.00000i
-0.70711 - 0.70711i  -0.70711 - 0.70711i  -0.70711 - 0.70711i
 0.00000 - 1.00000i   0.00000 - 1.00000i   0.00000 - 1.00000i
 0.70711 - 0.70711i   0.70711 - 0.70711i   0.70711 - 0.70711i
\end{outbit}
\end{script}

\subsection{Transformations}

In this category only two functions can be applied to complex
matrices, namely \cmd{cum} and \cmd{diff}.

\section{File input/output}

Complex matrices should be stored and retrieved correctly in the
XML serialization used for gretl session files (\texttt{*.gretl}).

The functions \cmd{mwrite} and \cmd{mread} work in two modes:
binary mode if the filename ends with ``\texttt{.bin}'' and text mode
otherwise. Both modes handle complex matrices correctly if both the
writing and the reading are to be done by gretl, but for exchange of
data with ``foreign'' programs text mode will \textit{not} work for
complex matrices as a whole. The options are:
\begin{itemize}
\item In text mode, use \cmd{mwrite} and \cmd{mread} on the two
  parts of a complex matrix separately, and reassemble the matrix in
  the target program.
\item Use binary mode (on the whole matrix), if this is supported for
  the given foreign program.
\end{itemize}

At present binary mode transfer of complex matrices is supported for
\textsf{octave}, \textsf{python} and \textsf{julia}.
Listing~\ref{cmplx-io} shows some examples: we export a complex matrix
to each of these programs in turn; calculate its inverse in the
foreign program; then verify that the result as imported back into
gretl is the same as that calculated in gretl.

\begin{script}[htbp]
  \scriptcaption{Exporting and importing complex matrices}
  \label{cmplx-io}
\begin{scode}
set seed 34756
matrix C = complex(mnormal(3,3), mnormal(3,3))
D = inv(C)

mwrite(C, "C.bin", 1)

foreign language=octave
  C = gretl_loadmat('C.bin');
  gretl_export(inv(C), 'oct_D.bin');
end foreign

oct_D = mread("oct_D.bin", 1)
eval D - oct_D

foreign language=python
   import numpy as np
   C = gretl_loadmat('C.bin')
   gretl_export(np.linalg.inv(C), 'py_D.bin')
end foreign

py_D = mread("py_D.bin", 1)
eval D - py_D

foreign language=julia
  C = gretl_loadmat("C.bin")
  gretl_export(inv(C), "jl_D.bin")
end foreign

jl_D = mread("jl_D.bin", 1)
eval D - jl_D
\end{scode}
\end{script}

\section{Backward compatibility}
\label{sec:cmplx-compat}

Compatibility issues arise in two contexts, both related to the fact
that gretl offered some degree of support for complex matrices before
they became full citizens of the hansl polity.

\begin{enumerate}
\item The functions \cmd{fft} (fast Fourier transform for real
  input) and \cmd{eigengen} (eigenvalues and/or eigenvectors of a
  non-symmetric real matrix) returned complex matrices in what we may
  call the ``legacy'' representation. In the case of \cmd{fft} and
  the eigenvalues from \cmd{eigengen} this took the form of a
  regular gretl matrix with real values in the first (or odd-numbered)
  column(s) and imaginary parts in the second (or even-numbered)
  column(s). Since calculating with such matrices using the standard
  matrix operators would result in nonsense, we provided the tailored
  functions \cmd{cmult} and \cmd{cdiv}.

  In the case of complex eigenvectors from \cmd{eigengen}---well,
  you probably don't want to know, but if you do, consult the help text
  for \cmd{eigengen}; they were not easy for a user to handle!
\item The function packages \texttt{cmatrix} and
  \texttt{ghosts}. These were developed to support frequency-domain
  analysis via gretl in the absence of built-in complex matrix
  functionality. The \texttt{cmatrix} package was needed as an
  dependency for \texttt{ghosts} (multivariate spectral analysis).
\end{enumerate}

So what happens with these functions and function packages under the
new regime? The resolution on the two built-in functions is this:
\begin{itemize}
\item \cmd{fft} and \cmd{eigengen} continue to behave exactly as
  before. They do not accept complex input and they produce old-style
  output. In the documentation they are marked as legacy functions,
  not for use in newly written hansl code.
\item We have added new counterpart functions, \cmd{fft2} and
  \cmd{eigen}. These accept either real or complex input and they
  produce new-style complex output in both cases.
\end{itemize}

On the affected packages: \texttt{cmatrix} is no longer required, and
is not be supported any more. But an updated version of
\texttt{ghosts}, which uses gretl's native complex functionality, is
available.

We might mention that the \cmd{ffti} function (inverse Fourier
transform) is backward compatible: the new functionality is just a
superset of the old. The input must be complex, and is accepted by
\cmd{ffti} in either the legacy or the new format. The output is
real if the input is Hermitian, which in this case means that the
first (zero-frequency) row is real and the remaining rows are
conjugate symmetrical about their mid-point. That's the only case that
can arise when \cmd{ffti} is used on output from the old
\cmd{fft} (which only accepted real input).  Otherwise
(non-Hermitian complex input, which can arise only under the new
scheme) the output will be complex, and in the new format.

%%% Local Variables:
%%% mode: latex
%%% TeX-master: "gretl-guide"
%%% End: