\chapter{Matrix Functions}

There are several functions\index{functions } provided for the
requirements of matrix algebra. The source code is provided so
you can expand this section with your own source code. There are
many popular sources of mathematical software in C which can be
translated to C++. Here is a list of authors I like to use for
code in C: Holub\cite{Ho:dd}, Baker\cite{Bk:ct}\cite{Bk:mt},
Kassab\cite{Ka:tc}, and Press\cite{Pr:nrc}. These sources range
from easy to hard. Most of the functions below are translations of C,
ALGOL, Pascal or FORTRAN programs I have in my library. I have also
produced some simple functions found in most matrix packages. Most of
the routines are contained in VIRTOP.CPP, and work the same way
under the in--ram and virtual memory models.


\section{The Functions}
There are five classes of functions in this package:
display control, IO, subsetting and concatenation,
elementwise operations, and math functions. Their
prototypes are contained below: 
\begin{verbatim}

typedef enum boolean { FFALSE, TTRUE } boolean;

int Setwid( int w );                      // set display width
int Setdec( int d );                      // set decimals
int Getwid( void  );                      // report display width
int Getdec( void  );                      // report decimals

////////////////////////// io

VMatrix& Reada( char *fid);               // read ascii file
VMatrix& Readb( char *fid);               // read binary file
void Writea( char *fid, VMatrix &mat);    // write ascii file
//void VMatrix::Writeb( char *fid, VMatrix &mat)   // write binary file

//////////// sorting and subsetting

VMatrix& MSort( VMatrix &a, int col, int order=0);
VMatrix& Submat( VMatrix &a, int r, int c, int lr=1, int lc=1);
VMatrix& Ch( VMatrix &a, VMatrix &b );
VMatrix& Cv( VMatrix &a, VMatrix &b );

/////////////// Specially patterned matrices

VMatrix& Tran(VMatrix &ROp );
VMatrix& Mexp(VMatrix &ROp );
VMatrix& Mabs(VMatrix &ROp );
VMatrix& Mlog(VMatrix &ROp );
VMatrix& Msqrt(VMatrix &ROp );
double Trace(VMatrix &ROp );
VMatrix& Ident(int n );
VMatrix& Helm( int n );
VMatrix& Fill(int r, int c, double d );
VMatrix& Index( int lower, int upper, int step = 1);

//////////////// Inversions and decompositions

VMatrix& Sweep( int k1, int k2, VMatrix& ROp); /* sweep out a row */
VMatrix& Inv( VMatrix &ROp ); /*matrix inversion using sweep */
VMatrix& Kron( VMatrix &a, VMatrix &b );
double Det( VMatrix &ROp );
VMatrix& Cholesky(VMatrix& ROp);
void QR( VMatrix& ROp, VMatrix& Q, VMatrix& R, boolean makeq = TTRUE);
int Svd( VMatrix &A, VMatrix &U, VMatrix &S, VMatrix &V,
          boolean makeu = TTRUE, boolean makev = TTRUE);
VMatrix& Ginv( VMatrix& ROp );
VMatrix& Fft( VMatrix &ROp, boolean INZEE = TTRUE);

/////////////// Reshaping functions

VMatrix& Vec( VMatrix& ROp );      // send matrix to vector
                                   // stack columns into an (r*c)x1
VMatrix& Vecdiag( VMatrix& ROp );  // extract diagonal into a vector
VMatrix& Diag( VMatrix& ROp );     // create a diagonal matrix from 
                                   // square matrix or column vector
VMatrix& Shape( VMatrix& ROp, int nrows = 1 );       // reshape into nrows

///////////////// accumulation 

typedef enum margins { ALL, ROWS, COLUMNS } margins; // summation margins

VMatrix& Sum( VMatrix& ROp, margins marg = ALL );    // sums over margins
VMatrix& Sumsq( VMatrix& ROp, margins marg = ALL );  // sum of squares
VMatrix& Cusum( VMatrix& ROp);                       // accumulate along rows

///////////////// max and min

VMatrix& Mmax( VMatrix& ROp, margins marg=ALL );   // matrix max
VMatrix& Mmin( VMatrix& ROp, margins marg=ALL );   //  matrix min

///////////////// elementary operations
///////////////// !!!!!!! modify input matrices

