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

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

Native support for complex matrices was added to gretl in version
2019d.\footnote{Prior to that release gretl offered improvised support
  for some complex functionality; see section~\ref{sec:cmplx-compat}
  for details.} 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 relevant 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 standard 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
non-zero if it names a complex matrix, 0 if it names a real matrix, or
\texttt{NA} otherwise. The non-zero return value is either 1 or 2,
with the following interpretation:
\begin{itemize}
\item 1 indicates that the matrix is ``nominally complex'' (each
  element is represented as having a real part and an imaginary part)
  but all imaginary parts are zero.
\item 2 indicates that at least one element has a non-zero imaginary
  part.
\end{itemize}
The following code snippet illustrates the point.
\begin{code}
matrix z1 = complex(1,0)
scalar a = iscomplex(z1)
matrix z2 = complex(1,1)
scalar b = iscomplex(z2)
printf "a = %d, b = %d\n", a, b
\end{code}
The code above gives
\begin{code}
a = 1, b = 2
\end{code}

\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{eigen}, \cmd{eigensym} (for Hermitian
matrices), \cmd{fft}, \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 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{ex:complex-modes} illustrates usage of \cmd{cmod}
and \cmd{carg}.

\begin{script}[htbp]
  \scriptinfo{complex-modes}{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}$.}
\begin{scode}
# 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{scode}

  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 are 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{ex:complex-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]
  \scriptinfo{complex-io}{Exporting and importing complex matrices}
\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 (in)compatibility}
\label{sec:cmplx-compat}

Prior to version 2019d gretl did not provide native support for
complex matrices. It did, however, offer an improvised representation
of such matrices for certain restricted purposes, taking the form of
an expanded regular gretl matrix with real values and imaginary parts
in odd- and even-numbered columns, respectively.  The functions
\cmd{fft}, \cmd{eigengen} and \cmd{polroots} returned matrices in this
special form, and the functions \cmd{cmult} and \cmd{cdiv} operated on
such matrices.

As of version 2022b, \cmd{fft} and \cmd{polroots} have been redefined
to work with ``proper'' complex matrices as described above. The other
affected functions are deprecated and will be removed or redefined in
a subsequent release. If you have any hansl code using the legacy
representation the following brief porting guide may be helpful.

\subsection{Porting old complex code}

\texttt{cmult} and \texttt{cdiv}: These functions performed
element-wise multiplication and division of complex column vectors in
the old two-column form. The statements
\begin{code}
# old element-wise operations
c1 = cmult(a1, b1)
d1 = cdiv(a1, b1)
\end{code}
can be updated as
\begin{code}
# new element-wise operations
c2 = a2 .* b2
d2 = a2 ./ b2
\end{code}
(where \texttt{a2} and \texttt{b2} are new-style complex vectors or
matrices). The following statements
\begin{code}
c3 = a2 * b2
d3 = a2 / b2
\end{code}
are also valid but have different effects, the first performing
standard (rather than element-wise) multiplication of matrices
(complex or real) and the second performing ``right division'',
equivalent to \texttt{a2 * inv(b2)}. Note that while the return value
from \texttt{cmult} and \texttt{cdiv} could be either a real vector or
a (two column) complex vector, the new-style operations yield a
nominally complex result if at least one of the operands is complex,
even if the result has an all-zero imaginary part.

A piece of code that appears in some contexts (such as calculation of
a periodogram) is as follows: given a complex vector, \texttt{v},
compute a vector \texttt{w} holding the squared moduli of the elements
of \texttt{v}. The old-style code to accomplish this was
\begin{code}
# legacy: v has two columns
w = sumr(v.^2)
\end{code}
and the new replacement is
\begin{code}
# current: v has a single complex column
w = abs(v).^2
\end{code}
where \texttt{abs} gives the complex modulus.

\texttt{eigengen}: Most uses of this legacy function simply retrieve
the eigenvalues of a general (that is, not symmetric) matrix, and do
not exploit the option of retrieving eigenvectors. In that context it
is straightforward to substitute a call to the new function
\texttt{eigen}. The only point to note is that \texttt{eigen} returns
a new-style complex vector; if you have need to convert this to the
legacy representation you can use the \texttt{cswitch} function, which
is documented in the \GCR. In brief, the following code gives you the
legacy equivalent of a new-style complex vector \texttt{v}:
\texttt{newvec}.
\begin{code}
if v[imag] == 0
   oldv = v[real]
else
   oldv = cswitch(v, 2)
endif
\end{code}

\texttt{polroots}: This function now returns a new-style complex
vector. As with \texttt{eigengen}, you can use \texttt{cswitch} to
convert the vector if necessary.

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