[Rcpp-commits] r4224 - in pkg/RcppEigen/vignettes: . jss835
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Jan 26 14:35:51 CET 2013
Author: edd
Date: 2013-01-26 14:35:50 +0100 (Sat, 26 Jan 2013)
New Revision: 4224
Added:
pkg/RcppEigen/vignettes/jss835/
pkg/RcppEigen/vignettes/jss835/jss835.zip
pkg/RcppEigen/vignettes/jss835/jss835zeileis.pdf
pkg/RcppEigen/vignettes/jss835/jss835zeileis.tex
Log:
housekeeping
Added: pkg/RcppEigen/vignettes/jss835/jss835.zip
===================================================================
(Binary files differ)
Property changes on: pkg/RcppEigen/vignettes/jss835/jss835.zip
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: pkg/RcppEigen/vignettes/jss835/jss835zeileis.pdf
===================================================================
(Binary files differ)
Property changes on: pkg/RcppEigen/vignettes/jss835/jss835zeileis.pdf
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: pkg/RcppEigen/vignettes/jss835/jss835zeileis.tex
===================================================================
--- pkg/RcppEigen/vignettes/jss835/jss835zeileis.tex (rev 0)
+++ pkg/RcppEigen/vignettes/jss835/jss835zeileis.tex 2013-01-26 13:35:50 UTC (rev 4224)
@@ -0,0 +1,1518 @@
+\documentclass[article]{jss}
+\usepackage{booktabs,bm,amsmath,thumbpdf}
+
+\newcommand{\argmin}{\operatorname{argmin}\displaylimits}
+\newcommand{\rank}{\operatorname{rank}}
+\newcommand{\hlstd}[1]{\textcolor[rgb]{0,0,0}{#1}}
+\newcommand{\hlnum}[1]{\textcolor[rgb]{0,0,0}{#1}}
+\newcommand{\hlopt}[1]{\textcolor[rgb]{0,0,0}{#1}}
+\newcommand{\hlesc}[1]{\textcolor[rgb]{0.74,0.55,0.55}{#1}}
+\newcommand{\hlstr}[1]{\textcolor[rgb]{0.90,0.15,0.15}{#1}}
+\newcommand{\hldstr}[1]{\textcolor[rgb]{0.74,0.55,0.55}{#1}}
+\newcommand{\hlslc}[1]{\textcolor[rgb]{0.67,0.13,0.13}{\it{#1}}}
+\newcommand{\hlcom}[1]{\textcolor[rgb]{0.67,0.13,0.13}{\it{#1}}}
+\newcommand{\hldir}[1]{\textcolor[rgb]{0,0,0}{#1}}
+\newcommand{\hlsym}[1]{\textcolor[rgb]{0,0,0}{#1}}
+\newcommand{\hlline}[1]{\textcolor[rgb]{0.33,0.33,0.33}{#1}}
+\newcommand{\hlkwa}[1]{\textcolor[rgb]{0.61,0.13,0.93}{\bf{#1}}}
+\newcommand{\hlkwb}[1]{\textcolor[rgb]{0.13,0.54,0.13}{#1}}
+\newcommand{\hlkwc}[1]{\textcolor[rgb]{0,0,1}{#1}}
+\newcommand{\hlkwd}[1]{\textcolor[rgb]{0,0,0}{#1}}
+\definecolor{bgcolor}{rgb}{1,1,1}
+
+
+\author{Douglas Bates\\University of Wisconsin-Madison \And Dirk Eddelbuettel\\Debian Project}
+\Plainauthor{Douglas Bates, Dirk Eddelbuettel}
+
+\title{Fast and Elegant Numerical Linear Algebra Using the \pkg{RcppEigen} Package}
+\Plaintitle{Fast and Elegant Numerical Linear Algebra Using the RcppEigen Package}
+\Shorttitle{Fast and Elegant Numerical Linear Algebra with \pkg{RcppEigen}}
+
+\Abstract{
+ The \pkg{RcppEigen} package provides access from \proglang{R}
+ \citep{R:Main} to the \pkg{Eigen} \citep*{Eigen:Web} \proglang{C++}
+ template library for numerical linear algebra. \pkg{Rcpp}
+ \citep{JSS:Rcpp,CRAN: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. Several introductory
+ examples are presented. This is followed by an in-depth discussion of various
+ available approaches for solving least-squares problems, including
+ rank-revealing methods, concluding with an empirical run-time
+ comparison. Last but not least, sparse matrix methods are discussed.
+}
+
+\Keywords{linear algebra, template programming, \proglang{R}, \proglang{C++}, \pkg{Rcpp}}
+\Plainkeywords{linear algebra, template programmig, R, C++, Rcpp}
+
+\Volume{52}
+\Issue{5}
+\Year{2013}
+\Month{January}
+\Submitdate{2011-11-14}
+\Acceptdate{2012-09-26}
+
+\Address{
+ Douglas Bates \\
+ Department of Statistics \\
+ University of Wisconsin-Madison \\
+ Madison, WI, United States of America \\
+ E-mail: \email{bates at stat.wisc.edu} \\
+ URL: \url{http://www.stat.wisc.edu/~bates/}\\
+
+ Dirk Eddelbuettel \\
+ Debian Project \\
+ River Forest, IL, United States of America\\
+ E-mail: \email{edd at debian.org}\\
+ URL: \url{http://dirk.eddelbuettel.com}\\
+}
+
+
+\begin{document}
+
+\section{Introduction}
+\label{sec:intro}
+
+Linear algebra is an essential building block of statistical
+computing. Operations such as matrix decompositions, linear program
+solvers, and eigenvalue/eigenvector computations are used in many
+estimation and analysis routines. As such, libraries supporting linear
+algebra have long been provided by statistical programmers for
+different programming languages and environments. Because it is
+object-oriented, \proglang{C++}, one of the central modern languages
+for numerical and statistical computing, is particularly effective at
+representing matrices, vectors and decompositions, and numerous class
+libraries providing linear algebra routines have been written over the
+years.
+
+As both the \proglang{C++} language and standards have evolved
+\citep{Meyers:2005:EffectiveC++,Meyers:1995:MoreEffectiveC++,cpp11},
+so have the compilers implementing the language. Relatively modern
+language constructs such as template meta-programming are particularly
+useful because they provide overloading of operations (allowing
+expressive code in the compiled language similar to what can be done
+in scripting languages) and can shift some of the computational burden
+from the run-time to the compile-time. (A more detailed discussion of
+template meta-programming is, however, beyond the scope of this
+paper). \cite{Veldhuizen:1998:Blitz} provided an early and influential
+implementation of numerical linear algebra classes that already
+demonstrated key features of this approach. Its usage was held back
+at the time by the limited availability of compilers implementing all the
+necessary features of the \proglang{C++} language.
+
+This situation has greatly improved over the last decade, and many more
+libraries have been made available. One such \proglang{C++} library is
+\pkg{Eigen} by \citet*{Eigen:Web} which started as a sub-project to
+KDE (a popular Linux desktop environment), initially focussing on fixed-size
+matrices to represent projections in a visualization application. \pkg{Eigen}
+grew from there and has over the course of about a decade produced three
+major releases with \pkg{Eigen}3 being the current major version. To
+check the minor and patch version numbers, load the \pkg{RcppEigen}
+package and call
+\begin{CodeInput}
+R> .Call("eigen_version", FALSE)
+\end{CodeInput}
+\begin{CodeOutput}
+major minor patch
+ 3 1 1
+\end{CodeOutput}
+\pkg{Eigen} is of interest as the \proglang{R} system for statistical
+computation and graphics \citep{R:Main} is itself easily extensible. This is
+particular true via the \proglang{C} language that most of \proglang{R}'s
+compiled core parts are written in, but also for the \proglang{C++} language
+which can interface with \proglang{C}-based systems rather easily. The manual
+``Writing \proglang{R} Extensions'' \citep{R:Extensions} is the basic reference for
+extending \proglang{R} with either \proglang{C} or \proglang{C++}.
+
+The \pkg{Rcpp} package by \citet{JSS:Rcpp,CRAN:Rcpp} facilitates extending
+\proglang{R} with \proglang{C++} code by providing seamless object mapping
+between both languages.
+%
+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 \code{Rcpp::as} (for
+ conversion of objects from \proglang{R} to \proglang{C++}) and
+ \code{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, along with implementations of
+\code{Rcpp::as} and \code{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. The next section introduces some
+of these classes and shows how to interface to them from \proglang{R}.
+
+\section[Eigen classes]{\pkg{Eigen} classes}
+\label{sec:eclasses}
+
+\pkg{Eigen} 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 \pkg{Eigen} is often more like
+\proglang{R} code, working on the ``whole object'', than like 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
+\code{typedef}s for common classes that correspond to \proglang{R} matrices and
+vectors, as shown in Table~\ref{tab:REigen}, and this paper will use these
+\code{typedef}s.
+
+\begin{table}[t!]
+ \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} \phantom{XXX} & \code{SparseMatrix<double>}\\
+ \hline
+ \end{tabular}
+ \caption{Correspondence between \proglang{R} matrix and vector types and classes in the \pkg{Eigen} namespace.
+ \label{tab:REigen}}
+\end{table}
+
+Here, \code{Vector} and \code{Matrix} describe the dimension of the
+object. The \code{X} signals that these are dynamically-sized objects (as opposed
+to fixed-size matrices such as $3 \times 3$ matrices also available in
+\pkg{Eigen}). Lastly, the trailing characters \code{i}, \code{d} and
+\code{cd} denote, respectively, storage types \code{integer}, \code{double} and
+\code{complex double}.
+
+The \proglang{C++} classes shown in Table~\ref{tab:REigen} are in the
+\pkg{Eigen} namespace, which means that they must be written as
+\code{Eigen::MatrixXd}. However, if one prefaces the use of these class
+names with a declaration like
+\begin{quote}
+ \noindent
+ \ttfamily
+ \hlstd{}\hlkwa{using\ }\hlstd{Eigen}\hlopt{::}\hlstd{MatrixXd}\hlopt{;}\hlstd{}\hspace*{\fill}\\
+ \mbox{}
+ \normalfont
+ \normalsize
+\end{quote}
+
+\vspace*{-0.4cm}
+
+then one can use these names without the namespace qualifier.
+
+\subsection[Mapped matrices in Eigen]{Mapped matrices in \pkg{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 \pkg{Eigen} class. For
+dense matrices one can use the \pkg{Eigen} templated class \code{Map}, and for
+sparse matrices one can deploy the \pkg{Eigen} templated class \code{MappedSparseMatrix}.
+
+One 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 {\ttfamily\hlkwb{const}\normalfont}.
+
+\subsection[Arrays in Eigen]{Arrays in \pkg{Eigen}}
+\label{sec:arrays}
+
+For matrix and vector classes \pkg{Eigen} overloads the \code{`*'}
+operator as matrix multiplication. Occasionally component-wise
+operations instead of matrix operations are desired, for which the
+\code{Array} templated classes are used in \pkg{Eigen}. Switching
+back and forth between matrices or vectors and arrays is usually done
+via the \code{array()} method for matrix or vector objects and the
+\code{matrix()} method for arrays.
+
+\subsection[Structured matrices in Eigen]{Structured matrices in \pkg{Eigen}}
+\label{sec:structured}
+
+\pkg{Eigen} provides 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 one must
+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 \pkg{Eigen} objects.
+\item Create the result.
+\item Return \code{Rcpp::wrap} applied to the result.
+\end{enumerate}
+
+An idiom for the first step is
+% using Eigen::Map;
+% using Eigen::MatrixXd;
+% using Rcpp::as;
+%
+% const Map<MatrixXd> A(as<Map<MatrixXd> >(AA));
+%\end{lstlisting}
+\begin{quote}
+ \noindent
+ \ttfamily
+ \hlstd{}\hlkwa{using\ }\hlstd{Eigen}\hlsym{::}\hlstd{Map}\hlsym{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwa{using\ }\hlstd{Eigen}\hlsym{::}\hlstd{MatrixXd}\hlsym{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwa{using\ }\hlstd{Rcpp}\hlsym{::}\hlstd{as}\hlsym{;}\hspace*{\fill}\\
+ \hlstd{}\hspace*{\fill}\\
+ \hlkwb{const\ }\hlstd{Map}\hlsym{$<$}\hlstd{MatrixXd}\hlsym{$>$}\hlstd{\ \ }\hlsym{}\hlstd{}\hlkwd{A}\hlstd{}\hlsym{(}\hlstd{as}\hlsym{$<$}\hlstd{Map}\hlsym{$<$}\hlstd{MatrixXd}\hlsym{$>$\ $>$(}\hlstd{AA}\hlsym{));}\hlstd{}\hspace*{\fill}\\
+ \mbox{}
+ \normalfont
+\end{quote}
+
+\vspace*{-0.4cm}
+
+where \code{AA} is the name of the \proglang{R} object (of type \code{SEXP} in
+\proglang{C} and \proglang{C++}) passed to the \proglang{C++} function.
+
+An alternative to the \code{using} declarations is to declare a \code{typedef} as in
+% typedef Eigen::Map<Eigen::MatrixXd> MapMatd;
+% const MapMatd A(Rcpp::as<MapMatd>(AA));
+\begin{quote}
+ \noindent
+ \ttfamily
+ \hlstd{}\hlkwc{typedef\ }\hlstd{Eigen}\hlopt{::}\hlstd{Map}\hlopt{$<$}\hlstd{Eigen}\hlopt{::}\hlstd{MatrixXd}\hlopt{$>$}\hlstd{\ \ \ }\hlopt{}\hlstd{MapMatd}\hlopt{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwb{const\ }\hlstd{MapMatd}\hlstd{\ \ \ \ \ \ \ \ \ \ }\hlstd{}\hlkwd{A}\hlstd{}\hlopt{(}\hlstd{Rcpp}\hlopt{::}\hlstd{as}\hlopt{$<$}\hlstd{MapMatd}\hlopt{$>$(}\hlstd{AA}\hlopt{));}\hlstd{}\hspace*{\fill}\\
+ \mbox{}
+ \normalfont
+ \normalsize
+\end{quote}
+
+\vspace*{-0.4cm}
+
+The \code{cxxfunction} function from the \pkg{inline} package \citep*{CRAN:inline} for
+\proglang{R} and its \pkg{RcppEigen} plugin provide a convenient method of
+developing and debugging the \proglang{C++} code. For actual production code
+one generally incorporates the \proglang{C++} source code files in a package
+and includes the line \code{LinkingTo: Rcpp, RcppEigen} in the package's
+\code{DESCRIPTION} file. The \code{RcppEigen.package.skeleton} function
+provides a quick way of generating the skeleton of a package that will use
+\pkg{RcppEigen}.
+
+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
+these declarations are omitted. It is important to remember that they are
+needed in actual \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}.
+Finally, the results from \pkg{Eigen} are compared to those from the direct
+\proglang{R} results.
+
+\subsection{Transpose of an integer matrix}
+\label{sec:transpose}
+
+The next \proglang{R} code snippet creates a simple matrix of integers
+\begin{CodeInput}
+R> (A <- matrix(1:6, ncol = 2))
+\end{CodeInput}
+\begin{CodeOutput}
+ [,1] [,2]
+[1,] 1 4
+[2,] 2 5
+[3,] 3 6
+\end{CodeOutput}
+\begin{CodeInput}
+R> str(A)
+\end{CodeInput}
+\begin{CodeOutput}
+ int [1:3, 1:2] 1 2 3 4 5 6
+\end{CodeOutput}
+and, in Figure~\ref{trans}, the \code{transpose()} method for the
+\code{Eigen::MatrixXi} class is used to return the transpose of the supplied matrix. The \proglang{R}
+matrix in the \code{SEXP} \code{AA} is first mapped to an
+\code{Eigen::MatrixXi} object, and then the matrix \code{At} is constructed
+from its transpose and returned to \proglang{R}.
+
+\begin{figure}[t!]
+ \hrule \smallskip
+ \noindent
+ \ttfamily
+ \hlstd{}\hlkwa{using\ }\hlstd{Eigen}\hlsym{::}\hlstd{Map}\hlsym{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwa{using\ }\hlstd{Eigen}\hlsym{::}\hlstd{MatrixXi}\hlsym{;}\hspace*{\fill}\\
+ \hlstd{}\hlstd{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\hlstd{}\hlslc{//\ Map\ the\ integer\ matrix\ AA\ from\ R}\hspace*{\fill}\\
+ \hlstd{}\hlkwb{const\ }\hlstd{Map}\hlsym{$<$}\hlstd{MatrixXi}\hlsym{$>$}\hlstd{\ \ }\hlsym{}\hlstd{}\hlkwd{A}\hlstd{}\hlsym{(}\hlstd{as}\hlsym{$<$}\hlstd{Map}\hlsym{$<$}\hlstd{MatrixXi}\hlsym{$>$\ $>$(}\hlstd{AA}\hlsym{));}\hspace*{\fill}\\
+ \hlstd{}\hlstd{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\hlstd{}\hlslc{//\ evaluate\ and\ return\ the\ transpose\ of\ A}\hspace*{\fill}\\
+ \hlstd{}\hlkwb{const\ }\hlstd{MatrixXi}\hlstd{\ \ \ \ \ \ }\hlstd{}\hlkwd{At}\hlstd{}\hlsym{(}\hlstd{A}\hlsym{.}\hlstd{}\hlkwd{transpose}\hlstd{}\hlsym{());}\hspace*{\fill}\\
+ \hlstd{}\hlkwa{return\ }\hlstd{}\hlkwd{wrap}\hlstd{}\hlsym{(}\hlstd{At}\hlsym{);}\hlstd{}\hspace*{\fill}
+ \normalfont
+ \hrule
+ \caption{\code{transCpp}: Transpose a matrix of integers.
+ \label{trans}}
+\end{figure}
+
+The \proglang{R} snippet below compiles and links the \proglang{C++} code
+segment. The actual work is done by the function \code{cxxfunction} from the \pkg{inline}
+package which compiles, links and loads code written in \proglang{C++} and
+supplied as a character variable. This frees the user from having to know about
+compiler and linker details and options, which makes ``exploratory
+programming'' much easier. Here the program piece to be compiled is stored
+as a character variable named \code{transCpp}, and \code{cxxfunction} creates
+an executable function which is assigned to \code{ftrans}. This new function
+is used on the matrix $\bm A$ shown above, and one can check that it works as intended
+by comparing the output to an explicit transpose of the matrix argument.
+\begin{CodeInput}
+R> ftrans <- cxxfunction(signature(AA = "matrix"), transCpp,
++ plugin = "RcppEigen")
+R> (At <- ftrans(A))
+\end{CodeInput}
+\begin{CodeOutput}
+ [,1] [,2] [,3]
+[1,] 1 2 3
+[2,] 4 5 6
+\end{CodeOutput}
+\begin{CodeInput}
+R> stopifnot(all.equal(At, t(A)))
+\end{CodeInput}
+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 the \code{adjoint()}
+method is used in what follows to create the transpose of numeric or
+integer matrices.
+
+\subsection{Products and cross-products}
+\label{sec:products}
+
+As mentioned in Section~\ref{sec:arrays}, the \code{`*'} operator is
+overloaded as matrix multiplication for the various matrix and vector
+classes in \pkg{Eigen}. The \proglang{C++} code in
+Figure~\ref{prod} produces a list containing both the product and
+cross-product (in the sense of the \proglang{R} function call
+\code{crossproduct(A, B)} evaluating $\bm A^\top\bm B$) of its two arguments
+%
+\begin{figure}[t!]
+ \hrule \smallskip
+ \noindent
+ \ttfamily
+ \hlstd{}\hlkwc{typedef\ }\hlstd{Eigen}\hlopt{::}\hlstd{Map}\hlopt{$<$}\hlstd{Eigen}\hlopt{::}\hlstd{MatrixXi}\hlopt{$>$}\hlstd{\ \ \ }\hlopt{}\hlstd{MapMati}\hlopt{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwb{const\ }\hlstd{MapMati}\hlstd{\ \ \ \ }\hlstd{}\hlkwd{B}\hlstd{}\hlopt{(}\hlstd{as}\hlopt{$<$}\hlstd{MapMati}\hlopt{$>$(}\hlstd{BB}\hlopt{));}\hspace*{\fill}\\
+ \hlstd{}\hlkwb{const\ }\hlstd{MapMati}\hlstd{\ \ \ \ }\hlstd{}\hlkwd{C}\hlstd{}\hlopt{(}\hlstd{as}\hlopt{$<$}\hlstd{MapMati}\hlopt{$>$(}\hlstd{CC}\hlopt{));}\hspace*{\fill}\\
+ \hlstd{}\hlkwa{return\ }\hlstd{List}\hlopt{::}\hlstd{}\hlkwd{create}\hlstd{}\hlopt{(}\hlstd{Named}\hlopt{{(}}\hlstd{}\hlstr{"B\ \%{*}\%\ C"}\hlstd{}\hlopt{{)}}\hlstd{\ \ \ \ \ \ \ \ \ }\hlopt{=\ }\hlstd{B\ }\hlopt{{*}\ }\hlstd{C}\hlopt{,}\hspace*{\fill}\\
+ \hlstd{}\hlstd{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\hlstd{Named}\hlopt{{(}}\hlstd{}\hlstr{"crossprod(B,\ C)"}\hlstd{}\hlopt{{)}\ =\ }\hlstd{B}\hlopt{.}\hlstd{}\hlkwd{adjoint}\hlstd{}\hlopt{()\ {*}\ }\hlstd{C}\hlopt{);}\hlstd{}\hspace*{\fill}\\
+ \mbox{}
+ \normalfont
+ \normalsize
+ \hrule
+ \caption{\code{prodCpp}: Product and cross-product of two matrices.
+ \label{prod}}
+\end{figure}
+%
+\begin{CodeInput}
+R> fprod <- cxxfunction(signature(BB = "matrix", CC = "matrix"),
++ prodCpp, "RcppEigen")
+R> B <- matrix(1:4, ncol = 2)
+R> C <- matrix(6:1, nrow = 2)
+R> str(fp <- fprod(B, C))
+\end{CodeInput}
+\begin{CodeOutput}
+List of 2
+ $ B %*% C : int [1:2, 1:3] 21 32 13 20 5 8
+ $ crossprod(B, C): int [1:2, 1:3] 16 38 10 24 4 10
+\end{CodeOutput}
+\begin{CodeInput}
+R> stopifnot(all.equal(fp[[1]], B %*% C), all.equal(fp[[2]], crossprod(B, C)))
+\end{CodeInput}
+%
+Note that the \code{create} class method for \code{Rcpp::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^\top\bm X$. One could, of
+course, calculate this product as
+\begin{verbatim}
+ t(X) %*% X
+\end{verbatim}
+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 \pkg{Eigen}, a
+\code{SelfAdjointView}---which is a dense matrix of which only one
+triangle is used, the other triangle being inferred from the
+symmetry---is created. (The characterization ``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^\top$ to the current symmetric matrix. The scalar
+multiple defaults to 1. The code in Figure~\ref{crossprod} produces
+both $\bm A^\top\bm A$ and $\bm A\bm A^\top$ from a matrix $\bm A$.
+
+\begin{figure}[t!]
+ \hrule \smallskip
+ \noindent
+ \ttfamily
+ \hlstd{}\hlkwa{using\ }\hlstd{Eigen}\hlopt{::}\hlstd{Map}\hlopt{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwa{using\ }\hlstd{Eigen}\hlopt{::}\hlstd{MatrixXi}\hlopt{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwa{using\ }\hlstd{Eigen}\hlopt{::}\hlstd{Lower}\hlopt{;}\hspace*{\fill}\\
+ \hlstd{}\hspace*{\fill}\\
+ \hlkwb{const\ }\hlstd{Map}\hlopt{$<$}\hlstd{MatrixXi}\hlopt{$>$\ }\hlstd{}\hlkwd{A}\hlstd{}\hlopt{(}\hlstd{as}\hlopt{$<$}\hlstd{Map}\hlopt{$<$}\hlstd{MatrixXi}\hlopt{$>$\ $>$(}\hlstd{AA}\hlopt{));}\hspace*{\fill}\\
+ \hlstd{}\hlkwb{const\ int}\hlstd{\ \ \ \ \ \ \ \ \ \ \ }\hlkwb{}\hlstd{}\hlkwd{m}\hlstd{}\hlopt{(}\hlstd{A}\hlopt{.}\hlstd{}\hlkwd{rows}\hlstd{}\hlopt{()),\ }\hlstd{}\hlkwd{n}\hlstd{}\hlopt{(}\hlstd{A}\hlopt{.}\hlstd{}\hlkwd{cols}\hlstd{}\hlopt{());}\hspace*{\fill}\\
+ \hlstd{MatrixXi}\hlstd{\ \ \ \ \ \ \ \ \ \ }\hlstd{}\hlkwd{AtA}\hlstd{}\hlopt{(}\hlstd{}\hlkwd{MatrixXi}\hlstd{}\hlopt{(}\hlstd{n}\hlopt{,\ }\hlstd{n}\hlopt{).}\hlstd{}\hlkwd{setZero}\hlstd{}\hlopt{().}\hspace*{\fill}\\
+ \hlstd{}\hlstd{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\hlstd{selfadjointView}\hlopt{$<$}\hlstd{Lower}\hlopt{$>$().}\hlstd{}\hlkwd{rankUpdate}\hlstd{}\hlopt{(}\hlstd{A}\hlopt{.}\hlstd{}\hlkwd{adjoint}\hlstd{}\hlopt{()));}\hspace*{\fill}\\
+ \hlstd{MatrixXi}\hlstd{\ \ \ \ \ \ \ \ \ \ }\hlstd{}\hlkwd{AAt}\hlstd{}\hlopt{(}\hlstd{}\hlkwd{MatrixXi}\hlstd{}\hlopt{(}\hlstd{m}\hlopt{,\ }\hlstd{m}\hlopt{).}\hlstd{}\hlkwd{setZero}\hlstd{}\hlopt{().}\hspace*{\fill}\\
+ \hlstd{}\hlstd{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\hlstd{selfadjointView}\hlopt{$<$}\hlstd{Lower}\hlopt{$>$().}\hlstd{}\hlkwd{rankUpdate}\hlstd{}\hlopt{(}\hlstd{A}\hlopt{));}\hspace*{\fill}\\
+ \hlstd{}\hspace*{\fill}\\
+ \hlkwa{return\ }\hlstd{List}\hlopt{::}\hlstd{}\hlkwd{create}\hlstd{}\hlopt{(}\hlstd{Named}\hlopt{{(}}\hlstd{}\hlstr{"crossprod(A)"}\hlstd{}\hlopt{{)}}\hlstd{\ \ }\hlopt{=\ }\hlstd{AtA}\hlopt{,}\hspace*{\fill}\\
+ \hlstd{}\hlstd{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\hlstd{Named}\hlopt{{(}}\hlstd{}\hlstr{"tcrossprod(A)"}\hlstd{}\hlopt{{)}\ =\ }\hlstd{AAt}\hlopt{);}\hlstd{}\hspace*{\fill}
+ \normalfont
+ \normalsize
+ \hrule
+ \caption{\code{crossprodCpp}: Cross-product and transposed cross-product of a single matrix.
+ \label{crossprod}}
+\end{figure}
+\begin{CodeInput}
+R> fcprd <- cxxfunction(signature(AA = "matrix"), crossprodCpp, "RcppEigen")
+R> str(crp <- fcprd(A))
+\end{CodeInput}
+\begin{CodeOutput}
+List of 2
+ $ crossprod(A) : int [1:2, 1:2] 14 32 32 77
+ $ tcrossprod(A): int [1:3, 1:3] 17 22 27 22 29 36 27 36 45
+\end{CodeOutput}
+\begin{CodeInput}
+R> stopifnot(all.equal(crp[[1]], crossprod(A)),
++ all.equal(crp[[2]], tcrossprod(A)))
+\end{CodeInput}
+To some, the expressions in Figure~\ref{crossprod} to construct
+\code{AtA} and \code{AAt} 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$), ensuring that all its
+elements are zero, regarding it as a self-adjoint (i.e., symmetric) matrix
+using the elements in the lower triangle, then adding $\bm A\bm A^\top$
+to it and converting back to a general matrix form (i.e., the strict lower
+triangle is copied into the strict upper triangle).
+
+In more detail:
+\begin{enumerate}
+\item \code{MatrixXi(n, n)} creates an $n\times n$ integer matrix with
+ arbitrary contents
+\item \code{.setZero()} zeros out the contents of the matrix
+\item \code{.selfAdjointView<Lower>()} causes what follows to treat
+ the matrix as a symmetric matrix in which only the lower triangle is
+ used, the strict upper triangle being inferred by symmetry
+\item \code{.rankUpdate(A)} forms the sum $\bm B+\bm A\bm A^\top$
+ where $\bm B$ is the symmetric matrix of zeros created in the
+ previous steps.
+\end{enumerate}
+The assignment of this symmetric matrix to the (general)
+\code{MatrixXi} object \code{AAt} causes the result to be symmetrized
+during the assignment.
+
+For these products one could define the symmetric matrix from either
+the lower triangle or the upper triangle as the result will be
+symmetrized before it is returned.
+
+To cut down on repetition of \code{using} statements we gather them in
+a character variable, \code{incl}, that will be given as the \code{includes} argument
+in the calls to \code{cxxfunction}. We also define a utility
+function, \code{AtA}, that returns the crossproduct matrix as shown in Figure~\ref{fig:incl}
+
+\begin{figure}[t!]
+ \hrule \smallskip
+ \noindent
+ \ttfamily
+ \hlstd{}\hlkwa{using}\hlstd{\ \ \ }\hlkwa{}\hlstd{Eigen}\hlopt{::}\hlstd{LLT}\hlopt{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwa{using}\hlstd{\ \ \ }\hlkwa{}\hlstd{Eigen}\hlopt{::}\hlstd{Lower}\hlopt{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwa{using}\hlstd{\ \ \ }\hlkwa{}\hlstd{Eigen}\hlopt{::}\hlstd{Map}\hlopt{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwa{using}\hlstd{\ \ \ }\hlkwa{}\hlstd{Eigen}\hlopt{::}\hlstd{MatrixXd}\hlopt{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwa{using}\hlstd{\ \ \ }\hlkwa{}\hlstd{Eigen}\hlopt{::}\hlstd{MatrixXi}\hlopt{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwa{using}\hlstd{\ \ \ }\hlkwa{}\hlstd{Eigen}\hlopt{::}\hlstd{Upper}\hlopt{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwa{using}\hlstd{\ \ \ }\hlkwa{}\hlstd{Eigen}\hlopt{::}\hlstd{VectorXd}\hlopt{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwc{typedef\ }\hlstd{Map}\hlopt{$<$}\hlstd{MatrixXd}\hlopt{$>$}\hlstd{\ \ }\hlopt{}\hlstd{MapMatd}\hlopt{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwc{typedef\ }\hlstd{Map}\hlopt{$<$}\hlstd{MatrixXi}\hlopt{$>$}\hlstd{\ \ }\hlopt{}\hlstd{MapMati}\hlopt{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwc{typedef\ }\hlstd{Map}\hlopt{$<$}\hlstd{VectorXd}\hlopt{$>$}\hlstd{\ \ }\hlopt{}\hlstd{MapVecd}\hlopt{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwc{inline\ }\hlstd{MatrixXd\ }\hlkwd{AtA}\hlstd{}\hlopt{(}\hlstd{}\hlkwb{const\ }\hlstd{MapMatd}\hlopt{\&\ }\hlstd{A}\hlopt{)\ \{}\hspace*{\fill}\\
+ \hlstd{}\hlstd{\ \ \ \ }\hlstd{}\hlkwb{int}\hlstd{\ \ \ \ }\hlkwb{}\hlstd{}\hlkwd{n}\hlstd{}\hlopt{(}\hlstd{A}\hlopt{.}\hlstd{}\hlkwd{cols}\hlstd{}\hlopt{());}\hspace*{\fill}\\
+ \hlstd{}\hlstd{\ \ \ \ }\hlstd{}\hlkwa{return}\hlstd{\ \ \ }\hlkwa{}\hlstd{}\hlkwd{MatrixXd}\hlstd{}\hlopt{(}\hlstd{n}\hlopt{,}\hlstd{n}\hlopt{).}\hlstd{}\hlkwd{setZero}\hlstd{}\hlopt{().}\hlstd{selfadjointView}\hlopt{$<$}\hlstd{Lower}\hlopt{$>$()}\hspace*{\fill}\\
+ \hlstd{}\hlstd{\ \ \ \ \ \ \ \ \ \ \ \ \ }\hlstd{}\hlopt{.}\hlstd{}\hlkwd{rankUpdate}\hlstd{}\hlopt{(}\hlstd{A}\hlopt{.}\hlstd{}\hlkwd{adjoint}\hlstd{}\hlopt{());}\hspace*{\fill}\\
+ \hlstd{}\hlopt{\}}\hlstd{}\hspace*{\fill}\\
+ \mbox{}
+ \normalfont
+ \normalsize
+ \hrule
+ \caption{The contents of the character vector, \code{incl}, that will preface \proglang{C++} code segments that follow.
+ \label{fig:incl}}
+\end{figure}
+
+\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^\top$ 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^\top$.
+Statisticians often write the decomposition as $\bm A=\bm R^\top\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} and the corresponding methods are
+\code{.llt()} and \code{.ldlt()}.
+Because the Cholesky decomposition involves taking square roots, we
+pass a numeric matrix, $\bm A$, not an integer matrix.
+%
+\begin{figure}[t!]
+ \hrule \smallskip
+ \noindent
+ \ttfamily
+ \hlstd{}\hlkwb{const}\hlstd{\ \ }\hlkwb{}\hlstd{LLT}\hlopt{$<$}\hlstd{MatrixXd}\hlopt{$>$\ }\hlstd{}\hlkwd{llt}\hlstd{}\hlopt{(}\hlstd{}\hlkwd{AtA}\hlstd{}\hlopt{(}\hlstd{as}\hlopt{$<$}\hlstd{MapMatd}\hlopt{$>$(}\hlstd{AA}\hlopt{)));}\hspace*{\fill}\\
+ \hlstd{}\hlkwa{return\ }\hlstd{List}\hlopt{::}\hlstd{}\hlkwd{create}\hlstd{}\hlopt{(}\hlstd{Named}\hlopt{{(}}\hlstd{}\hlstr{"L"}\hlstd{}\hlopt{{)}\ =\ }\hlstd{}\hlkwd{MatrixXd}\hlstd{}\hlopt{(}\hlstd{llt}\hlopt{.}\hlstd{}\hlkwd{matrixL}\hlstd{}\hlopt{()),}\hspace*{\fill}\\
+ \hlstd{}\hlstd{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\hlstd{Named}\hlopt{{(}}\hlstd{}\hlstr{"R"}\hlstd{}\hlopt{{)}\ =\ }\hlstd{}\hlkwd{MatrixXd}\hlstd{}\hlopt{(}\hlstd{llt}\hlopt{.}\hlstd{}\hlkwd{matrixU}\hlstd{}\hlopt{()));}\hlstd{}\hspace*{\fill}\\
+ \mbox{}
+ \normalfont
+ \normalsize
+ \hrule
+ \caption{\code{cholCpp}: Cholesky decomposition of a cross-product.
+ \label{chol}}
+\end{figure}
+%
+\begin{CodeInput}
+R> storage.mode(A) <- "double"
+R> fchol <- cxxfunction(signature(AA = "matrix"), cholCpp, "RcppEigen", incl)
+R> (ll <- fchol(A))
+\end{CodeInput}
+\begin{CodeOutput}
+$L
+ [,1] [,2]
+[1,] 3.74166 0.00000
+[2,] 8.55236 1.96396
+
+$R
+ [,1] [,2]
+[1,] 3.74166 8.55236
+[2,] 0.00000 1.96396
+\end{CodeOutput}
+\begin{CodeInput}
+R> stopifnot(all.equal(ll[[2]], chol(crossprod(A))))
+\end{CodeInput}
+
+
+
+\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^\top\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^\top=\bm X^\top\bm X$, is
+the product of its diagonal elements, as is the case for any
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rcpp -r 4224
More information about the Rcpp-commits
mailing list