void Crow( VMatrix& ROp, int row=0, double c=1.0);
void Srow( VMatrix &ROp, int row1 = 0, int row2 = 0);
void Lrow( VMatrix &ROp, int row1 = 0, int row2 = 0, double c=0.0);
void Ccol( VMatrix& ROp, int col=0, double c=1.0);
void Scol( VMatrix &ROp, int col1 = 0, int col2 = 0);
void Lcol( VMatrix &ROp, int col1 = 0, int col2 = 0, double c=0.0);

\end{verbatim}

A boolean type is provided. The values are FFALSE and
TTRUE. They are misspelled to avoid conflicts with other
compilers which may have also defined a boolean type.

\subsection{Display Control}
Matrix elements are displayed in fields with "width"
characters, and "dec" digits to the right of the decimal
place. These parameters may be set using
Setwid()\index{functions!Setwid()} and
Setdec()\index{functions!Setdec()}. Two other functions are
provided for finding the current width and dec.
Getwid()\index{functions!Getwid()} and
Getdec()\index{functions!Getdec()} return integers for
width and dec. You should rarely require the last two
functions. They store width and dec in static variables.
They are also called by other functions needing these
values. 

\subsection{Matrix IO}
There are two formats for matrix files used in this
program: ASCII and binary. Each format have read and write
functions. Reada()\index{functions!Reada()} reads ASCII
matrices, and Writea()\index{functions!Writea()} writes
them. The format for an ASCII file is reasonably simple:
\begin{enumerate}
  \item Line 1: The first line is the matrix name. It may
        include blanks, may be no longer than 80 columns,
        and is terminated by a carriage return.
  \item Line 2: Two integers, the first is the number of
        rows and the second is the number of columns. 
  \item Lines 3...r+2 : The next r rows are the rows of the
        matrix. Columns are separated by spaces.
\end{enumerate}
There is no facility for missing data.  The data is written
to the file using "width" characters, and "dec" decimal
places, so some accuracy may be lost if you do not set the
width and decimal places before writing the matrix in ASCII
format. The argument to Reada() is a file name character
string. The file is opened and read. Reada() returns a
VMatrix reference. The arguments to Writea() are a file
name character string, and the VMatrix to be written.

Readb()\index{functions!Readb()} and
Writeb()\index{functions!Writeb()} read and write binary
matrices. Writeb() creates a binary file storing the
matrix. There is no loss in accuracy when using a binary
matrix since the binary form of the elements is written.
The file format for a binary file is:
\begin{enumerate}
  \item Bytes 1-4: An integer counting the number of
        characters in the matrix name. This is called strlen.
  \item Bytes 5...strlen+5: The matrix name string. 
  \item The next two fields are integers containing the
        numbers of rows and columns of the matrix.
  \item The remaining fields: The matrix elements stored in
        row major order.
\end{enumerate}
The arguments for Readb() and Writeb() are the same as
their ASCII counterparts. Note however that Writeb() is a
member function of the VMatrix class. This is required so
that Writeb() may be more efficient. It needs access to the
hdr pointer in its vdoub parent class. Writeb() sends data
out in blocks of size BLK\_SIZE bytes. This simplifies the
paging of the matrix blocks in vmm storage. A block is read
from the vmm and written. However, Readb() is not as
efficient. It reads elements one at a time. To use
Writeb(), you must use it as a member functions, eg.
\verb+ a.Writeb("junk.bin",a)+ where a is a VMatrix. 

Another function for loading a matrix is
LoadMat()\index{functions!LoadMat()}. LoadMat() is another
member function of VMatrix. It lets you enter matrix elements
interactively. 

\subsection{Sorting, Subsetting, and Concatenation}
MSort()\index{functions!MSort()} sorts the rows of a
matrix according to one of the columns of the matrix. The
column elements are sorted by a heap sort, then the rows of
the new matrix are stored in the order of the sorted column
elements. If the optional parameter order is supplied, and
is not zero, then columns are sorted by the elements of a
row. Note that a heap sort is used in this routine rather
than the system qsort(). Recall that qsort() only works on
vectors stored completely in RAM. 

Submat()\index{functions!Submat()} extracts a submatrix
from a matrix. The order of the arguments is a little
weird, but I wanted to take advantage of default
arguments. The bottom right corner is supplied first by the
row number, r, and the column number c. The default top
left corner is (1,1). You may supply other topmost row in
the variable lr, and left most column in the variable lc.
Submat() returns a VMatrix reference. The submatrix has top
left indices of (1,1).

