[Rcpp-commits] r3288 - pkg/RcppEigen/vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Nov 6 04:59:55 CET 2011
Author: edd
Date: 2011-11-06 04:59:52 +0100 (Sun, 06 Nov 2011)
New Revision: 3288
Added:
pkg/RcppEigen/vignettes/code.R
Modified:
pkg/RcppEigen/vignettes/RcppEigen-intro-jss.tex
Log:
changes up to end of section 3
Modified: pkg/RcppEigen/vignettes/RcppEigen-intro-jss.tex
===================================================================
--- pkg/RcppEigen/vignettes/RcppEigen-intro-jss.tex 2011-11-05 10:43:37 UTC (rev 3287)
+++ pkg/RcppEigen/vignettes/RcppEigen-intro-jss.tex 2011-11-06 03:59:52 UTC (rev 3288)
@@ -4,7 +4,7 @@
%\VignetteKeywords{linear algebra, template programming, C++, R, Rcpp}
%\VignettePackage{RcppEigen}
-\usepackage{booktabs,flafter}
+\usepackage{booktabs,flafter,bm}
\newcommand{\R}{\proglang{R}\ } % NB forces a space so not good before % fullstop etc.
\newcommand{\Cpp}{\proglang{C++}\ }
@@ -342,9 +342,15 @@
\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
+The next \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 compiled, links and loads code written in \proglang{C++} and
+supplie as text 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 text in a variable named \code{transCpp}), and \code{cxxfunction} createan
+an executable function which is assigned to \code{ftrans}. This new function
+is used on the matrix $A$ used above, and one can check that it works as intended
by compariing the output to an explicit transpose of the matrix argument.
\begin{verbatim}
> ftrans <- cxxfunction(signature(AA="matrix"), transCpp, plugin="RcppEigen")
@@ -355,7 +361,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
@@ -364,8 +369,191 @@
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{}\hlkwa{using\ }\hlstd{Eigen}\hlopt{::}\hlstd{Map}\hlopt{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwa{using\ }\hlstd{Eigen}\hlopt{::}\hlstd{MatrixXi}\hlopt{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwb{const\ }\hlstd{Map}\hlopt{$<$}\hlstd{MatrixXi}\hlopt{$>$}\hlstd{\ \ \ \ }\hlopt{}\hlstd{}\hlkwd{B}\hlstd{}\hlopt{(}\hlstd{as}\hlopt{$<$}\hlstd{Map}\hlopt{$<$}\hlstd{MatrixXi}\hlopt{$>$\ $>$(}\hlstd{BB}\hlopt{));}\hspace*{\fill}\\
+ \hlstd{}\hlkwb{const\ }\hlstd{Map}\hlopt{$<$}\hlstd{MatrixXi}\hlopt{$>$}\hlstd{\ \ \ \ }\hlopt{}\hlstd{}\hlkwd{C}\hlstd{}\hlopt{(}\hlstd{as}\hlopt{$<$}\hlstd{Map}\hlopt{$<$}\hlstd{MatrixXi}\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}
+ \normalfont
+ \normalsize
+ \end{quote}
+ \caption{\textbf{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}
+
+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
+\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 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
+
+\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}
+\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 to construct \code{AtA} and \code{AAt} in
+that code fragment are compact and elegant. To others they are
+hopelessly confusing. If you find yourself in the latter group, you
+just need to read the expression left to right. So, for example, we
+construct \code{AAt} by creating a general integer matrix of size
+$m\times m$ (where $\bm A$ is $m\times n$), ensure that all its
+elements are zero, regard it as a self-adjoint (i.e. symmetric) matrix
+using the elements in the lower triangle, then add $\bm A\bm A^\prime$
+to it and convert back to a general matrix form (i.e. the strict
+lower triangle is copied into the strict upper triangle).
+
+For these products one could use either the lower triangle or the upper
+triangle as the result will be symmetrized before it is returned.
+
+
+
+\subsection{Cholesky decomposition of the crossprod}
+\label{sec:chol}
+
+The Cholesky decomposition of the positive-definite, symmetric matrix,
+$\bm A$, can be written in several forms. Numerical analysts define
+the ``LLt'' form as the lower triangular matrix, $\bm L$, such that
+$\bm A=\bm L\bm L^\prime$ and the ``LDLt'' form as a unit lower
+triangular matrix $\bm L$ and a diagonal matrix $\bm D$ with positive
+diagonal elements such that $\bm A=\bm L\bm D\bm L^\prime$.
+Statisticians often write the decomposition as $\bm A=\bm R^\prime\bm
+R$ where $\bm R$ is an upper triangular matrix. Of course, this $\bm
+R$ is simply the transpose of $\bm L$ from the ``LLt'' form.
+
+The templated \pkg{Eigen} classes for the LLt and LDLt forms are
+called \code{LLT} and \code{LDLT}. In general, 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
+\begin{verbatim}
+storage.mode(A) <- "double"
+\end{verbatim}
+before applying the code in Figure~\ref{chol}.
+
+\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{MatrixXd}\hlopt{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwa{using\ }\hlstd{Eigen}\hlopt{::}\hlstd{LLT}\hlopt{;}\hspace*{\fill}\\
+ \hlstd{}\hlkwa{using\ }\hlstd{Eigen}\hlopt{::}\hlstd{Lower}\hlopt{;}\hspace*{\fill}\\
+ \hlstd{}\hspace*{\fill}\\
+ \hlkwb{const\ }\hlstd{Map}\hlopt{$<$}\hlstd{MatrixXd}\hlopt{$>$}\hlstd{\ \ \ }\hlopt{}\hlstd{}\hlkwd{A}\hlstd{}\hlopt{(}\hlstd{as}\hlopt{$<$}\hlstd{Map}\hlopt{$<$}\hlstd{MatrixXd}\hlopt{$>$\ $>$(}\hlstd{AA}\hlopt{));}\hspace*{\fill}\\
+ \hlstd{}\hlkwb{const\ int}\hlstd{\ \ \ \ \ \ \ \ \ \ \ \ \ }\hlkwb{}\hlstd{}\hlkwd{n}\hlstd{}\hlopt{(}\hlstd{A}\hlopt{.}\hlstd{}\hlkwd{cols}\hlstd{}\hlopt{());}\hspace*{\fill}\\
+ \hlstd{}\hlkwb{const\ }\hlstd{LLT}\hlopt{$<$}\hlstd{MatrixXd}\hlopt{$>$\ }\hlstd{}\hlkwd{llt}\hlstd{}\hlopt{(}\hlstd{}\hlkwd{MatrixXd}\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{}\hspace*{\fill}\\
+ \hlkwa{return\ }\hlstd{List}\hlopt{::}\hlstd{}\hlkwd{create}\hlstd{}\hlopt{(}\hlstd{\textunderscore }\hlopt{{[}}\hlstd{}\hlstr{"L"}\hlstd{}\hlopt{{]}\ =\ }\hlstd{}\hlkwd{MatrixXd}\hlstd{}\hlopt{(}\hlstd{llt}\hlopt{.}\hlstd{}\hlkwd{matrixL}\hlstd{}\hlopt{()),}\hspace*{\fill}\\
+ \hlstd{}\hlstd{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }\hlstd{\textunderscore }\hlopt{{[}}\hlstd{}\hlstr{"R"}\hlstd{}\hlopt{{]}\ =\ }\hlstd{}\hlkwd{MatrixXd}\hlstd{}\hlopt{(}\hlstd{llt}\hlopt{.}\hlstd{}\hlkwd{matrixU}\hlstd{}\hlopt{()));}\hlstd{}\hspace*{\fill}
+ \normalfont
+ \normalsize
+ \end{quote}
+ \caption{\textbf{cholCpp}: Cholesky decomposition of a cross-product}
+ \label{chol}
+\end{figure}
+
+\begin{verbatim}
+> fchol <- cxxfunction(signature(AA = "matrix"), cholCpp, "RcppEigen")
+> (ll <- fchol(A))
+> stopifnot(all.equal(ll[[2]], chol(crossprod(A))))
+>
+\end{verbatim}
+
\bibliography{Rcpp}
\end{document}
Added: pkg/RcppEigen/vignettes/code.R
===================================================================
--- pkg/RcppEigen/vignettes/code.R (rev 0)
+++ pkg/RcppEigen/vignettes/code.R 2011-11-06 03:59:52 UTC (rev 3288)
@@ -0,0 +1,113 @@
+
+library(inline)
+library(RcppEigen)
+
+
+## section 3.1
+(A <- matrix(1:6, ncol=2))
+str(A)
+
+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);
+'
+
+ftrans <- cxxfunction(signature(AA="matrix"), transCpp, plugin="RcppEigen")
+(At <- ftrans(A))
+stopifnot(all.equal(At, t(A)))
+
+
+
+## section 3.2
+prodCpp <- '
+using Eigen::Map;
+using Eigen::MatrixXi;
+const Map<MatrixXi> B(as<Map<MatrixXi> >(BB));
+const Map<MatrixXi> C(as<Map<MatrixXi> >(CC));
+return List::create(_["B %*% C"] = B * C,
+ _["crossprod(B, C)"] = B.adjoint() * C);
+'
+
+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)))
+
+
+
+## section 3.3
+
+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);
+'
+fcprd <- cxxfunction(signature(AA = "matrix"), crossprodCpp, "RcppEigen")
+str(crp <- fcprd(A))
+stopifnot(all.equal(crp[[1]], crossprod(A)),
+ all.equal(crp[[2]], tcrossprod(A)))
+
+
+
+## section 3.4
+
+cholCpp <- '
+using Eigen::Map;
+using Eigen::MatrixXd;
+using Eigen::LLT;
+using Eigen::Lower;
+
+const Map<MatrixXd> A(as<Map<MatrixXd> >(AA));
+const int n(A.cols());
+const LLT<MatrixXd> llt(MatrixXd(n, n).setZero().
+ selfadjointView<Lower>().rankUpdate(A.adjoint()));
+
+return List::create(_["L"] = MatrixXd(llt.matrixL()),
+ _["R"] = MatrixXd(llt.matrixU()));
+'
+
+fchol <- cxxfunction(signature(AA = "matrix"), cholCpp, "RcppEigen")
+(ll <- fchol(A))
+stopifnot(all.equal(ll[[2]], chol(crossprod(A))))
+
+
+# section 3.5
+
+cholDetCpp <- '
+using Eigen::Lower;
+using Eigen::Map;
+using Eigen::MatrixXd;
+using Eigen::VectorXd;
+
+const Map<MatrixXd> A(as<Map<MatrixXd> >(AA));
+const int n(A.cols());
+const MatrixXd AtA(MatrixXd(n, n).setZero().
+ selfadjointView<Lower>().rankUpdate(A.adjoint()));
+const MatrixXd Lmat(AtA.llt().matrixL());
+const double detL(Lmat.diagonal().prod());
+const VectorXd Dvec(AtA.ldlt().vectorD());
+
+return List::create(_["d1"] = detL * detL,
+ _["d2"] = Dvec.prod(),
+ _["ld"] = Dvec.array().log().sum());
+'
+
+fdet <- cxxfunction(signature(AA = "matrix"), cholDetCpp, "RcppEigen")
+unlist(ll <- fdet(A))
+
More information about the Rcpp-commits
mailing list