[Rcpp-commits] r3355 - in pkg/RcppEigen: inst/doc vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Nov 14 21:07:50 CET 2011
Author: dmbates
Date: 2011-11-14 21:07:50 +0100 (Mon, 14 Nov 2011)
New Revision: 3355
Added:
pkg/RcppEigen/vignettes/RcppEigen-intro-nojss.Rnw
Removed:
pkg/RcppEigen/inst/doc/RcppEigen-intro-nojss.tex
Log:
Now moving things back into the vignettes directory, sigh.
Deleted: pkg/RcppEigen/inst/doc/RcppEigen-intro-nojss.tex
===================================================================
--- pkg/RcppEigen/inst/doc/RcppEigen-intro-nojss.tex 2011-11-14 19:52:08 UTC (rev 3354)
+++ pkg/RcppEigen/inst/doc/RcppEigen-intro-nojss.tex 2011-11-14 20:07:50 UTC (rev 3355)
@@ -1,1557 +0,0 @@
-
-\documentclass[shortnames,article,nojss]{jss}
-%\VignetteIndexEntry{RcppEigen-intro}
-%\VignetteKeywords{linear algebra, template programming, C++, R, Rcpp}
-%\VignettePackage{RcppEigen}
-
-%\usepackage{booktabs,flafter,bm,amsmath}
-\usepackage{booktabs,bm,amsmath}
-
-\newcommand{\R}{\proglang{R}\ } % NB forces a space so not good before % fullstop etc.
-\newcommand{\Cpp}{\proglang{C++}\ }
-
-\author{Douglas Bates\\Univ.~of Wisconsin/Madison \And Dirk Eddelbuettel\\Debian Project} % \And Romain Fran\c{c}ois\\R Enthusiasts}
-\title{Fast and Elegant Numerical Linear Algebra\\ Using the \pkg{RcppEigen} Package}
-
-\Plainauthor{Douglas Bates, Dirk Eddelbuettel} %, Romain Fran\c{c}ois}
-\Plaintitle{Fast and Elegant Numerical Linear Algebra Using the RcppEigen Package}
-\Shorttitle{Fast and Elegant Numerical Linear Algebra with RcppEigen}
-
-\Abstract{
- \noindent
- 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}} %% at least one keyword must be supplied
-\Plainkeywords{Linear algebra, template programmig, R, C++, Rcpp} %% without formatting
-
-\Address{
- Douglas Bates \\
- Department of Statistics \\
- University of Wisconsin - Madison \\
- Madison, WI, USA \\
- E-mail: \email{bates at stat.wisc.edu} \\
- URL: \url{http://www.stat.wisc.edu/~bates/}\\
-
- Dirk Eddelbuettel \\
- Debian Project \\
- River Forest, IL, USA\\
- E-mail: \email{edd at debian.org}\\
- URL: \url{http://dirk.eddelbuettel.com}\\
-
- % Romain Fran\c{c}ois\\
- % Professional R Enthusiast\\
- % 1 rue du Puits du Temple, 34 000 Montpellier\\
- % FRANCE \\
- % E-mail: \email{romain at r-enthusiasts.com}\\
- % URL: \url{http://romainfrancois.blog.free.fr}
-}
-
-%% need no \usepackage{Sweave.sty}
-
-\newcommand{\rank}{\operatorname{rank}}
-
-%% highlights macros
-%% Style definition file generated by highlight 2.7, http://www.andre-simon.de/
-% Highlighting theme definition:
-\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.74,0.55,0.55}{#1}}
-\newcommand{\hlstr}[1]{\textcolor[rgb]{0.90,0.15,0.15}{#1}}
-%green: \newcommand{\hlstr}[1]{\textcolor[rgb]{0.13,0.67,0.13}{#1}} % 0.74 -> % 0.90; 0.55 -> 0.25
-\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}
-
-
-% ------------------------------------------------------------------------
-
-\begin{document}
-
-%\SweaveOpts{engine=R,eps=FALSE}
-
-\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.
-% Could cite Eddelbuettel (1996) here, but a real survey would be better.
-
-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 ``Eigen3'' being the current version.
-
-\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 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 \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, along with 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. 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}[tb]
- \centering
- \begin{tabular}{l l}
- \toprule
- \multicolumn{1}{c}{\proglang{R} object type} & \multicolumn{1}{c}{\pkg{Eigen} class typedef}\\
- \midrule
- 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>}\\
- \bottomrule
- \end{tabular}
- \caption{Correspondence between 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
-
-%% Alternatively, use 'highlight --enclose-pre --no-doc --latex --style=emacs --syntax=C++'
-%% as the command invoked from C-u M-|
-%% For version 3.5 of highlight this should be
-%% highlight --enclose-pre --no-doc --out-format=latex --syntax=C++
-%%
-%% keep one copy to redo later
-%%
-%% using Eigen::MatrixXd;
-%%
-\begin{quote}
- \noindent
- \ttfamily
- \hlstd{}\hlkwa{using\ }\hlstd{Eigen}\hlopt{::}\hlstd{MatrixXd}\hlopt{;}\hlstd{}\hspace*{\fill}\\
- \mbox{}
- \normalfont
- \normalsize
-\end{quote}
-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 \texttt{`*'}
-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}
-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}
-
-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{verbatim}
-> (A <- matrix(1:6, ncol=2))
- [,1] [,2]
-[1,] 1 4
-[2,] 2 5
-[3,] 3 6
-> str(A)
- int [1:3, 1:2] 1 2 3 4 5 6
-\end{verbatim}
-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}[htb]
- %\begin{quote}
- \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
- %\end{quote}
- \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{verbatim}
-> ftrans <- cxxfunction(signature(AA="matrix"), transCpp, plugin="RcppEigen")
-> (At <- ftrans(A))
- [,1] [,2] [,3]
-[1,] 1 2 3
-[2,] 4 5 6
-> stopifnot(all.equal(At, t(A)))
-\end{verbatim}
-
-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 Sec.~\ref{sec:arrays}, the \code{`*'} operator
-performs matrix multiplication on \code{Eigen::Matrix} or
-\code{Eigen::Vector} objects. The \proglang{C++} code in
-Figure~\ref{prod} produces a list containing both the product and
-cross-product of its two arguments.
-
-\begin{figure}[htb]
- %\begin{quote}
- \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
- %\end{quote}
- \caption{\code{prodCpp}: Product and cross-product of two matrices.}
- \label{prod}
-\end{figure}
-
-\begin{verbatim}
-> 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))
-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
-> stopifnot(all.equal(fp[[1]], B %*% C), all.equal(fp[[2]], crossprod(B, C)))
-\end{verbatim}
-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^\prime\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^\prime$ to the current symmetric matrix. The scalar
-multiple defaults to 1. The code in Figure~\ref{crossprod} produces
-
-\begin{figure}[htb]
- %\begin{quote}
- \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
- %\end{quote}
- \caption{\code{crossprodCpp}: Cross-product and transposed cross-product of a single matrix.}
- \label{crossprod}
-\end{figure}
-\begin{verbatim}
-> fcprd <- cxxfunction(signature(AA = "matrix"), crossprodCpp, "RcppEigen")
-> str(crp <- fcprd(A))
-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
-> stopifnot(all.equal(crp[[1]], crossprod(A)),
-+ all.equal(crp[[2]], tcrossprod(A)))
-\end{verbatim}
-
-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 (\textit{i.e.} symmetric) matrix
-using the elements in the lower triangle, then adding $\bm A\bm A^\prime$
-to it and converting back to a general matrix form (\textit{i.e.} the strict lower
-triangle is copied into the strict upper triangle).
-
-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}[htb]
- %\begin{quote}
- \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
-%\end{quote}
- \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^\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} 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}[htb]
- %\begin{quote}
- \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
- %\end{quote}
- \caption{\code{cholCpp}: Cholesky decomposition of a cross-product.}
- \label{chol}
-\end{figure}
-
-\begin{verbatim}
-> storage.mode(A) <- "double"
-> fchol <- cxxfunction(signature(AA = "matrix"), cholCpp, "RcppEigen", incl)
-> (ll <- fchol(A))
-$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
-
-> stopifnot(all.equal(ll[[2]], chol(crossprod(A))))
-\end{verbatim}
-
-
-
-\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 using 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 it is known that the diagonal
-elements of $\bm D$ must be non-negative, one often evaluates the
-logarithm of the determinant as the sum of the logarithms of the
-diagonal elements of $\bm D$. Several options are shown in
-Figure~\ref{cholDet}.
-
-\begin{figure}[htb]
- %\begin{quote}
- \noindent
- \ttfamily
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rcpp -r 3355
More information about the Rcpp-commits
mailing list