library(inline)
library(RcppEigen)
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 MatrixXd& A) {
int n(A.cols());
return MatrixXd(n,n).setZero().selfadjointView<Lower>()
.rankUpdate(A.adjoint());
}
inline MatrixXd AAt(const MatrixXd& A) {
int n(A.cols());
return MatrixXd(n,n).setZero().selfadjointView<Lower>()
.rankUpdate(A);
}
'
## 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 <- '
typedef Eigen::Map<Eigen::MatrixXi> MapMati;
const MapMati B(as<MapMati>(BB));
const MapMati C(as<MapMati>(CC));
return List::create(Named("B %*% C") = B * C,
Named("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(Named("crossprod(A)") = AtA,
Named("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
storage.mode(A) <- "double"
cholCpp <- '
const LLT<MatrixXd> llt(AtA(as<MapMatd>(AA)));
return List::create(Named("L") = MatrixXd(llt.matrixL()),
Named("R") = MatrixXd(llt.matrixU()));
'
fchol <- cxxfunction(signature(AA = "matrix"), cholCpp, "RcppEigen", incl)
(ll <- fchol(A))
stopifnot(all.equal(ll[[2]], chol(crossprod(A))))
# section 3.5
cholDetCpp <- '
const MatrixXd ata(AtA(as<MapMatd>(AA)));
const double detL(MatrixXd(ata.llt().matrixL()).diagonal().prod());
const VectorXd Dvec(ata.ldlt().vectorD());
return List::create(Named("d1") = detL * detL,
Named("d2") = Dvec.prod(),
Named("ld") = Dvec.array().log().sum());
'
fdet <- cxxfunction(signature(AA = "matrix"), cholDetCpp, "RcppEigen", incl)
unlist(ll <- fdet(A))
## section 4.1
lltLSCpp <- '
const MapMatd X(as<MapMatd>(XX));
const MapVecd y(as<MapVecd>(yy));
const int n(X.rows()), p(X.cols());
const LLT<MatrixXd> llt(AtA(X));
const VectorXd betahat(llt.solve(X.adjoint() * y));
const VectorXd fitted(X * betahat);
const VectorXd resid(y - fitted);
const int df(n - p);
const double ssq(resid.squaredNorm() / double(df));
const MatrixXd vcov(ssq * llt.solve(MatrixXd::Identity(p, p)));
return List::create(Named("coefficients") = betahat,
Named("fitted.values") = fitted,
Named("residuals") = resid,
Named("s") = sqrt(ssq),
Named("df.residual") = df,
Named("rank") = p,
Named("vcov") = vcov);
'
lltLS <- cxxfunction(signature(XX = "matrix", yy = "numeric"),
lltLSCpp, "RcppEigen", incl)
data(trees, package="datasets")
str(lltFit <- with(trees, lltLS(cbind(1, log(Girth)), log(Volume))))
str(lmFit <- with(trees, lm.fit(cbind(1, log(Girth)), log(Volume))))
for (nm in c("coefficients", "residuals", "fitted.values", "rank"))
stopifnot(all.equal(lltFit[[nm]], unname(lmFit[[nm]])))
stopifnot(all.equal(unname(vcov(lm(log(Volume) ~ log(Girth), trees))),
lltFit$vcov))
## section 4.3
dd <- data.frame(f1 = gl(4, 6, labels = LETTERS[1:4]),
f2 = gl(3, 2, labels = letters[1:3]))[-(7:8), ]
xtabs(~ f2 + f1, dd) # one missing cell
mm <- model.matrix(~ f1 * f2, dd)
kappa(mm) # large condition number, indicating rank deficiency
rcond(mm) # alternative evaluation, the reciprocal condition number
(c(rank=qr(mm)$rank, p=ncol(mm))) # rank as computed in R's qr function
set.seed(1)
dd$y <- mm %*% seq_len(ncol(mm)) + rnorm(nrow(mm), sd = 0.1)
# lm detects the rank deficiency
fm1 <- lm(y ~ f1 * f2, dd)
writeLines(capture.output(print(summary(fm1), signif.stars=FALSE))[9:22])
## section 4.6
print(summary(fmPQR <- fastLm(y ~ f1 * f2, dd)), signif.stars=FALSE)
all.equal(coef(fm1), coef(fmPQR))
all.equal(unname(fitted(fm1)), fitted(fmPQR))
all.equal(unname(residuals(fm1)), residuals(fmPQR))
print(summary(fmSVD <- fastLm(y ~ f1 * f2, dd, method=4L)), signif.stars=FALSE)
all.equal(coef(fm1), coef(fmSVD))
all.equal(unname(fitted(fm1)), fitted(fmSVD))
all.equal(unname(residuals(fm1)), residuals(fmSVD))
fmVLV <- fastLm(y ~ f1 * f2, dd, method=5L)
all.equal(coef(fmSVD), coef(fmVLV))
## section 5
badtransCpp <- '
const MapMati A(as<MapMati>(AA));
return wrap(A.transpose());
'
Ai <- matrix(1:6, ncol=2L)
ftrans2 <- cxxfunction(signature(AA = "matrix"), badtransCpp, "RcppEigen", incl)
(At <- ftrans2(Ai))
all.equal(At, t(Ai))
## section 6
sparseProdCpp <- '
using Eigen::MappedSparseMatrix;
using Eigen::SparseMatrix;
const MappedSparseMatrix<double> A(as<MappedSparseMatrix<double> >(AA));
const MapVecd y(as<MapVecd>(yy));
const SparseMatrix<double> At(A.adjoint());
return List::create(Named("At") = At,
Named("Aty") = At * y);
'
sparse1 <- cxxfunction(signature(AA = "dgCMatrix", yy = "numeric"),
sparseProdCpp, "RcppEigen", incl)
data(KNex, package="Matrix")
rr <- sparse1(KNex$mm, KNex$y)
stopifnot(all.equal(rr$At, t(KNex$mm)),
all.equal(rr$Aty, as.vector(crossprod(KNex$mm, KNex$y))))
sparseLSCpp <- '
typedef Eigen::MappedSparseMatrix<double> MSpMat;
typedef Eigen::SparseMatrix<double> SpMat;
typedef Eigen::SimplicialLDLT<SpMat> SpChol;
const SpMat At(as<MSpMat>(AA).adjoint());
const VectorXd Aty(At * as<MapVecd>(yy));
const SpChol Ch(At * At.adjoint());
if (Ch.info() != Eigen::Success) return R_NilValue;
return List::create(Named("betahat") = Ch.solve(Aty),
Named("perm") = Ch.permutationP().indices());
'
sparse2 <- cxxfunction(signature(AA = "dgCMatrix", yy = "numeric"),
sparseLSCpp, "RcppEigen", incl)
str(rr <- sparse2(KNex$mm, KNex$y))
res <- as.vector(solve(Ch <- Cholesky(crossprod(KNex$mm)),
crossprod(KNex$mm, KNex$y)))
stopifnot(all.equal(rr$betahat, res))
all(rr$perm == Ch@perm) # fill-reducing permutations are different
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.