\documentclass[10pt]{article}
%\VignetteIndexEntry{RcppEigen-intro}
\usepackage[margin=1in]{geometry}
\usepackage{color, alltt, bm, amsmath, listings}
\lstset{language=C++,basicstyle=\small}
\usepackage[authoryear,round,longnamesfirst]{natbib}
\usepackage[colorlinks]{hyperref}
\definecolor{link}{rgb}{0,0,0.3}	%% next few lines courtesy of RJournal.sty
\hypersetup{
    colorlinks,%
    citecolor=link,%
    filecolor=link,%
    linkcolor=link,%
    urlcolor=link
}
\newcommand{\proglang}[1]{\textsf{#1}}
\newcommand{\pkg}[1]{{\fontseries{b}\selectfont #1}}
\newcommand{\code}[1]{\texttt{#1}}
\newcommand{\rank}{\operatorname{rank}}
<<version,echo=FALSE,print=FALSE>>=
prettyVersion <- packageDescription("RcppEigen")$Version
prettyDate <- format(Sys.Date(), "%B %e, %Y")
@
\author{Douglas Bates}
\title{An Introduction to \pkg{RcppEigen}}
\date{\pkg{RcppEigen} version \Sexpr{prettyVersion} as of \Sexpr{prettyDate}}

<<preliminaries,echo=FALSE>>=
link <- function( f, package, text = f, root = "http://finzi.psych.upenn.edu/R/library/" ){
	h <- if( missing(package) ) {
		as.character( help( f ) )
	} else {
		as.character( help( f, package = paste( package, sep = "" ) ) )
	}
	if( ! length(h) ){
		sprintf( "\\\\textbf{%s}", f )
	} else {
		rx <- "^.*/([^/]*?)/help/(.*?)$"
		package <- sub( rx, "\\1", h, perl = TRUE )
		page <- sub( rx, "\\2", h, perl = TRUE )
		sprintf( "\\\\href{%s%s/html/%s.html}{\\\\texttt{%s}}", root, package, page, text )
	}
}
linkS4class <- function( cl, package, text = cl, root = "http://finzi.psych.upenn.edu/R/library/" ){
	link( sprintf("%s-class", cl), package, text, root )
}
require( inline )
require( RcppEigen )
@
\begin{document}
\maketitle

\abstract{
  \noindent
  The \pkg{RcppEigen} package provides access from \proglang{R} to the
  \pkg{Eigen} \proglang{C++} template library for numerical linear
  algebra. \pkg{Rcpp} \citep{JSS:Rcpp} classes and specializations of
  the \proglang{C++} templated functions \code{as} and \code{wrap}
  from \pkg{Rcpp} provide the ``glue'' for passing objects from
  \proglang{R} to \proglang{C++} and back.  }

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

As stated in the \pkg{Rcpp} \citep{CRAN:Rcpp} vignette, ``Extending \pkg{Rcpp}''
\begin{quote}
  \pkg{Rcpp} facilitates data interchange between \proglang{R} and
  \proglang{C++} through the templated functions \texttt{Rcpp::as} (for
  conversion of objects from \proglang{R} to \proglang{C++}) and
  \texttt{Rcpp::wrap} (for conversion from \proglang{C++} to \proglang{R}).
\end{quote}
The \pkg{RcppEigen} package provides the header files composing the
\pkg{Eigen} \proglang{C++} template library and implementations of
\texttt{Rcpp::as} and \texttt{Rcpp::wrap} for the \proglang{C++}
classes defined in \pkg{Eigen}.

The \pkg{Eigen} classes themselves provide high-performance,
versatile and comprehensive representations of dense and sparse
matrices and vectors, as well as decompositions and other functions
to be applied to these objects.  In the next section we introduce some
of these classes and show how to interface to them from R.

\section{Eigen classes}
\label{sec:eclasses}

\pkg{Eigen} (\url{http://eigen.tuxfamily.org}) is a \proglang{C++} template
library providing classes for many forms of matrices, vectors, arrays
and decompositions.  These classes are flexible and comprehensive
allowing for both high performance and well structured code
representing high-level operations. \proglang{C++} code based on Eigen
is often more like \proglang{R} code, working on the ``whole object'',
than compiled code in other languages where operations often must be
coded in loops.

As in many \proglang{C++} template libraries using template
meta-programming
\citep{Abrahams+Gurtovoy:2004:TemplateMetaprogramming}, the templates
themselves can be very complicated.  However, \pkg{Eigen} provides
typedef's for common classes that correspond to \proglang{R} matrices and
vectors, as shown in Table~\ref{tab:REigen}. We will use these
typedef's throughout this document.
\begin{table}[tb]
  \caption{Correspondence between R matrix and vector types and classes in the \code{Eigen} namespace.}
  \label{tab:REigen}
  \centering 
  \begin{tabular}{l l}
    \hline
    \multicolumn{1}{c}{\proglang{R} object type} & \multicolumn{1}{c}{\pkg{Eigen} class typedef}\\
    \hline
    numeric matrix & \code{MatrixXd}\\
    integer matrix & \code{MatrixXi}\\
    complex matrix & \code{MatrixXcd}\\
    numeric vector & \code{VectorXd}\\
    integer vector & \code{VectorXi}\\
    complex vector & \code{VectorXcd}\\
    \code{Matrix::dgCMatrix} & \code{SparseMatrix<double>}\\
    \hline
  \end{tabular}
\end{table}

The \proglang{C++} classes shown in Table~\ref{tab:REigen} are in the
\code{Eigen} namespace, which means that they must be written as
\code{Eigen::MatrixXd}.  However, if we preface our use of these class
names with a declaration like
\begin{lstlisting}[language=C++]
using Eigen::MatrixXd;
\end{lstlisting}
we can use these names without the qualifier.

\subsection{Mapped matrices in Eigen}
\label{sec:mapped}

Storage for the contents of matrices from the classes shown in
Table~\ref{tab:REigen} is allocated and controlled by the class
constructors and destructors.  Creating an instance of such a class
from an \proglang{R} object involves copying its contents.  An
alternative is to have the contents of the \proglang{R} matrix or
vector mapped to the contents of the object from the Eigen class.  For
dense matrices we use the Eigen templated class \code{Map}.  For
sparse matrices we use the Eigen templated class \code{MappedSparseMatrix}.

We must, of course, be careful not to modify the contents of the
\proglang{R} object in the \proglang{C++} code.  A recommended
practice is always to declare mapped objects as \lstinline!const!.

\subsection{Arrays in Eigen}
\label{sec:arrays}

For matrix and vector classes \pkg{Eigen} overloads the \texttt{`*'}
operator to indicate matrix multiplication.  Occasionally we want
component-wise operations instead of matrix operations.  The
\code{Array} templated classes are used in \pkg{Eigen} for
component-wise operations.  Most often we use the \code{array()} method
for Matrix or Vector objects to create the array.  On those occasions
when we wish to convert an array to a matrix or vector object we use
the \code{matrix()} method.

\subsection{Structured matrices in \pkg{Eigen}}
\label{sec:structured}

There are \pkg{Eigen} classes for matrices with special structure such
as symmetric matrices, triangular matrices and banded matrices.  For
dense matrices, these special structures are described as ``views'',
meaning that the full dense matrix is stored but only part of the
matrix is used in operations.  For a symmetric matrix we need to
specify whether the lower triangle or the upper triangle is to be used as
the contents, with the other triangle defined by the implicit symmetry.


\section{Some simple examples}
\label{sec:simple}

\proglang{C++} functions to perform simple operations on matrices or
vectors can follow a pattern of:
\begin{enumerate}
\item Map the \proglang{R} objects passed as arguments into Eigen objects.
\item Create the result.
\item Return \code{Rcpp::wrap} applied to the result.
\end{enumerate}

An idiom for the first step is
\begin{lstlisting}[language=C++]
using Eigen::Map;
using Eigen::MatrixXd;
using Rcpp::as;

const Map<MatrixXd>  A(as<Map<MatrixXd> >(AA));  
\end{lstlisting}
where \code{AA} is the name of the R object (called an \code{SEXP} in
\proglang{C} and \proglang{C++}) passed to the \proglang{C++} function.

The \Sexpr{link("cxxfunction")} from the \pkg{inline}
\citep{CRAN:inline} package for \proglang{R} and its \pkg{RcppEigen}
plugin provide a convenient way of developing and debugging the
\proglang{C++} code.  For actual production code we generally
incorporate the \proglang{C++} source code files in a package and
include the line \code{LinkingTo: Rcpp, RcppEigen} in the package's
\code{DESCRIPTION} file.  The
\Sexpr{link("RcppEigen.package.skeleton")} function provides a quick
way of generating the skeleton of a package using \pkg{RcppEigen}
facilities.

The \code{cxxfunction} with the \code{"Rcpp"} or \code{"RcppEigen"}
plugins has the \code{as} and \code{wrap} functions already defined as
\code{Rcpp::as} and \code{Rcpp::wrap}.  In the examples below we will
omit these declarations.  Do remember that you will need them in
\proglang{C++} source code for a package.

The first few examples are simply for illustration as the operations
shown could be more effectively performed directly in \proglang{R}.
We do compare the results from \pkg{Eigen} to those from the direct
\proglang{R} results.

\subsection{Transpose of an integer matrix}
\label{sec:transpose}

We create a simple matrix of integers
<<Adef>>=
(A <- matrix(1:6, ncol=2))
str(A)
@ 
and, in Listing~\ref{trans}, use the \code{transpose()} method for the
\code{Eigen::MatrixXi} class to return its transpose. The \proglang{R}
matrix in the \code{SEXP} \code{AA} is mapped to an
\code{Eigen::MatrixXi} object then the matrix \code{At} is constructed
from its transpose and returned to \proglang{R}.  We check that it
works as intended.

<<transCpp,echo=FALSE>>=
transCpp <-'
using Eigen::Map;
using Eigen::MatrixXi;
                 // Map the integer matrix AA from R
const Map<MatrixXi>  A(as<Map<MatrixXi> >(AA));
                 // evaluate and return the transpose of A
const MatrixXi      At(A.transpose());
return wrap(At);
'
@
\begin{lstlisting}[frame=tb,float,caption={transCpp: Transpose a matrix of integers},label=trans]
<<transCppLst,results=tex,echo=FALSE>>=
cat(transCpp, "\n")
@ 
\end{lstlisting}
<<ftrans>>=
ftrans <- cxxfunction(signature(AA="matrix"), transCpp, plugin="RcppEigen")
(At <- ftrans(A))
stopifnot(all.equal(At, t(A)))
@
For numeric or integer matrices the \code{adjoint()} method is
equivalent to the \code{transpose()} method.  For complex matrices, the
adjoint is the conjugate of the transpose.  In keeping with the
conventions in the \pkg{Eigen} documentation we will, in what follows,
use the \code{adjoint()} method to create the transpose of numeric or
integer matrices.

\subsection{Products and cross-products}
\label{sec:products}

As mentioned in Sec.~\ref{sec:arrays}, the \code{`*'} operator
performs matrix multiplication on \code{Eigen::Matrix} or
\code{Eigen::Vector} objects. The \proglang{C++} code in
Listing~\ref{prod} produces
<<prodCpp,echo=FALSE>>=
prodCpp <- '
using Eigen::Map;
using Eigen::MatrixXi;
const Map<MatrixXi>    B(as<Map<MatrixXi> >(BB));
const Map<MatrixXi>    C(as<Map<MatrixXi> >(CC));
return List::create(_["B %*% C"]         = B * C,
                    _["crossprod(B, C)"] = B.adjoint() * C);
'
@
\begin{lstlisting}[frame=tb,float,caption={prodCpp: Product and cross-product of two matrices},label=prod]
<<prodCppLst,results=tex,echo=FALSE>>=
cat(prodCpp, "\n")
@ 
\end{lstlisting}
<<prod>>=
fprod <- cxxfunction(signature(BB = "matrix", CC = "matrix"), prodCpp, "RcppEigen")
B <- matrix(1:4, ncol=2); C <- matrix(6:1, nrow=2)
str(fp <- fprod(B, C))
stopifnot(all.equal(fp[[1]], B %*% C), all.equal(fp[[2]], crossprod(B, C)))
@

Notice that the \code{create} method for the \pkg{Rcpp} class
\code{List} implicitly applies \code{Rcpp::wrap} to its arguments.

\subsection{Crossproduct of a single matrix}
\label{sec:crossproduct}

As shown in the last example, the \proglang{R} function
\code{crossprod} calculates the product of the transpose of its first
argument with its second argument.  The single argument form,
\code{crossprod(X)}, evaluates $\bm X^\prime\bm X$.  We could, of
course, calculate this product as
<<eval=FALSE>>=
t(X) %*% X
@
but \code{crossprod(X)} is roughly twice as fast because the result is
known to be symmetric and only one triangle needs to be calculated.
The function \code{tcrossprod} evaluates \code{crossprod(t(X))}
without actually forming the transpose.

To express these calculations in Eigen we create a \code{SelfAdjointView},
which is a dense matrix of which only one triangle is used, the other
triangle being inferred from the symmetry.  (``Self-adjoint'' is
equivalent to symmetric for non-complex matrices.)

The \pkg{Eigen} class name is \code{SelfAdjointView}.  The method for
general matrices that produces such a view is called
\code{selfadjointView}.  Both require specification of either the
\code{Lower} or \code{Upper} triangle.

For triangular matrices the class is \code{TriangularView} and the
method is \code{triangularView}.  The triangle can be specified as
\code{Lower}, \code{UnitLower}, \code{StrictlyLower}, \code{Upper},
\code{UnitUpper} or \code{StrictlyUpper}.

For self-adjoint views the \code{rankUpdate} method adds a scalar multiple
of $\bm A\bm A^\prime$ to the current symmetric matrix.  The scalar
multiple defaults to 1.  The code in Listing~\ref{crossprod} produces

<<crossprod,echo=FALSE>>=
crossprodCpp <- '
using Eigen::Map;
using Eigen::MatrixXi;
using Eigen::Lower;

const Map<MatrixXi> A(as<Map<MatrixXi> >(AA));
const int           m(A.rows()), n(A.cols());                        
MatrixXi          AtA(MatrixXi(n, n).setZero().
                      selfadjointView<Lower>().rankUpdate(A.adjoint()));
MatrixXi          AAt(MatrixXi(m, m).setZero().
                      selfadjointView<Lower>().rankUpdate(A));

return List::create(_["crossprod(A)"]  = AtA,
                    _["tcrossprod(A)"] = AAt);
'
@
\begin{lstlisting}[frame=tb,float,caption={crossprodCpp: Cross-product and transposed cross-product of a single matrix},label=crossprod]
<<crossprodCppLst,results=tex,echo=FALSE>>=
cat(crossprodCpp, "\n")
@ 
\end{lstlisting}
<<>>=
fcprd <- cxxfunction(signature(AA = "matrix"), crossprodCpp, "RcppEigen")
str(crp <- fcprd(A))
stopifnot(all.equal(crp[[1]], crossprod(A)), all.equal(crp[[2]], tcrossprod(A)))
@

To some, the expressions to construct \code{AtA} and \code{AAt} in
that code fragment are compact and elegant.  To others they are
hopelessly confusing.  If you find yourself in the latter group, you
just need to read the expression left to right.  So, for example, we
construct \code{AAt} by creating a general integer matrix of size
$m\times m$ (where $\bm A$ is $m\times n$), ensure that all its
elements are zero, regard it as a self-adjoint (i.e. symmetric) matrix
using the elements in the lower triangle, then add $\bm A\bm A^\prime$
to it and convert back to a general matrix form (i.e. the strict
lower triangle is copied into the strict upper triangle).

For these products we could use either the lower triangle or the upper
triangle as the result will be symmetrized before it is returned.

\subsection{Cholesky decomposition of the crossprod}
\label{sec:chol}

The Cholesky decomposition of the positive-definite, symmetric matrix,
$\bm A$, can be written in several forms.  Numerical analysts define
the ``LLt'' form as the lower triangular matrix, $\bm L$, such that
$\bm A=\bm L\bm L^\prime$ and the ``LDLt'' form as a unit lower
triangular matrix $\bm L$ and a diagonal matrix $\bm D$ with positive
diagonal elements such that $\bm A=\bm L\bm D\bm L^\prime$.
Statisticians often write the decomposition as $\bm A=\bm R^\prime\bm
R$ where $\bm R$ is an upper triangular matrix.  Of course, this $\bm
R$ is simply the transpose of $\bm L$ from the ``LLt'' form.

The templated \pkg{Eigen} classes for the LLt and LDLt forms are
called \code{LLT} and \code{LDLT}.  In general we would preserve the
objects from these classes so that we could use them for solutions of
linear systems.  For illustration we simply return the matrix $\bm L$
from the ``LLt'' form.

Because the Cholesky decomposition involves taking square roots we
switch to numeric matrices
<<storage>>=
storage.mode(A) <- "double"
@ 
before applying the code in Listing~\ref{chol}.
<<cholCpp,echo=FALSE>>=
cholCpp <- '
using Eigen::Map;
using Eigen::MatrixXd;
using Eigen::LLT;
using Eigen::Lower;

const Map<MatrixXd>   A(as<Map<MatrixXd> >(AA));
const int             n(A.cols());
const LLT<MatrixXd> llt(MatrixXd(n, n).setZero().
                        selfadjointView<Lower>().rankUpdate(A.adjoint()));

return List::create(_["L"] = MatrixXd(llt.matrixL()),
                    _["R"] = MatrixXd(llt.matrixU()));
'
@ 
\begin{lstlisting}[frame=tb,float,caption={cholCpp: Cholesky decomposition of a cross-product},label=chol]
<<cholCppLst,results=tex,echo=FALSE>>=
cat(cholCpp, "\n")
@ 
\end{lstlisting}
<<fchol>>=
fchol <- cxxfunction(signature(AA = "matrix"), cholCpp, "RcppEigen")
(ll <- fchol(A))
stopifnot(all.equal(ll[[2]], chol(crossprod(A))))
@

\subsection{Determinant of the cross-product matrix}
\label{sec:determinant}

The ``D-optimal'' criterion for experimental design chooses the design
that maximizes the determinant, $|\bm X^\prime\bm X|$, for the
$n\times p$ model matrix (or Jacobian matrix), $\bm X$.  The
determinant, $|\bm L|$, of the $p\times p$ lower Cholesky factor
$\bm L$, defined so that $\bm L\bm L^\prime=\bm X^\prime\bm X$, is
the product of its diagonal elements, as is the case for any
triangular matrix.  By the properties of determinants,
\begin{displaymath}
  |\bm X^\prime\bm X|=|\bm L\bm L^\prime|=|\bm L|\,|\bm L^\prime|=|\bm L|^2
\end{displaymath}

Alternatively, if we use the ``LDLt'' decomposition, $\bm L\bm D\bm
L^\prime=\bm X^\prime\bm X$ where $\bm L$ is unit lower triangular and
$\bm D$ is diagonal then $|\bm X^\prime\bm X|$ is the product of the
diagonal elements of $\bm D$.  Because we know that the diagonals of
$\bm D$ must be non-negative, we often evaluate the logarithm of the
determinant as the sum of the logarithms of the diagonal elements of
$\bm D$.  Several options are shown in Listing~\ref{cholDet}.

<<echo=FALSE>>=
cholDetCpp <- '
using Eigen::Lower;
using Eigen::Map;
using Eigen::MatrixXd;
using Eigen::VectorXd;

const Map<MatrixXd>   A(as<Map<MatrixXd> >(AA));
const int             n(A.cols());
const MatrixXd      AtA(MatrixXd(n, n).setZero().
                        selfadjointView<Lower>().rankUpdate(A.adjoint()));
const MatrixXd     Lmat(AtA.llt().matrixL());
const double       detL(Lmat.diagonal().prod());
const VectorXd     Dvec(AtA.ldlt().vectorD());

return List::create(_["d1"] = detL * detL,
                    _["d2"] = Dvec.prod(),
                    _["ld"] = Dvec.array().log().sum());
'
@ 
\begin{lstlisting}[frame=tb,float,caption={cholDetCpp: Determinant of a cross-product using the Cholesky decomposition},label=cholDet]
<<cholDetCppLst,results=tex,echo=FALSE>>=
cat(cholDetCpp, "\n")
@ 
\end{lstlisting}
<<fdet>>=
fdet <- cxxfunction(signature(AA = "matrix"), cholDetCpp, "RcppEigen")
unlist(ll <- fdet(A))
@

Note the use of the \code{array()} method in the calculation of the
log-determinant.  Because the \code{log()} method applies to arrays,
not to vectors or matrices, we must create an array from \code{Dvec}
before applying the \code{log()} method.

\section{Least squares solutions}
\label{sec:leastSquares}

A common operation in statistical computing is calculating a least
squares solution, $\widehat{\bm\beta}$, defined as
\begin{displaymath}
  \widehat{\bm\beta}=\arg\min_{\beta}\|\bm y-\bm X\bm\beta\|^2
\end{displaymath}
where the model matrix, $\bm X$, is $n\times p$ ($n\ge p$) and $\bm y$
is an $n$-dimensional response vector.  There are several ways based
on matrix decompositions, to determine such a solution.  We have
already seen two forms of the Cholesky decomposition: ``LLt'' and
``LDLt'', that can be used to solve for $\widehat{\bm\beta}$.  Other
decompositions that can be used are the QR decomposition, with or
without column pivoting, the singular value decomposition and the
eigendecomposition of a symmetric matrix.

Determining a least squares solution is relatively straightforward.
However, in statistical computing we often require additional information,
such as the standard errors of the coefficient estimates.  Calculating
these involves evaluating the diagonal elements of $\left(\bm
  X^\prime\bm X\right)^{-1}$ and the residual sum of squares, $\|\bm
y-\bm X\widehat{\bm\beta}\|^2$.

\subsection{Least squares using the ``LLt'' Cholesky}
\label{sec:LLtLeastSquares}

Listing~\ref{lltLS}
<<lltLSCpp,echo=FALSE>>=
lltLSCpp <- '
using Eigen::LLT;
using Eigen::Lower;
using Eigen::Map;
using Eigen::MatrixXd;
using Eigen::VectorXd;

const Map<MatrixXd>   X(as<Map<MatrixXd> >(XX));
const Map<VectorXd>   y(as<Map<VectorXd> >(yy));
const int             n(X.rows()), p(X.cols());
const LLT<MatrixXd> llt(MatrixXd(p, p).setZero().
                        selfadjointView<Lower>().rankUpdate(X.adjoint()));
const VectorXd  betahat(llt.solve(X.adjoint() * y));
const VectorXd   fitted(X * betahat);
const VectorXd    resid(y - fitted);
const int            df(n - p);
const double          s(resid.norm() / std::sqrt(double(df)));
const VectorXd       se(s * llt.matrixL().solve(MatrixXd::Identity(p, p)).
                        colwise().norm());
return     List::create(_["coefficients"]   = betahat,
                        _["fitted.values"]  = fitted,
                        _["residuals"]      = resid,
                        _["s"]              = s,
                        _["df.residual"]    = df,
                        _["rank"]           = p,
                        _["Std. Error"]     = se);
'
@ 
\begin{lstlisting}[frame=tb,float,caption={lltLSCpp: Least squares using the Cholesky decomposition},label=lltLS]
<<lltLSCppLst,results=tex,echo=FALSE>>=
cat(lltLSCpp, "\n")
@ 
\end{lstlisting}
shows a calculation of the least squares coefficient estimates
(\code{betahat}) and the standard errors (\code{se}) through an
``LLt'' Cholesky decomposition of the crossproduct of the model
matrix, $\bm X$.  We check that the results from this calculation do
correspond to those from the \code{lm.fit} function in \proglang{R}
(\code{lm.fit} is the workhorse function called by \code{lm} once the
model matrix and response have been evaluated).
<<lltLS>>=
lltLS <- cxxfunction(signature(XX = "matrix", yy = "numeric"), lltLSCpp, "RcppEigen")
data(trees, package="datasets")
str(lltFit <- with(trees, lltLS(cbind(1, log(Girth)), log(Volume))))
str(lmFit <- with(trees, lm.fit(cbind(1, log(Girth)), log(Volume))))
for (nm in c("coefficients", "residuals", "fitted.values", "rank"))
    stopifnot(all.equal(lltFit[[nm]], unname(lmFit[[nm]])))
stopifnot(all.equal(lltFit[["Std. Error"]],
                    unname(coef(summary(lm(log(Volume) ~ log(Girth), trees)))[,2])))
@

There are several aspects of the \proglang{C++} code in
Listing~\ref{lltLS} worth mentioning.  The \code{solve} method for the
\code{LLT} object evaluates, in this case, $\left(\bm X^\prime\bm
  X\right)^{-1}\bm X^\prime\bm y$ but without actually evaluating the
inverse.  The calculation of the residuals, $\bm y-\widehat{\bm y}$,
can be written, as in \proglang{R}, as \code{y - fitted}. (But note
that \pkg{Eigen} classes do not have a ``recycling rule as in
\proglang{R}.  That is, the two vector operands must have the same
length.)  The \code{norm()} method evaluates the square root of the
sum of squares of the elements of a vector.  Although we don't
explicitly evaluate $\left(\bm X^\prime\bm X\right)^{-1}$ we do
evaluate $\bm L^{-1}$ to obtain the standard errors.  Note also the
use of the \code{colwise()} method in the evaluation of the standard
errors.  It applies a method to the columns of a matrix, returning a
vector.  The \pkg{Eigen} \code{colwise()} and \code{rowwise()} methods
are similar in effect to the \code{apply} function in \proglang{R}.

In the descriptions of other methods for solving least squares
problems, much of the code parallels that shown in
Listing~\ref{lltLS}.  We will omit the redundant parts and show only
the evaluation of the coefficients, the rank and the standard errors.
Actually, we only calculate the standard errors up to the scalar
multiple of $s$, the residual standard error, in these code fragments.
The calculation of the residuals and $s$ and the scaling of the
coefficient standard errors is the same for all methods.  (See the
files \code{fastLm.h} and \code{fastLm.cpp} in the \pkg{RcppEigen}
source package for details.)

\subsection{Least squares using the unpivoted QR decomposition}
\label{sec:QR}

A QR decomposition has the form
\begin{displaymath}
  \bm X=\bm Q\bm R=\bm Q_1\bm R_1
\end{displaymath}
where $\bm Q$ is an $n\times n$ orthogonal matrix, which means that
$\bm Q^\prime\bm Q=\bm Q\bm Q^\prime=\bm I_n$, and the $n\times p$
matrix $\bm R$ is zero below the main diagonal.  The $n\times p$
matrix $\bm Q_1$ is the first $p$ columns of $\bm Q$ and the $p\times
p$ upper triangular matrix $\bm R_1$ is the top $p$ rows of $\bm R$.
There are three \pkg{Eigen} classes for the QR decomposition:
\code{HouseholderQR} provides the basic QR decomposition using
Householder transformations, \code{ColPivHouseholderQR} incorporates
column pivots and \code{FullPivHouseholderQR} incorporates both row
and column pivots.

Listing~\ref{QRLS} shows a least squares solution using the unpivoted
QR decomposition.
\begin{lstlisting}[frame=tb,float,caption={QRLSCpp: Least squares using the unpivoted QR decomposition},label=QRLS]
using Eigen::HouseholderQR;

const HouseholderQR<MatrixXd> QR(X);
const VectorXd           betahat(QR.solve(y));
const VectorXd            fitted(X * betahat);
const int                     df(n - p);
const VectorXd                se(QR.matrixQR().topRows(p).triangularView<Upper>().
                                 solve(MatrixXd::Identity(p,p)).rowwise().norm());
\end{lstlisting}
The calculations in Listing~\ref{QRLS} are quite similar to those in
Listing~\ref{lltLS}.  In fact, if we had extracted the upper
triangular factor (the \code{matrixU()} method) from the \code{LLT}
object in Listing~\ref{lltLS}, the rest of the code would be nearly
identical.

\subsection{Handling the rank-deficient case}
\label{sec:rankdeficient}

One important consideration when determining least squares solutions
is whether $\rank(\bm X)$ is $p$, a situation we describe by saying
that $\bm X$ has ``full column rank''.   When $\bm X$ does not have
full column rank we say it is ``rank deficient''.  

Although the theoretical rank of a matrix is well-defined, its 
evaluation in practice is not.  At best we can compute an effective
rank according to some tolerance.  We refer to decompositions that
allow us to estimate the rank of the matrix in this way as
``rank-revealing''. 

Because the \code{model.matrix} function in \proglang{R} does a
considerable amount of symbolic analysis behind the scenes, we usually
end up with full-rank model matrices.  The common cases of
rank-deficiency, such as incorporating both a constant term and a full
set of indicators columns for the levels of a factor, are eliminated.
Other, more subtle, situations will not be detected at this stage,
however.  A simple example occurs when there is a ``missing cell'' in a
two-way layout and the interaction of the two factors is included in
the model.

<<missingcell>>=
dd <- data.frame(f1 = gl(4, 6, labels = LETTERS[1:4]),
                 f2 = gl(3, 2, labels = letters[1:3]))[-(7:8), ]
xtabs(~ f2 + f1, dd)                    # one missing cell
mm <- model.matrix(~ f1 * f2, dd)
kappa(mm)         # large condition number, indicating rank deficiency
rcond(mm)         # alternative evaluation, the reciprocal condition number
(c(rank=qr(mm)$rank, p=ncol(mm))) # rank as computed in R's qr function
set.seed(1)
dd$y <- mm %*% seq_len(ncol(mm)) + rnorm(nrow(mm), sd = 0.1)
                         # lm detects the rank deficiency
fm1 <- lm(y ~ f1 * f2, dd)
writeLines(capture.output(print(summary(fm1), signif.stars=FALSE))[9:22])
@ 

The \code{lm} function for fitting linear models in \proglang{R} uses
a rank-revealing form of the QR decomposition.  When the model matrix
is determined to be rank deficient, according to the threshold used in
\proglang{R}'s QR decomposition, the model matrix is reduced to
$\rank{(\bm X)}$ columns by pivoting selected columns (those that are
apparently linearly dependent on columns to their left) to the right
hand side of the matrix.  A solution for this reduced model matrix is
determined and the coefficients and standard errors for the redundant
columns are flagged as missing.

An alternative approach is to evaluate the ``pseudo-inverse'' of $\bm
X$ from the singular value decomposition (SVD) of $\bm X$ or the
eigendecomposition of $\bm X^\prime\bm X$.  The SVD is of the form
\begin{displaymath}
  \bm X=\bm U\bm D\bm V^\prime=\bm U_1\bm D_1\bm V^\prime
\end{displaymath}
where $\bm U$ is an orthogonal $n\times n$ matrix and $\bm U_1$ is its
leftmost $p$ columns, $\bm D$ is $n\times p$ and zero off the main
diagonal so that $\bm D_1$ is a $p\times p$ diagonal matrix with
non-decreasing non-negative diagonal elements, and $\bm V$ is a $p\times
p$ orthogonal matrix.  The pseudo-inverse of $\bm D_1$, written $\bm
D_1^+$ is a $p\times p$ diagonal matrix whose first $r=\rank(\bm X)$
diagonal elements are the inverses of the corresponding diagonal
elements of $\bm D_1$ and whose last $p-r$ diagonal elements are zero.

The tolerance for determining if an element of the diagonal of $\bm D$
is considered to be (effectively) zero is a multiple of the largest
singular value (i.e. the $(1,1)$ element of $\bm D$).

In Listing~\ref{Dplus} we define a utility function, \code{Dplus}, to
return the pseudo-inverse as a diagonal matrix, given the singular
values (the diagonal of $\bm D$) and the apparent rank.  To be able to
use this function with the eigendecomposition where the eigenvalues
are in increasing order we include a Boolean argument \code{rev}
indicating whether the order is reversed.

\begin{lstlisting}[frame=tb,float,caption={DplusCpp: Create the
    diagonal matrix $\bm D^+$ from the array of singular values $\bm d$},label=Dplus]
using Eigen::DiagonalMatrix;
using Eigen::Dynamic;

inline DiagonalMatrix<double, Dynamic> Dplus(const ArrayXd& D,
                                             int r, bool rev=false) {
    VectorXd   Di(VectorXd::Constant(D.size(), 0.));
    if (rev) Di.tail(r)  = D.tail(r).inverse();
    else Di.head(r)      = D.head(r).inverse();
    return DiagonalMatrix<double, Dynamic>(Di);
}
\end{lstlisting}

\subsection{Least squares using the SVD}
\label{sec:SVDls}

With these definitions the code for least squares using the singular
value decomposition can be written as in Listing~\ref{SVDLS}.
\begin{lstlisting}[frame=tb,float,caption={SVDLSCpp: Least squares using the SVD},label=SVDLS]
using Eigen::JacobiSVD;

const JacobiSVD<MatrixXd> UDV(X.jacobiSvd(Eigen::ComputeThinU|Eigen::ComputeThinV));
const ArrayXd               D(UDV.singularValues());
const int                   r((D > D[0] * threshold()).count());
const MatrixXd            VDp(UDV.matrixV() * Dplus(D, r));
const VectorXd        betahat(VDp * UDV.matrixU().adjoint() * y);
const int                  df(n - r);
const VectorXd             se(s * VDp.rowwise().norm());
\end{lstlisting}
In the rank-deficient case this code will produce a complete set of
coefficients and their standard errors.  It is up to the user to note
that the rank is less than $p$, the number of columns in $\bm X$, and
hence that the estimated coefficients are just one of an infinite
number of coefficient vectors that produce the same fitted values.  It
happens that this solution is the minimum norm solution.

The standard errors of the coefficient estimates in the rank-deficient
case must be interpreted carefully.  The solution with one or more missing
coefficients, as returned by the \code{lm.fit} function in
\proglang{R} and the column-pivoted QR decomposition described in
Section~\ref{sec:colPivQR} does not provide standard errors for the
missing coefficients.  That is, both the coefficient and its standard
error are returned as \code{NA} because the least squares solution is
performed on a reduced model matrix.  It is also true that the
solution returned by the SVD method is with respect to a reduced model
matrix but the $p$ coefficient estimates and their $p$ standard errors
don't show this.  They are, in fact, linear combinations of a set of
$r$ coefficient estimates and their standard errors.

\subsection{Least squares using the eigendecomposition}
\label{sec:eigendecomp}

The eigendecomposition of $\bm X^\prime\bm X$ is defined as
\begin{displaymath}
  \bm X^\prime\bm X=\bm V\bm\Lambda\bm V^\prime
\end{displaymath}
where $\bm V$, the matrix of eigenvectors, is a $p\times p$ orthogonal
matrix and $\bm\Lambda$ is a $p\times p$ diagonal matrix with
non-increasing, non-negative diagonal elements, called the eigenvalues
of $\bm X^\prime\bm X$.  When the eigenvalues are distinct this $\bm
V$ is the same as that in the SVD.  Also the eigenvalues of $\bm
X^\prime\bm X$ are the squares of the singular values of $\bm X$.

With these definitions we can adapt much of the code from the SVD
method for the eigendecomposition, as shown in Listing~\ref{SymmEigLS}.
\begin{lstlisting}[frame=tb,float,caption={SymmEigLSCpp: Least squares using the eigendecomposition},label=SymmEigLS]
using Eigen::SelfAdjointEigenSolver;

const SelfAdjointEigenSolver<MatrixXd>
    VLV(MatrixXd(p, p).setZero().selfadjointView<Lower>.rankUpdate(X.adjoint()));
const ArrayXd               D(eig.eigenvalues());
const int                   r((D > D[p - 1] * threshold()).count());
const MatrixXd            VDp(VLV.eigenvectors() * Dplus(D.sqrt(), r, true));
const VectorXd        betahat(VDp * VDp.adjoint() * X.adjoint() * y);
const VectorXd             se(s * VDp.rowwise().norm());
\end{lstlisting}

\subsection{Least squares using the column-pivoted QR decomposition}
\label{sec:colPivQR}

The column-pivoted QR decomposition provides results similar to those
from \proglang{R} in both the full-rank and the rank-deficient cases.
The decomposition is of the form
\begin{displaymath}
  \bm X\bm P=\bm Q\bm R=\bm Q_1\bm R_1
\end{displaymath}
where, as before, $\bm Q$ is $n\times n$ and orthogonal and $\bm R$ is
$n\times p$ and upper triangular.  The $p\times p$ matrix $\bm P$ is a
permutation matrix.  That is, its columns are a permutation of the
columns of $\bm I_p$.  It serves to reorder the columns of $\bm X$ so
that the diagonal elements of $\bm R$ are non-increasing in magnitude.

An instance of the class \code{Eigen::ColPivHouseholderQR} has a
\code{rank()} method returning the computational rank of the matrix.
When $\bm X$ is of full rank we can use essentially the same code as
in the unpivoted decomposition except that we must reorder the
standard errors.  When $\bm X$ is rank-deficient we evaluate the
coefficients and standard errors for the leading $r$ columns of $\bm
X\bm P$ only.

In the rank-deficient case the straightforward calculation of the
fitted values, as $\bm X\widehat{\bm\beta}$, cannot be used.  We 
could do some complicated rearrangement of the columns of X and the
coefficient estimates but it is conceptually (and computationally)
easier to employ the relationship
\begin{displaymath}
  \widehat{\bm y} = \bm Q_1\bm Q_1^\prime\bm y=\bm Q
  \begin{bmatrix}
    \bm I_r & \bm 0\\
    \bm 0   & \bm 0
  \end{bmatrix}
  \bm Q^\prime\bm y
\end{displaymath}
The vector $\bm Q^\prime\bm y$ is called the ``effects'' vector in \proglang{R}.

\begin{lstlisting}[frame=tb,float,caption={ColPivQRLSCpp: Least squares using the pivoted QR decomposition},label=ColPivQRLS]
using Eigen::ColPivHouseholderQR;
typedef ColPivHouseholderQR<MatrixXd>::PermutationType  Permutation;

const ColPivHouseholderQR<MatrixXd> PQR(X);
const Permutation                   Pmat(PQR.colsPermutation());
const int                              r(PQR.rank());
VectorXd                       betahat, fitted, se;
if (r == X.cols()) {	// full rank case
    betahat  = PQR.solve(y);
    fitted   = X * betahat;
    se       = Pmat * PQR.matrixQR().topRows(p).triangularView<Upper>().
	       solve(MatrixXd::Identity(p, p)).rowwise().norm();
} else {
    MatrixXd                      Rinv(PQR.matrixQR().topLeftCorner(r, r).
				       triangularView<Upper>().
				       solve(MatrixXd::Identity(r, r)));
    VectorXd                   effects(PQR.householderQ().adjoint() * y);
    betahat.head(r)                    = Rinv * effects.head(r);
    betahat                            = Pmat * betahat;
			// create fitted values from effects
			// (cannot use X * betahat when X is rank-deficient)
    effects.tail(X.rows() - r).setZero();
    fitted                             = PQR.householderQ() * effects;
    se.head(r)                         = Rinv.rowwise().norm();
    se                                 = Pmat * se;
}
\end{lstlisting}

Just to check that the code in Listing~\ref{ColPivQRLS} does indeed provide the desired answer
<<rankdeficientPQR>>=
print(summary(fmPQR <- fastLm(y ~ f1 * f2, dd)), signif.stars=FALSE)
all.equal(coef(fm1), coef(fmPQR))
all.equal(unname(fitted(fm1)), fitted(fmPQR))
all.equal(unname(residuals(fm1)), residuals(fmPQR))
@ 

The rank-revealing SVD method produces the same fitted
values but not the same coefficients.
<<rankdeficientSVD>>=
print(summary(fmSVD <- fastLm(y ~ f1 * f2, dd, method=4L)), signif.stars=FALSE)
all.equal(coef(fm1), coef(fmSVD))
all.equal(unname(fitted(fm1)), fitted(fmSVD))
all.equal(unname(residuals(fm1)), residuals(fmSVD))
@ 
The coefficients from the symmetric eigendecomposition method are the same as those from the SVD
<<rankdeficientVLV>>=
print(summary(fmVLV <- fastLm(y ~ f1 * f2, dd, method=5L)), signif.stars=FALSE)
all.equal(coef(fmSVD), coef(fmVLV))
all.equal(unname(fitted(fm1)), fitted(fmSVD))
all.equal(unname(residuals(fm1)), residuals(fmSVD))
@ 

\subsection{Comparative speed}

In the \pkg{RcppEigen} package the \proglang{R} function to fit linear
models using the methods described above is called \code{fastLm}.  The
natural question to ask is, ``Is it indeed fast to use these methods
based on \pkg{Eigen}?''.  We have provided benchmarking code for these
methods, \proglang{R}'s \code{lm.fit} function and the \code{fastLm}
implementations in the \pkg{RcppArmadillo} \citep{CRAN:RcppArmadillo}
and \pkg{RcppGSL} \citep{CRAN:RcppGSL} packages, if they are
installed.  The benchmark code, which uses the \pkg{rbenchmark}
\citep{CRAN:rbenchmark} package, is in a file named
\code{lmBenchmark.R} in the \code{examples} subdirectory of the
installed \pkg{RcppEigen} package.


It can be run as 
<<benchmark,eval=FALSE>>=
source(system.file("examples", "lmBenchmark.R", package="RcppEigen"))
@
Results will vary according to the speed of the processor, the
number of cores and the implementation of the BLAS (Basic Linear
Algebra Subroutines) used.  (\pkg{Eigen} methods do not use the BLAS
but the other methods do.)

Results obtained on a desktop computer, circa 2010, are shown in
Table~\ref{tab:lmRes}
\begin{table}[tb]
  \caption{\code{lmBenchmark} results on a desktop computer for the
    default size, $100,000\times 40$, full-rank model matrix running
    20 repetitions for each method.  Times (Elapsed, User and Sys) are
    in seconds.  The BLAS in use is a single-threaded version of Atlas
    (Automatically Tuned Linear Algebra System).}
  \label{tab:lmRes}
  \centering
  \begin{tabular}{r r r r r}
    \hline
    \multicolumn{1}{c}{Method} & \multicolumn{1}{c}{Relative} &
    \multicolumn{1}{c}{Elapsed} & \multicolumn{1}{c}{User} &
    \multicolumn{1}{c}{Sys}\\
    \hline
     LLt &   1.000000 &   1.227 &     1.228 &    0.000 \\
    LDLt &   1.037490 &   1.273 &     1.272 &    0.000 \\
 SymmEig &   2.895681 &   3.553 &     2.972 &    0.572 \\
      QR &   7.828036 &   9.605 &     8.968 &    0.620 \\
   PivQR &   7.953545 &   9.759 &     9.120 &    0.624 \\
    arma &   8.383048 &  10.286 &    10.277 &    0.000 \\
  lm.fit &  13.782396 &  16.911 &    15.521 &    1.368 \\
     SVD &  54.829666 &  67.276 &    66.321 &    0.912 \\
     GSL & 157.531377 & 193.291 &   192.568 &    0.640 \\
     \hline
  \end{tabular}
\end{table}

These results indicate that methods based on forming and decomposing
$\bm X^\prime\bm X$, (i.e. LDLt, LLt and SymmEig) are considerably
faster than the others.  The SymmEig method, using a rank-revealing
decomposition, would be preferred, although the LDLt method could
probably be modified to be rank-revealing.  Do bear in mind that the
dimensions of the problem will influence the comparative results.
Because there are 100,000 rows in $\bm X$, methods that decompose the
whole $\bm X$ matrix (all the methods except those named above) will
be at a disadvantage.

The pivoted QR method is 1.6 times faster than R's \Sexpr{link("lm.fit")} on
this test and provides nearly the same information as \Sexpr{link("lm.fit")}.
Methods based on the singular value decomposition (SVD and GSL) are
much slower but, as mentioned above, this is caused in part by $\bm X$
having many more rows than columns.  The GSL method from the GNU
Scientific Library uses an older algorithm for the SVD and is clearly
out of contention.

An SVD method using the Lapack SVD subroutine, \code{dgesv}, may be
faster than the native \pkg{Eigen} implementation of the SVD, which is
not a particularly fast method.

\section{Delayed evaluation}
\label{sec:delayed}

A form of delayed evaluation is used in \pkg{Eigen}.  That is, many
operators and methods do not force the evaluation of the object but
instead return an ``expression object'' that is evaluated when
needed.  As an example, even though we write the $\bm X^\prime\bm X$
evaluation using \code{.rankUpdate(X.adjoint())} the
\code{X.adjoint()} part is not evaluated immediately.  The
\code{rankUpdate} method detects that it has been passed a matrix
that is to be used in its transposed form and evaluates the update
by taking inner products of columns of $\bm X$ instead of rows of $\bm
X^\prime$.

Occasionally the method for \code{Rcpp::wrap} will not force an
evaluation when it should.  This is at least what Bill Venables calls
an ``infelicity'' in the code, if not an outright bug.  In the code
for the transpose of an integer matrix shown in Listing~\ref{trans} we
assigned the transpose as a \code{MatrixXi} before returning it with
\code{wrap}.  The assignment forces the evaluation.  If we skip this
step, as in Listing~\ref{badtrans} we get an answer with the correct
shape but incorrect contents.

<<badtransCpp,echo=FALSE>>=
badtransCpp <- '
using Eigen::Map;
using Eigen::MatrixXi;
const Map<MatrixXi>  A(as<Map<MatrixXi> >(AA));
return wrap(A.transpose());
'
@
\begin{lstlisting}[frame=tb,float,caption={badtransCpp: Transpose producing incorrect results},label=badtrans]
<<badtransCppLst,results=tex,echo=FALSE>>=
cat(badtransCpp, "\n")
@ 
\end{lstlisting}
<<ftrans2>>=
Ai <- matrix(1:6, ncol=2L)
ftrans2 <- cxxfunction(signature(AA = "matrix"), badtransCpp, "RcppEigen")
(At <- ftrans2(Ai))
all.equal(At, t(Ai))
@

Another recommended practice is to assign objects before wrapping them
for return to \proglang{R}.

\section{Sparse matrices}
\label{sec:sparse}

\pkg{Eigen} provides sparse matrix classes.  An \proglang{R} object of
class \Sexpr{linkS4class("dgCMatrix")} (from the \pkg{Matrix}
\citep{CRAN:Matrix} package) can be mapped as in Listing~\ref{sparseProd}.
<<echo=FALSE>>=
sparseProdCpp <- '
using Eigen::Map;
using Eigen::MappedSparseMatrix;
using Eigen::SparseMatrix;
using Eigen::VectorXd;

const MappedSparseMatrix<double>  A(as<MappedSparseMatrix<double> >(AA));
const Map<VectorXd>               y(as<Map<VectorXd> >(yy));
const SparseMatrix<double>       At(A.adjoint());
return List::create(_["At"]  = At,
                    _["Aty"] = At * y);
'
@
\begin{lstlisting}[frame=tb,float,
  caption={sparseProdCpp: Transpose and product with sparse matrices},label=sparseProd]
<<sparseProdCppLst,results=tex,echo=FALSE>>=
cat(sparseProdCpp, "\n")
@ 
\end{lstlisting}
<<>>=
sparse1 <- cxxfunction(signature(AA = "dgCMatrix", yy = "numeric"),
                       sparseProdCpp, "RcppEigen")
data(KNex, package="Matrix")
rr <- sparse1(KNex$mm, KNex$y)
stopifnot(all.equal(rr$At, t(KNex$mm)),
          all.equal(rr$Aty, as.vector(crossprod(KNex$mm, KNex$y))))
@

A sparse Cholesky decomposition is provided in \pkg{Eigen} as the
\code{SimplicialCholesky} class. There are also linkages to the
\pkg{CHOLMOD} code from the \pkg{Matrix} package.  At present, both of
these are regarded as experimental.

\bibliographystyle{plainnat}
\bibliography{Rcpp}

\end{document}

\subsection{Transpose and adjoint of a complex matrix}
\label{sec:complex}
%% this example exposes a problem in converting Rcomplex* to std::complex<double>
<<echo=FALSE,eval=FALSE>>=
code <- '
using Eigen::Map;
using Eigen::MatrixXcd;
  // Map the complex matrix A_ from R
const Map<MatrixXcd>    A(as<Map<MatrixXcd> >(A_));
return List::create(_["transpose"] = A.transpose(),
                    _["adjoint"]   = A.adjoint());
'
writeLines( code, "code.cpp" )
@
<<echo=FALSE,eval=FALSE,results=tex>>=
ex_highlight( "code.cpp" )
@

<<eval=FALSE>>=
fadj <-
    cxxfunction(signature( A_ = "matrix"),
                paste(readLines( "code.cpp" ), collapse = "\n"),
                plugin = "RcppEigen")
A <- matrix(1:6 + 1i*(6:1), nc=2)
fadj(A)
@