Vertical and horizontal concatenation is achieved through
the functions Cv()\index{functions!Cv()} and
Ch()\index{functions!Ch()}. Both functions take two
matrices as arguments. The matrices must have the same
number of columns to be vertically concatenated, and the
same number of rows to be horizontally concatenated. Both
functions return a VMatrix reference. 

\subsection{Elementwise Operators}
Four elementwise functions are provided: Mexp(),
Mabs(), Msqrt(), and
Mlog().\index{functions!Mexp(),Mabs(),Msqrt(),Mlog()}.
These take the elementwise exponent, absolute value, square
root, and natural logarithm. Argument ranges are checked in
the square root and logarithm functions. If an invalid
argument is given, the program terminates through the error
function, Dispatch-$>$Nrerror(). This function was defined a
returning a double so that it could be used in the ternary
operator as an error branch. The functions return
references to a VMatrix.

\section{Math Functions}
\subsection{Tran()}
Tran()\index{functions!Tran()} returns the transpose of a
matrix.\footnote{This is one of my disappointments with C++. I
would have preferred to overload a postfix operator, say `,
for the transpose, but such an operator does not exist in
C++. So what, I never intended to write APL.}

\subsection{Trace()}
Trace()\index{functions!Trace()} returns the trace of a
square matrix, i.e. the sum of its diagonal elements.

\subsection{Ident(), Helm(), Fill(), Index()}
These three functions create specially patterned matrices.
Ident()\index{functions!Ident()} creates an nxn identity
matrix. 

Helm()\index{functions!Helm()} creates an nxn Helmert
matrix. It is useful in mixed linear models. It is an
orthonormal matrix in the sense that \verb+ H'H=HH'=I+. The
elements of an nxn Helmert matrix are given by
$$
 \left[
\begin{array}{ccccc}
{1 \over \sqrt{n}}&{1 \over \sqrt{2}}&{1\over \sqrt{3\cdot 2}}&\ldots&
   {1 \over \sqrt{ n(n-1) }} \\
{1 \over \sqrt{n}}&{-1\over \sqrt{2}}&{1\over \sqrt{3\cdot 2}}&\ldots&
   {1 \over \sqrt{ n(n-1) }} \\
{1 \over \sqrt{n}}& 0 &{-2\over \sqrt{3\cdot 2}}&\ldots&{1 \over \sqrt{ n(n-1)
}} \\
\vdots & \vdots & 0 & \ldots & {1\over \sqrt{n(n-1)}}\\
\vdots & \vdots & \vdots & \ldots & {-(n-1)\over \sqrt{n(n-1)}}
\end{array}
\right]
$$

Fill()\index{functions!Fill()} creates an $rxc$ matrix of
elements of $d$. It returns a VMatrix reference. The
Helmert matrix is a set of eigenvectors of the matrix 

\begin{eqnarray*}
V & = & \phi_0 \hbox{Ident}(n) + \hbox{Fill}(n,n,\phi_1) \\
  & = & \hbox{Helm}(n) 
        \hbox{Diag}_n(\phi_0 + n\phi_1, \phi_0, \ldots, \phi_0)
        \hbox{Tran}( \hbox{Helm} (n)) 
\end{eqnarray*}

Index()\index{functions!Index()} creates index matrices. It
simplifies making a column whose elements range from \verb+lower+ to
\verb+upper+ by increments of \verb+step+. If \verb+lower+ is
greater than \verb+upper+, then \verb+step+ should be negative.
The default step size is 1.

\subsection{Sweep() and Inv()}
Sweep()\index{functions!Sweep()} may be used for creating a
generalized inverse of a square matrix using Gaussian
elimination. A generalized inverse of the matrix $A$, say $A^+$
satisfies the properties, $A^+AA^+=A^+$ and $AA^+A=A$.
Also, the quantity $A^+A$ is an idempotent matrix of zeros
and ones, so its rank is its trace. These are handy
properties when doing linear models. 

The arguments of Sweep() are the rows to be eliminated,
from k1 to k2, and the matrix to operate on. If all rows
are eliminated, and the matrix has full rank, then the
matrix is inverted. The function,
Inv()\index{functions!Inv()}, performs matrix inversion by
eliminating all rows and columns. If any pivot is less than
1.0E-8, then the row and column being eliminated is set to
zero. This procedure creates $A^+$. 

Sweep() also simplifies the calculation of regression
coefficients and the residual mean square. Consider the
regression of $X_{Nxp}$ on $Y_{Nx1}$. The steps for
calculating the regression coefficients are as follows:
\begin{enumerate}
\item Concatenate $X$ and $Y$ using the statement
      VMatrix XY = Ch(X,Y); This produces an Nx(p+1) matrix.
\item Calculate the matrix of crossproducts using
      VMatrix xypxy = Tran(XY)*XY;
\item Calculate the coefficients by sweeping out the first p
      rows of xypxy, VMatrix ixypxy = Sweep(1,p,xypxy);
\item The regression coefficients are in the right most
      column: VMatrix betahat = Submat(ixypxy,p,p+1,1,p+1);
\item The mean square error is calculated from the bottom
      right corner: double mse = ixypxy.m(p+1,p+1) /
      ((double)(X.r-X.c));
\end{enumerate}
This results from sweeping out the columns of $X'X$. The
original matrix is partitioned as
$$
xypxy = 
\left[
\begin{array}{cc}
X'X & X'Y \\
Y'X & Y'Y
\end{array}
\right]
$$
The swept matrix is given by
$$
\hbox{Sweep}(1,p,xpyxy) =
\left[
\begin{array}{cc}
(X'X)^{-1} & (X'X)^{-1}X'Y \\
-Y'X(X'X)^{-1} & Y'(I-X(X'X)^{-1}X')Y
\end{array}
\right]
$$
The top right matrix is the solution to the least squares
normal equations. More information about the sweep operator
can be found in Kennedy and Gentle\cite{KG:sc},
Heiberger\cite{HE:ad} and Hocking\cite{HK:am}.

\subsection{Kron()}
Kron()\index{functions!Kron()} calculates the Kroniker
product of two matrices, $C=A\otimes B$. The (i,j) block of
$C$ is $a_{i,j}B$. This is a very useful operator in
balanced mixed linear models. See Hocking\cite{HK:am} for
these uses, and Gram\cite{Gm:kp} for its applications in
matrix calculus. 

\subsection{Det()}
Det()\index{functions!Det()} calculates the determinant of a
square matrix. It uses pivoting to move small pivots to the
right as they are encountered. It uses Gaussian elimination to
convert the matrix to an upper triangular matrix. The
determinant is the product of the diagonal elements.

\subsection{Cholesky(), QR()}
Cholesky()\index{functions!Cholesky()} converts a square
symmetric matrix, $A$ into aan upper triangular matrix, $S$, such
that $A=S'S$. The function does not pivot, and stops if a zero
pivot is found. This version is taken from Hocking\cite{HK:am}. 

QR()\index{functions!QR()} converts an $rxc$ matrix $A$ into the
product $QR$, where $Q$ is $rxc$ of orthogonal columns, and $R$
is an upper triangular $cxc$ matrix. The routine stops if a zero
diagonal element of $R$ is produced. The function produces $Q$
by default, but supplying makeq=FFALSE turns this option off.
Typically the QR reduction is used to reduce a $X$ matrix in
regression to avoid round off error in calculating $X'X$. The
$X$ matrix has at least as many rows as columns, so there is
some savings in calculating $Q$ as $Q = XR^{-1}$. You should
turn off the calculation of $Q$ if finding $R^{-1}$ is
prohibitive, or if only $R$ is needed. Its use in regression may
be improved by pivoting zero diagonal elements of R. (See
\cite{KBH:nm} for more). The program is a
translation of an Algol program found in \cite{SB:na}.
 
\subsection{Svd() and Ginv()}
Svd()\index{functions!Svd()} provides a singular valued
decomposition. It is a direct translation of the EISPACK routine
with the same name (\cite{SBG:ei}). The routine decomposes an
$rxc$ matrix $A$ as $A=UDV'$, where $U'U=UU'=I_r$,
$V'V=VV'=I_c$, and a vector $S$ containing the diagonal elements
of $D$, which are the square roots of the eigenvalues of $A'A$.
The default is to produce $U$ and $V$, but this may be turned
off by supplying makev=FFALSE or makeu=FFALSE.

Ginv()\index{functions!Ginv()} produces a generalized inverse of
$A$ using the singular valued decomposition: $A^+ = VD^+U'$. The
matrix $D^+$ is a diagonal matrix with the reciprocals of the
nonzero singular values of $A$.

\subsection{An Almost Fft()}
Fft()\index{functions!Fft()} calculates the Fourier transform
and its inverse for a vector of $n$ reals, or a complex vector
stored in a nx2 matrix. This Fft does not require n to be a
power of 2. However it slows down for n$>$1000, or for n with
a smallest prime factor greater than 19. In such cases it may be
better to resort to a $2^n$ fft and zero padding (see
\cite{Pr:nrc}). You can calculate the inverse fft by supplying
the variable INZEE=FFALSE. See the example program \verb+
testreg.cpp+ for calculating the first 5 serial correlations
using the fft method.

The program is a translation of a FORTRAN program found in Conte
and de Boor( see \cite{CB:en} ). The translation was produced by
R. G. Davies in his matrix class, NEWMAT currently available in
the Borland C++ library of Compuserve. I have quite unashamedly
translated it for use in YAMP. I have, however, added the
ability to do the inverse fft, so I don't feel like a total
thief. My version is slower since it calculates offsets in
vectors rather than incrementing pointers across the heap. My
version is also slower in the virtual memory model since it must
page memory almost constantly.  My guess is that explicit
evaluation of the discrete Fourier transform is faster since it
requires less memory paging. I eventually need to see what can
be done about this.

\subsection{Reshaping functions: Vec(), Vecdiag(), Diag(), and
Shape()} These functions reshape matrices.
Vec()\index{functions!Vec()} sends a matrix to a vector. The
columns of the vector are stacked into a $rcx1$ vector.
Vecdiag()\index{functions!Vecdiag()} extracts the principal
diagonal of a square matrix into a vector.
Diag()\index{functions!Diag()} sends a vector into the principal
diagonal of a square zero matrix. It also zeros the off-diagonal
elements if the input matrix is square.

Shape()\index{functions!Shape()} takes an $rxc$ matrix to a new matrix
with $nrows$ rows. The value of $nrows$ must divide $rc$ without a
remainder or the program will stop. 

\subsection{Summation Functions: Sum(), Sumsq(), Cusum()}
These functions form different kinds of sums of the elements of
a matrix. Sum()\index{functions!Sum()} sums the elements of a
matrix. The default is to sum all of the elements and return a
1x1 matrix. The optional parameter marg may be set to ROWS or
COLUMNS, to sum across rows or columns respectively. Setting
marg=ROWS produces a 1xc matrix of column sums. Setting
marg=COLUMNS produces an rx1 matrix of row sums.

Sumsq()\index{functions!Sumsq()} functions as does Sum(), but
produces sums of squared elements.

Cusum()\index{functions!Cusum()} forms a running total across
rows. It procedes row-wise through a matrix and replaces each
element by the current accumulated total.

\subsection{Mmin(), Mmax()}
Mmin() and Mmax()\index{functions!Mmin(), Mmax()} calculate the
min and max of matrices. They also return the index of values:
\begin{verbatim}
   if marg = ALL      return 3x1   [ maxr, maxc, max element ]
   if marg = ROWS     return 3xc   Tran([ maxr, col=1..c, max element ])
   if marg = COLUMNS  return rx3   [ row=1..r, maxc, max element]
\end{verbatim}
where maxr is the index of the row max, maxc is the index of the
column max, and max element is the maximum. Mmin() functions
similarly,  but returns mininums and their indices.

\subsection{Elementary Row and Column Operations}
The elementary row and column operations operate directly on a
matrix. That is, they alter the input matrix, instead of
returning a new matrix that results from the operations. The
functions are:
\begin{verbatim}
       Crow( VMatrix& ROp, int row=0, double c=1.0);
       Srow( VMatrix &ROp, int row1 = 0, int row2 = 0);
       Lrow( VMatrix &ROp, int row1 = 0, int row2 = 0, double c=0.0);
       Ccol( VMatrix& ROp, int col=0, double c=1.0);
       Scol( VMatrix &ROp, int col1 = 0, int col2 = 0);
       Lcol( VMatrix &ROp, int col1 = 0, int col2 = 0, double c=0.0);
\end{verbatim}
Crow()(Ccol())\index{functions!Crow(), Ccol()} multiply the
row(col) by the constant c.
Srow()(Scol())\index{functions!Srow(),Scol()} swap row1 with
row2 (col1 with col2).
Lrow()(Lcol())\index{functions!Lrow(),Lcol()} add a constant
multiple of row1 to row2.

Note that the functions supply out of range default arguments.
This forces you to supply all of the values correctly. The default
constants perform no action, so they must also be supplied.
