[Rcpp-commits] r3347 - pkg/RcppEigen/vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Nov 14 20:28:06 CET 2011
Author: dmbates
Date: 2011-11-14 20:28:06 +0100 (Mon, 14 Nov 2011)
New Revision: 3347
Removed:
pkg/RcppEigen/vignettes/RcppEigen-Intro.Rnw
Log:
Remove the older .Rnw file
Deleted: pkg/RcppEigen/vignettes/RcppEigen-Intro.Rnw
===================================================================
--- pkg/RcppEigen/vignettes/RcppEigen-Intro.Rnw 2011-11-14 17:47:43 UTC (rev 3346)
+++ pkg/RcppEigen/vignettes/RcppEigen-Intro.Rnw 2011-11-14 19:28:06 UTC (rev 3347)
@@ -1,1532 +0,0 @@
-\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{booktabs,flafter,thumbpdf} % booktabs for nice tables, flafter for float control, thumbpdf is a tip from Achim
-\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}}
-
-%% 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}
-
-
-<<version,echo=FALSE,print=FALSE>>=
-prettyVersion <- packageDescription("RcppEigen")$Version
-prettyDate <- format(Sys.Date(), "%B %e, %Y")
-options(prompt="> ") # to overcome DE's default of "R> "
-@
-\author{Douglas Bates \and Dirk Eddelbuettel \and Romain Fran\c{c}ois}
-\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}
-
-%% Below is what DE currently uses somewhere else -- I am not religuous about
-%% the grey background etc -- comment out and see if you like any of it
-%% Just like DB below, DE also uses frame=tb for top+bottom frame lines
-\definecolor{lightgray}{rgb}{0.975,0.975,0.975}
-%\lstset{backgroundcolor=\color{lightgray}}
-% \lstset{numbers=left, numberstyle=\tiny, stepnumber=2, numbersep=5pt}
-%\lstset{keywordstyle=\color{black}\bfseries\tt}
-% \lstset{ %
-% %language=Octave, % the language of the code
-% %basicstyle=\footnotesize, % the size of the fonts that are used for the code
-% basicstyle=\small, % the size of the fonts that are used for the code
-% numbers=left, % where to put the line-numbers
-% %numberstyle=\footnotesize, % the size of the fonts that are used for the line-numbers
-% stepnumber=2, % the step between two line-numbers. If it's 1, each line
-% % will be numbered
-% numbersep=5pt, % how far the line-numbers are from the code
-% %backgroundcolor=\color{white}, % choose the background color. You must add \usepackage{color}
-% showspaces=false, % show spaces adding particular underscores
-% showstringspaces=false, % underline spaces within strings
-% showtabs=false, % show tabs within strings adding particular underscores
-% %frame=single, % adds a frame around the code
-% tabsize=2, % sets default tabsize to 2 spaces
-% captionpos=b, % sets the caption-position to bottom
-% breaklines=true, % sets automatic line breaking
-% breakatwhitespace=false % sets if automatic breaks should only happen at whitespace
-% %title=\lstname, % show the filename of files included with \lstinputlisting;
-% % also try caption instead of title
-% %escapeinside={\%*}{*)}, % if you want to add a comment within your code
-% %morekeywords={*,...} % if you want to add more keywords to the set
-% }
-
-
-\maketitle
-
-\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.
-}
-
-\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++}, 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 (though a more
-detailed discussions of template programming is however beyond this
-paper). \cite{Veldhuizen:1998:Blitz} provided an early and influential
-implementation that already demonstrated key features of this approach. Its
-usage however was held back at the time by the somewhat limited availability
-of compilers implementing all necessary features of the \proglang{C++}
-language.
-
-This situation has greatly improved over the last decade, and many more such
-libraries have been contributed. One such \proglang{C++} library is
-\pkg{Eigen} by \citet*{Eigen:Web}. \pkg{Eigen} 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 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. The next section introduces some
-of these classes and shows how to interface to them from \proglang{R}.
-
-\section{Eigen classes}
-\label{sec:eclasses}
-
-\pkg{Eigen} \citep*{Eigen:Web} 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 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]
- \caption{Correspondence between R matrix and vector types and classes in the \code{Eigen} namespace.}
- \label{tab:REigen}
- \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}
-\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 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}
-\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 one can use the Eigen templated class \code{Map}, and for
-sparse matrices one can deploy the 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}
-\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 \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 one needs 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
-% 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 R object (called an \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::MatrixXi> MapMati;
-% const MapMati A(Rcpp::as<MapMati>(AA));
-\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{A}\hlstd{}\hlopt{(}\hlstd{Rcpp}\hlopt{::}\hlstd{as}\hlopt{$<$}\hlstd{MapMati}\hlopt{$>$(}\hlstd{AA}\hlopt{));}\hlstd{}\hspace*{\fill}\\
- \mbox{}
- \normalfont
- \normalsize
-\end{quote}
-
-The \Sexpr{link("cxxfunction")} 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
-\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
-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
-<<Adef>>=
-(A <- matrix(1:6, ncol=2))
-str(A)
-@
-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 mapped to an
-\code{Eigen::MatrixXi} object then the matrix \code{At} is constructed
-from its transpose and returned to \proglang{R}.
-
-<<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{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{\textbf{transCpp}: Transpose a matrix of integers}
- \label{trans}
-\end{figure}
-
-The next \proglang{R} snippet compiles and links the \proglang{C++} code
-segment (stored as text in a variable named \code{transCpp}) into an
-executable function \code{ftrans} and then checks that it works as intended
-by comparing the output to an explicit transpose of the matrix argument.
-<<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, in what follows,
-the \code{adjoint()} method is used 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.
-
-<<prodCpp,echo=FALSE>>=
-prodCpp <- '
-typedef Eigen::Map<Eigen::MatrixXi> MapMati;
-const MapMati B(as<MapMati>(BB));
-const MapMati C(as<MapMati>(CC));
-return List::create(_["B %*% C"] = B * C,
- _["crossprod(B, C)"] = B.adjoint() * C);
-'
-@
-\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{\textunderscore }\hlopt{{[}}\hlstd{}\hlstr{"B\ \%{*}\%\ C"}\hlstd{}\hlopt{{]}}\hlstd{\ \ \ \ \ \ \ \ \ }\hlopt{=\ }\hlstd{B\ }\hlopt{{*}\ }\hlstd{C}\hlopt{,}\hspace*{\fill}\\
- \hlstd{}\hlstd{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\hlstd{\textunderscore }\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{\textbf{prodCpp}: Product and cross-product of two matrices}
- \label{prod}
-\end{figure}
-<<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$. One 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, a \code{SelfAdjointView} is created,
-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 Figure~\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{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{\textunderscore }\hlopt{{[}}\hlstd{}\hlstr{"crossprod(A)"}\hlstd{}\hlopt{{]}}\hlstd{\ \ }\hlopt{=\ }\hlstd{AtA}\hlopt{,}\hspace*{\fill}\\
- \hlstd{}\hlstd{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\hlstd{\textunderscore }\hlopt{{[}}\hlstd{}\hlstr{"tcrossprod(A)"}\hlstd{}\hlopt{{]}\ =\ }\hlstd{AAt}\hlopt{);}\hlstd{}\hspace*{\fill}
- \normalfont
- \normalsize
- \end{quote}
- \caption{\textbf{crossprodCpp}: Cross-product and transposed cross-product of a single matrix}
- \label{crossprod}
-\end{figure}
-<<>>=
-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 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$), 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 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 string 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.
-<<AtA,echo=FALSE>>=
-incl <- '
-using Eigen::LLT;
-using Eigen::Lower;
-using Eigen::Map;
-using Eigen::MatrixXd;
-using Eigen::MatrixXi;
-using Eigen::Upper;
-using Eigen::VectorXd;
-typedef Map<MatrixXd> MapMatd;
-typedef Map<MatrixXi> MapMati;
-typedef Map<VectorXd> MapVecd;
-inline MatrixXd AtA(const MapMatd& A) {
- int n(A.cols());
- return MatrixXd(n,n).setZero().selfadjointView<Lower>()
- .rankUpdate(A.adjoint());
-}
-'
-@
-\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}
-\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, one would preserve the
-objects from these classes in order to re-use them for solutions of
-linear systems. For a simple illustration, the matrix $\bm L$
-from the ``LLt'' form is returned.
-
-Because the Cholesky decomposition involves taking square roots, the internal
-representation is switched to numeric matrices
-<<storage>>=
-storage.mode(A) <- "double"
-@
-before applying the code in Figure~\ref{chol}.
-<<cholCpp,echo=FALSE>>=
-cholCpp <- '
-const LLT<MatrixXd> llt(AtA(as<MapMatd>(AA)));
-return List::create(_["L"] = MatrixXd(llt.matrixL()),
- _["R"] = MatrixXd(llt.matrixU()));
-'
-@
-\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}\\
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rcpp -r 3347
More information about the Rcpp-commits
mailing list