From 7f3f61285b18ab3c29aca2b2f7c7cd8ecc51bba1 Mon Sep 17 00:00:00 2001 From: Andreas Alfons Date: Wed, 17 Oct 2012 21:26:34 +0200 Subject: [PATCH 01/14] added code for sparse maximum correlation (k = 1 is forced) --- NAMESPACE | 2 + R/cca.R | 51 +++++++- man/ccaGrid.Rd | 14 +++ src/cca.cpp | 329 +++++++++++++++++++++++++++++++++++++++++++++++-- 4 files changed, 382 insertions(+), 14 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 12067cc..e4236cd 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -15,6 +15,8 @@ export(fastMedian) export(maxCorGrid) export(maxCorProj) export(permTest) +export(sCCAgrid) +export(sccaGrid) import(Rcpp) import(RcppArmadillo) import(parallel) diff --git a/R/cca.R b/R/cca.R index 4708c51..0af810a 100644 --- a/R/cca.R +++ b/R/cca.R @@ -142,7 +142,34 @@ ccaGrid <- function(x, y, k = 1, cca } -## wrapper function for more compatibility with package pcaPP +#' @rdname ccaGrid +#' @export +sccaGrid <- function(x, y, k = 1, lambda = 0, + method = c("spearman", "kendall", "quadrant", "M", "pearson"), + control = list(...), nIterations = 10, nAlternate = 10, nGrid = 25, + select = NULL, tol = 1e-06, fallback = FALSE, seed = NULL, ...) { + ## initializations + matchedCall <- match.call() + ## define list of control arguments for algorithm + lambda <- rep(as.numeric(lambda), length.out=2) + lambda[is.na(lambda)] <- formals()$lambda + nIterations <- as.integer(nIterations) + nAlternate <- as.integer(nAlternate) + nGrid <- as.integer(nGrid) + tol <- as.numeric(tol) + ppControl <- list(lambda=lambda, nIterations=nIterations, + nAlternate=nAlternate, nGrid=nGrid, select=select, tol=tol) + ## call workhorse function + cca <- ccaPP(x, y, k, method=method, corControl=control, algorithm="sparse", + ppControl=ppControl, fallback=fallback, seed=seed) + cca$lambda <- lambda + cca$call <- matchedCall + cca +} + + +## wrapper functions for more compatibility with package pcaPP + #' @rdname ccaGrid #' @export @@ -160,6 +187,23 @@ CCAgrid <- function(x, y, k = 1, cca } +#' @rdname ccaGrid +#' @export + +sCCAgrid <- function(x, y, k = 1, lambda = 0, + method = c("spearman", "kendall", "quadrant", "M", "pearson"), + maxiter = 10, maxalter = 10, splitcircle = 25, select=NULL, + zero.tol = 1e-06, fallback = FALSE, seed = NULL, ...) { + ## initializations + matchedCall <- match.call() + ## call ccaGrid() + cca <- sccaGrid(x, y, k=k, lambda=lambda, method=method, nIterations=maxiter, + nAlternate=maxalter, nGrid=splitcircle, select=select, tol=zero.tol, + fallback=fallback, seed=seed, ...) + cca$call <- matchedCall + cca +} + #' (Robust) CCA via projections through the data points #' @@ -278,7 +322,8 @@ CCAproj <- function(x, y, k = 1, ## workhorse function ccaPP <- function(x, y, k = 1, method = c("spearman", "kendall", "quadrant", "M", "pearson"), - corControl, forceConsistency = TRUE, algorithm = c("grid", "proj"), + corControl, forceConsistency = TRUE, + algorithm = c("grid", "sparse", "proj"), ppControl, fallback = FALSE, seed = NULL) { ## initializations x <- as.matrix(x) @@ -303,7 +348,7 @@ ccaPP <- function(x, y, k = 1, method <- match.arg(method) corControl <- getCorControl(method, corControl, forceConsistency) # additional checks for grid search algorithm - if(algorithm == "grid") { + if(algorithm != "proj") { # check subset of variables to be used for determining the order of # the variables from the respective other data set select <- ppControl$select diff --git a/man/ccaGrid.Rd b/man/ccaGrid.Rd index 94b112d..1b72f81 100644 --- a/man/ccaGrid.Rd +++ b/man/ccaGrid.Rd @@ -2,6 +2,8 @@ \alias{ccaGrid} \alias{CCAgrid} \alias{print.cca} +\alias{sccaGrid} +\alias{sCCAgrid} \title{(Robust) CCA via alternating series of grid searches} \usage{ ccaGrid(x, y, k = 1, @@ -10,11 +12,23 @@ nGrid = 25, select = NULL, tol = 1e-06, fallback = FALSE, seed = NULL, ...) + sccaGrid(x, y, k = 1, lambda = 0, + method = c("spearman", "kendall", "quadrant", "M", "pearson"), + control = list(...), nIterations = 10, nAlternate = 10, + nGrid = 25, select = NULL, tol = 1e-06, + fallback = FALSE, seed = NULL, ...) + CCAgrid(x, y, k = 1, method = c("spearman", "kendall", "quadrant", "M", "pearson"), maxiter = 10, maxalter = 10, splitcircle = 25, select = NULL, zero.tol = 1e-06, fallback = FALSE, seed = NULL, ...) + + sCCAgrid(x, y, k = 1, lambda = 0, + method = c("spearman", "kendall", "quadrant", "M", "pearson"), + maxiter = 10, maxalter = 10, splitcircle = 25, + select = NULL, zero.tol = 1e-06, fallback = FALSE, + seed = NULL, ...) } \arguments{ \item{x,y}{each can be a numeric vector, matrix or data diff --git a/src/cca.cpp b/src/cca.cpp index fb7a9dd..a94a807 100644 --- a/src/cca.cpp +++ b/src/cca.cpp @@ -420,10 +420,11 @@ vec GridControl::getGrid(const uword& i) { // grid search to update one weighting vector // x ............ data matrix for which to update the weighting vector +// orderX ......... order in which to browse through the variables // y ............ linear combination of the other data matrix according to the // other weighting vector, which is kept fixed // corControl ... control object to compute correlation -// grid ........ grid points to be used for grid search +// grid ......... grid points to be used for grid search // maxCor ....... maximum correlation to be updated // a ............ weighting vector to be updated template @@ -443,7 +444,7 @@ void GridControl::gridSearch(const mat& x, const uvec& orderX, const vec& y, vec currentA = cos(angle) * a + sin(angle) * ej; corY(k) = abs(corControl.cor(x * currentA, y)); } - // find grid point that maximizes correlation functional and keep + // find grid point that maximizes the correlation functional and keep // maximum correlation of current grid search uword whichMax; double currentMaxCor = corY.max(whichMax); @@ -490,7 +491,12 @@ double GridControl::maxCor(const mat& x, const mat& y, CorControl corControl, // both data sets are univariate a.ones(p); b.ones(q); vec xx = x.unsafe_col(0), yy = y.unsafe_col(0); // reuse memory - maxCor = abs(corControl.cor(xx, yy)); // compute correlation + maxCor = corControl.cor(xx, yy); // compute correlation + // check sign of correlation + if(maxCor < 0) { + maxCor = -maxCor; + b = -b; + } } else { double previousMaxCor = R_NegInf; if((p > 1) && (q == 1)) { @@ -580,17 +586,292 @@ double GridControl::maxCor(const mat& x, const mat& y, CorControl corControl, } else { return NA_REAL; // should never happen } + // ensure that norm of weighting vectors is 1 + a = a / norm(a, 2); + b = b / norm(b, 2); + // check direction + double r = corControl.cor(x * a, y * b); + if(r < 0) { + if((p > 1) && (q == 1)) { + a = -a; + } else { + b = -b; + } + } } - // ensure that norm of weighting vectors is 1 - a = a / norm(a, 2); - b = b / norm(b, 2); - // check direction - double r = corControl.cor(x * a, y * b); - if(r < 0) { + // return maximum correlation + return maxCor; +} + + +// ------------------------------------------------ +// control class for sparse alternate grid searches +// ------------------------------------------------ + +class SparseGridControl: public GridControl { +public: + double lambdaX, lambdaY; // penalty parameters + // constructors + SparseGridControl(); + SparseGridControl(List&); + // grid search to update one weighting vector + template + void gridSearch(const mat&, const uvec&, const double&, const vec&, + CorControl, vec&, double&, vec&); + // grid search to update one weighting vector for multivariate x and y + template + void gridSearch(const mat&, const uvec&, const double&, const vec&, + const double&, CorControl, vec&, double&, vec&, double&, double&); + // maximum correlation between multivariate data sets x and y + template + double maxCor(const mat&, const mat&, CorControl, vec&, vec&); +}; + +// constructors +inline SparseGridControl::SparseGridControl() { + lambdaX = 0; + lambdaY = 0; +} +inline SparseGridControl::SparseGridControl(List& control) +: GridControl(control) { + SEXP R_lambda = control["lambda"]; + NumericVector lambda(R_lambda); + lambdaX = lambda[0]; + lambdaY = lambda[1]; +} + +// sparse grid search to update one weighting vector +// x .............. data matrix for which to update the weighting vector +// orderX ......... order in which to browse through the variables +// lambda ......... penalty parameter for the weighting vector to be updated +// y .............. linear combination of the other data matrix according to +// the other weighting vector, which is kept fixed +// corControl ..... control object to compute correlation +// grid ........... grid points to be used for grid search +// maxObjective ... maximum correlation to be updated +// a .............. weighting vector to be updated +template +void SparseGridControl::gridSearch(const mat& x, const uvec& orderX, + const double& lambda, const vec& y, CorControl corControl, + vec& grid, double& maxObjective, vec& a) { + // initializations + const uword p = x.n_cols, nGrid = grid.n_elem; + // perform grid searches for each canonical basis vector + for(uword j = 0; j < p; j++) { + // define current canonical basis vector according to order of columns + vec ej = zeros(p); + ej(orderX(j)) = 1; + // perform grid search for the current canonical basis vector + vec objective(nGrid); + for(uword k = 0; k < nGrid; k++) { + double angle = grid(k); + vec currentA = cos(angle) * a + sin(angle) * ej; + currentA = currentA / norm(currentA, 2); + double currentCorY = abs(corControl.cor(x * currentA, y)); + objective(k) = currentCorY - lambda * norm(currentA, 1); + } + // find grid point that maximizes the penalized correlation functional + // and keep maximum of current grid search + uword whichMax; + double currentMaxObjective = objective.max(whichMax); + // update maximum and weighting vector + // if 0 degree angle is not part of the grid, the penalized maximum + // correlation of the current grid search may be smaller than the + // previous maximum + if(currentMaxObjective > maxObjective) { + maxObjective = currentMaxObjective; + double optAngle = grid(whichMax); + a = cos(optAngle) * a + sin(optAngle) * ej; + } + } +} + +// sparse grid search to update one weighting vector for multivariate x and y +// x .............. data matrix for which to update the weighting vector +// orderX ......... order in which to browse through the variables +// lambda ......... penalty parameter for the weighting vector to be updated +// y .............. linear combination of the other data matrix according to +// the other weighting vector, which is kept fixed +// penaltyY ....... penalty for the weighting vector to kept fixed +// corControl ..... control object to compute correlation +// grid ........... grid points to be used for grid search +// maxCor ......... maximum correlation to be updated +// maxObjective ... maximum value of the objective function to be updated +// a .............. weighting vector to be updated +// penaltyX ....... penalty for the weighting vector to be updated +template +void SparseGridControl::gridSearch(const mat& x, const uvec& orderX, + const double& lambda, const vec& y, const double& penaltyY, + CorControl corControl, vec& grid, double& maxCor, vec& a, + double& penaltyX, double& maxObjective) { + // initializations + const uword p = x.n_cols, nGrid = grid.n_elem; + // perform grid searches for each canonical basis vector + for(uword j = 0; j < p; j++) { + // define current canonical basis vector according to order of columns + vec ej = zeros(p); + ej(orderX(j)) = 1; + // perform grid search for the current canonical basis vector + vec corY(nGrid), objective(nGrid), penalties(nGrid); + for(uword k = 0; k < nGrid; k++) { + double angle = grid(k); + vec currentA = cos(angle) * a + sin(angle) * ej; + currentA = currentA / norm(currentA, 2); + corY(k) = abs(corControl.cor(x * currentA, y)); + penalties(k) = lambda * norm(currentA, 1); + objective(k) = corY(k) - penalties(k) - penaltyY; + } + // find grid point that maximizes the penalized correlation functional + // and keep maximum of current grid search + uword whichMax; + double currentMaxObjective = objective.max(whichMax); + // update maximum and weighting vector + // if 0 degree angle is not part of the grid, the penalized maximum + // correlation of the current grid search may be smaller than the + // previous maximum + if(currentMaxObjective > maxObjective) { + maxCor = corY(whichMax); + double optAngle = grid(whichMax); + a = cos(optAngle) * a + sin(optAngle) * ej; + penaltyX = penalties(whichMax); + maxObjective = currentMaxObjective; + } + } +} + +// maximum correlation between multivariate data sets x and y based on +// sparse alternate grid searches +// x ............. first data matrix +// y ............. second data matrix +// corControl .... control object to compute correlation +// a ............. first weighting vector to be updated +// b ............. second weighting vector to be updated +template +double SparseGridControl::maxCor(const mat& x, const mat& y, + CorControl corControl, vec& a, vec& b) { + // initializations + uword p = x.n_cols, q = y.n_cols; + // perform alternate grid searches if both data sets are multivariate + // if one data set is univariate, alternate grid searches are not necessary + double maxCor; // initialize maximum correlation + if((p == 1) && (q == 1)) { + // both data sets are univariate + a.ones(p); b.ones(q); + vec xx = x.unsafe_col(0), yy = y.unsafe_col(0); // reuse memory + maxCor = corControl.cor(xx, yy); // compute correlation + // check sign of correlation + if(maxCor < 0) { + maxCor = -maxCor; + b = -b; + } + } else { + double maxObjective, previousMaxObjective = R_NegInf; if((p > 1) && (q == 1)) { - a = -a; + // x is multivariate, y is univariate + vec yy = y.unsafe_col(0); // reuse memory + uvec orderX(p); + a.zeros(p); b.ones(q); + findOrder(x, yy, corControl, orderX, maxCor, a); // column order + maxObjective = maxCor - lambdaX; // L1 norm of a is 1 + // stop if there are two consecutive iterations without improvement + uword i = 0, convCounter = 0; + while((i < nIterations) && (convCounter < 2)) { + previousMaxObjective = maxObjective; + vec grid = getGrid(i+1); // define vector of grid points + gridSearch(x, orderX, lambdaX, yy, corControl, + grid, maxObjective, a); // grid search + i++; + setCounter(convCounter, maxObjective, previousMaxObjective); + } + } else if((p == 1) && (q > 1)) { + // x is univariate, y is multivariate + vec xx = x.unsafe_col(0); // reuse memory + uvec orderY(q); + a.ones(p); b.zeros(q); + findOrder(y, xx, corControl, orderY, maxCor, b); // column order + maxObjective = maxCor - lambdaY; // L1 norm of b is 1 + // stop if there are two consecutive iterations without improvement + uword i = 0, convCounter = 0; + while((i < nIterations) && (convCounter < 2)) { + previousMaxObjective = maxObjective; + vec grid = getGrid(i+1); // define vector of grid points + gridSearch(y, orderY, lambdaY, xx, corControl, + grid, maxObjective, b); // grid search + i++; + setCounter(convCounter, maxObjective, previousMaxObjective); + } + } else if((p > 1) && (q > 1)) { + // both data sets are multivariate + // compute correlations between variables in x and y + uvec orderX(p), orderY(q); + a.zeros(p); b.zeros(q); + bool startWithX; + findOrder(x, y, corControl, orderX, orderY, + maxCor, a, b, startWithX); // column orders + double penaltyX = lambdaX, penaltyY = lambdaY; // L1 norms are 1 + maxObjective = maxCor - penaltyX - penaltyY; + // perform alternate grid searches + if(startWithX) { + // start with grid search for x + // stop if there are two consecutive iterations without improvement + uword i = 0, convCounter = 0; + while((i < nIterations) && (convCounter < 2)) { + previousMaxObjective = maxObjective; + vec grid = getGrid(i+1); // define vector of grid points + uword j = 0; + double altMaxObjective = R_NegInf; + while((j < nAlternate) && ((maxObjective - altMaxObjective) > tol)) { + altMaxObjective = maxObjective; + // maximize correlation functional over a keeping b fixed + vec yb = y * b; // linear combination of columns of y + gridSearch(x, orderX, lambdaX, yb, penaltyY, corControl, + grid, maxCor, a, penaltyX, maxObjective); + // maximize correlation functional over b keeping a fixed + vec xa = x * a; // linear combination of columns of x + gridSearch(y, orderY, lambdaY, xa, penaltyX, corControl, + grid, maxCor, b, penaltyY, maxObjective); + j++; + } + i++; + setCounter(convCounter, maxObjective, previousMaxObjective); + } + } else { + // start with grid search for y + // stop if there are two consecutive iterations without improvement + uword i = 0, convCounter = 0; + while((i < nIterations) && (convCounter < 2)) { + previousMaxObjective = maxObjective; + vec grid = getGrid(i+1); // define vector of grid points + uword j = 0; + double altMaxObjective = R_NegInf; + while((j < nAlternate) && ((maxObjective - altMaxObjective) > tol)) { + altMaxObjective = maxObjective; + // maximize correlation functional over b keeping a fixed + vec xa = x * a; // linear combination of columns of x + gridSearch(y, orderY, lambdaY, xa, penaltyX, corControl, + grid, maxCor, b, penaltyY, maxObjective); + // maximize correlation functional over a keeping b fixed + vec yb = y * b; // linear combination of columns of y + gridSearch(x, orderX, lambdaX, yb, penaltyY, corControl, + grid, maxCor, a, penaltyX, maxObjective); + j++; + } + i++; + setCounter(convCounter, maxObjective, previousMaxObjective); + } + } } else { - b = -b; + return NA_REAL; // should never happen + } + // check direction + maxCor = corControl.cor(x * a, y * b); + if(maxCor < 0) { + maxCor = -maxCor; + if((p > 1) && (q == 1)) { + a = -a; + } else { + b = -b; + } } } // return maximum correlation @@ -924,6 +1205,32 @@ SEXP R_ccaPP(SEXP R_x, SEXP R_y, SEXP R_k, SEXP R_method, SEXP R_corControl, } else { error("method not available"); } + } else if(algorithm == "sparse") { + // TODO: algorithm currently only works for k = 1 + // larger values of k require backtransformation of weighting vectors + // before computing the penalty + k = 1; + // define control object for alternate grid searches + SparseGridControl ppControl(Rcpp_ppControl); + // define control object for the correlations and call the arma version + if(method == "spearman") { + CorSpearmanControl corControl(Rcpp_corControl); + r = ccaPP(x, y, k, corControl, ppControl, true, fallback, A, B); + } else if(method == "kendall") { + CorKendallControl corControl(Rcpp_corControl); + r = ccaPP(x, y, k, corControl, ppControl, true, fallback, A, B); + } else if(method == "quadrant") { + CorQuadrantControl corControl(Rcpp_corControl); + r = ccaPP(x, y, k, corControl, ppControl, true, fallback, A, B); + } else if(method == "M") { + CorMControl corControl(Rcpp_corControl); + r = ccaPP(x, y, k, corControl, ppControl, true, fallback, A, B); + } else if(method == "pearson") { + CorPearsonControl corControl; + r = ccaPP(x, y, k, corControl, ppControl, false, false, A, B); + } else { + error("method not available"); + } } else if(algorithm == "proj") { // define control object for projections through data points ProjControl ppControl(Rcpp_ppControl); From 5d0e03d0679c6c13dc878c2c3f3f3150834db7bb Mon Sep 17 00:00:00 2001 From: Andreas Alfons Date: Thu, 18 Oct 2012 00:45:54 +0200 Subject: [PATCH 02/14] added function sMaxCorGrid() and commented out ccaGrid() since only k = 1 is implemented yet --- NAMESPACE | 3 +- NEWS | 2 ++ R/cca.R | 82 +++++++++++++++++++++++------------------------ R/maxCor.R | 36 +++++++++++++++++++++ man/ccaGrid.Rd | 14 -------- man/maxCorGrid.Rd | 7 ++++ src/cca.cpp | 5 ++- 7 files changed, 89 insertions(+), 60 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index e4236cd..38e6bc5 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -15,8 +15,7 @@ export(fastMedian) export(maxCorGrid) export(maxCorProj) export(permTest) -export(sCCAgrid) -export(sccaGrid) +export(sMaxCorGrid) import(Rcpp) import(RcppArmadillo) import(parallel) diff --git a/NEWS b/NEWS index d79de46..72a9b79 100644 --- a/NEWS +++ b/NEWS @@ -1,5 +1,7 @@ Changes in ccaPP version 0.2.1 + + New function sMaxCorGrid() for sparse maximum correlation estimators. + + Bugfix in corM(): robust starting values now work with dummy variables. diff --git a/R/cca.R b/R/cca.R index 0af810a..ccd2971 100644 --- a/R/cca.R +++ b/R/cca.R @@ -142,33 +142,32 @@ ccaGrid <- function(x, y, k = 1, cca } -#' @rdname ccaGrid -#' @export -sccaGrid <- function(x, y, k = 1, lambda = 0, - method = c("spearman", "kendall", "quadrant", "M", "pearson"), - control = list(...), nIterations = 10, nAlternate = 10, nGrid = 25, - select = NULL, tol = 1e-06, fallback = FALSE, seed = NULL, ...) { - ## initializations - matchedCall <- match.call() - ## define list of control arguments for algorithm - lambda <- rep(as.numeric(lambda), length.out=2) - lambda[is.na(lambda)] <- formals()$lambda - nIterations <- as.integer(nIterations) - nAlternate <- as.integer(nAlternate) - nGrid <- as.integer(nGrid) - tol <- as.numeric(tol) - ppControl <- list(lambda=lambda, nIterations=nIterations, - nAlternate=nAlternate, nGrid=nGrid, select=select, tol=tol) - ## call workhorse function - cca <- ccaPP(x, y, k, method=method, corControl=control, algorithm="sparse", - ppControl=ppControl, fallback=fallback, seed=seed) - cca$lambda <- lambda - cca$call <- matchedCall - cca -} +# @rdname ccaGrid +# @export + +#sccaGrid <- function(x, y, k = 1, lambda = 0, +# method = c("spearman", "kendall", "quadrant", "M", "pearson"), +# control = list(...), nIterations = 10, nAlternate = 10, nGrid = 25, +# select = NULL, tol = 1e-06, fallback = FALSE, seed = NULL, ...) { +# ## initializations +# matchedCall <- match.call() +# ## define list of control arguments for algorithm +# lambda <- rep(as.numeric(lambda), length.out=2) +# lambda[is.na(lambda)] <- formals()$lambda +# nIterations <- as.integer(nIterations) +# nAlternate <- as.integer(nAlternate) +# nGrid <- as.integer(nGrid) +# tol <- as.numeric(tol) +# ppControl <- list(lambda=lambda, nIterations=nIterations, +# nAlternate=nAlternate, nGrid=nGrid, select=select, tol=tol) +# ## call workhorse function +# cca <- ccaPP(x, y, k, method=method, corControl=control, algorithm="sparse", +# ppControl=ppControl, fallback=fallback, seed=seed) +# cca$call <- matchedCall +# cca +#} -## wrapper functions for more compatibility with package pcaPP #' @rdname ccaGrid #' @export @@ -187,22 +186,23 @@ CCAgrid <- function(x, y, k = 1, cca } -#' @rdname ccaGrid -#' @export -sCCAgrid <- function(x, y, k = 1, lambda = 0, - method = c("spearman", "kendall", "quadrant", "M", "pearson"), - maxiter = 10, maxalter = 10, splitcircle = 25, select=NULL, - zero.tol = 1e-06, fallback = FALSE, seed = NULL, ...) { - ## initializations - matchedCall <- match.call() - ## call ccaGrid() - cca <- sccaGrid(x, y, k=k, lambda=lambda, method=method, nIterations=maxiter, - nAlternate=maxalter, nGrid=splitcircle, select=select, tol=zero.tol, - fallback=fallback, seed=seed, ...) - cca$call <- matchedCall - cca -} +# @rdname ccaGrid +# @export + +#sCCAgrid <- function(x, y, k = 1, lambda = 0, +# method = c("spearman", "kendall", "quadrant", "M", "pearson"), +# maxiter = 10, maxalter = 10, splitcircle = 25, select=NULL, +# zero.tol = 1e-06, fallback = FALSE, seed = NULL, ...) { +# ## initializations +# matchedCall <- match.call() +# ## call ccaGrid() +# cca <- sccaGrid(x, y, k=k, lambda=lambda, method=method, nIterations=maxiter, +# nAlternate=maxalter, nGrid=splitcircle, select=select, tol=zero.tol, +# fallback=fallback, seed=seed, ...) +# cca$call <- matchedCall +# cca +#} #' (Robust) CCA via projections through the data points @@ -302,7 +302,7 @@ ccaProj <- function(x, y, k = 1, cca } -## wrapper function for more compatibility with package pcaPP + #' @rdname ccaProj #' @export diff --git a/R/maxCor.R b/R/maxCor.R index e4f9928..b2eb2fe 100644 --- a/R/maxCor.R +++ b/R/maxCor.R @@ -134,6 +134,42 @@ maxCorGrid <- function(x, y, } +#' @rdname maxCorGrid +#' @export + +sMaxCorGrid <- function(x, y, lambda = 0, + method = c("spearman", "kendall", "quadrant", "M", "pearson"), + control = list(...), nIterations = 10, nAlternate = 10, nGrid = 25, + select = NULL, tol = 1e-06, fallback = FALSE, seed = NULL, ...) { + ## initializations + matchedCall <- match.call() + ## define list of control arguments for algorithm + lambda <- rep(as.numeric(lambda), length.out=2) + lambda[is.na(lambda)] <- formals()$lambda + nIterations <- as.integer(nIterations) + nAlternate <- as.integer(nAlternate) + nGrid <- as.integer(nGrid) + tol <- as.numeric(tol) + ppControl <- list(lambda=lambda, nIterations=nIterations, + nAlternate=nAlternate, nGrid=nGrid, select=select, tol=tol) + ## call workhorse function + maxCor <- maxCorPP(x, y, method=method, corControl=control, + algorithm="sparse", ppControl=ppControl, fallback=fallback, + seed=seed) + ## recompute penalty + a <- maxCor$a + b <- maxCor$b + if(length(a) < 2) lambda[1] <- 0 + if(length(b) < 2) lambda[2] <- 0 + penalty <- lambda[1] * sum(abs(a)) + lambda[2] * sum(abs(b)) + ## add information on penalty and objective function + maxCor$lambda <- lambda + maxCor$objective <- maxCor$cor - penalty + maxCor$call <- matchedCall + maxCor +} + + #' (Robust) maximum correlation via projections through the data points #' #' Compute the maximum correlation between two data sets via projection pursuit diff --git a/man/ccaGrid.Rd b/man/ccaGrid.Rd index 1b72f81..94b112d 100644 --- a/man/ccaGrid.Rd +++ b/man/ccaGrid.Rd @@ -2,8 +2,6 @@ \alias{ccaGrid} \alias{CCAgrid} \alias{print.cca} -\alias{sccaGrid} -\alias{sCCAgrid} \title{(Robust) CCA via alternating series of grid searches} \usage{ ccaGrid(x, y, k = 1, @@ -12,23 +10,11 @@ nGrid = 25, select = NULL, tol = 1e-06, fallback = FALSE, seed = NULL, ...) - sccaGrid(x, y, k = 1, lambda = 0, - method = c("spearman", "kendall", "quadrant", "M", "pearson"), - control = list(...), nIterations = 10, nAlternate = 10, - nGrid = 25, select = NULL, tol = 1e-06, - fallback = FALSE, seed = NULL, ...) - CCAgrid(x, y, k = 1, method = c("spearman", "kendall", "quadrant", "M", "pearson"), maxiter = 10, maxalter = 10, splitcircle = 25, select = NULL, zero.tol = 1e-06, fallback = FALSE, seed = NULL, ...) - - sCCAgrid(x, y, k = 1, lambda = 0, - method = c("spearman", "kendall", "quadrant", "M", "pearson"), - maxiter = 10, maxalter = 10, splitcircle = 25, - select = NULL, zero.tol = 1e-06, fallback = FALSE, - seed = NULL, ...) } \arguments{ \item{x,y}{each can be a numeric vector, matrix or data diff --git a/man/maxCorGrid.Rd b/man/maxCorGrid.Rd index 2d003a2..873fab8 100644 --- a/man/maxCorGrid.Rd +++ b/man/maxCorGrid.Rd @@ -1,6 +1,7 @@ \name{maxCorGrid} \alias{maxCorGrid} \alias{print.maxCor} +\alias{sMaxCorGrid} \title{(Robust) maximum correlation via alternating series of grid searches} \usage{ maxCorGrid(x, y, @@ -8,6 +9,12 @@ control = list(...), nIterations = 10, nAlternate = 10, nGrid = 25, select = NULL, tol = 1e-06, fallback = FALSE, seed = NULL, ...) + + sMaxCorGrid(x, y, lambda = 0, + method = c("spearman", "kendall", "quadrant", "M", "pearson"), + control = list(...), nIterations = 10, nAlternate = 10, + nGrid = 25, select = NULL, tol = 1e-06, + fallback = FALSE, seed = NULL, ...) } \arguments{ \item{x,y}{each can be a numeric vector, matrix or data diff --git a/src/cca.cpp b/src/cca.cpp index a94a807..4ed4059 100644 --- a/src/cca.cpp +++ b/src/cca.cpp @@ -1206,9 +1206,8 @@ SEXP R_ccaPP(SEXP R_x, SEXP R_y, SEXP R_k, SEXP R_method, SEXP R_corControl, error("method not available"); } } else if(algorithm == "sparse") { - // TODO: algorithm currently only works for k = 1 - // larger values of k require backtransformation of weighting vectors - // before computing the penalty + // algorithm currently only works for k = 1; larger values of k require + // backtransformation of weighting vectors before computing the penalty k = 1; // define control object for alternate grid searches SparseGridControl ppControl(Rcpp_ppControl); From 671ad13060a7b3ce3a23e1414768f1293f07534b Mon Sep 17 00:00:00 2001 From: Andreas Alfons Date: Thu, 18 Oct 2012 15:09:08 +0200 Subject: [PATCH 03/14] faster normalization of the weighting vectors in sparse maximum correlation (closes #13) --- DESCRIPTION | 4 ++-- NEWS | 8 +++++-- man/ccaPP-package.Rd | 4 ++-- src/cca.cpp | 56 +++++++++++++++++++++++++++----------------- 4 files changed, 45 insertions(+), 27 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 60f9747..0a9d11e 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: ccaPP Type: Package Title: (Robust) canonical correlation analysis via projection pursuit -Version: 0.2.1 -Date: 2012-10-15 +Version: 0.3.0 +Date: 2012-10-18 Depends: Rcpp (>= 0.9.10), parallel, diff --git a/NEWS b/NEWS index 72a9b79..d89e7cf 100644 --- a/NEWS +++ b/NEWS @@ -1,7 +1,11 @@ -Changes in ccaPP version 0.2.1 +Changes in ccaPP version 0.3.0 + New function sMaxCorGrid() for sparse maximum correlation estimators. - + + + +Changes in ccaPP version 0.2.1 + + Bugfix in corM(): robust starting values now work with dummy variables. diff --git a/man/ccaPP-package.Rd b/man/ccaPP-package.Rd index 672ef69..80f64e7 100644 --- a/man/ccaPP-package.Rd +++ b/man/ccaPP-package.Rd @@ -14,8 +14,8 @@ Canonical correlation analysis and maximum correlation via \tabular{ll}{ Package: \tab ccaPP\cr Type: \tab Package\cr -Version: \tab 0.2.1\cr -Date: \tab 2012-10-15\cr +Version: \tab 0.3.0\cr +Date: \tab 2012-10-18\cr Depends: \tab Rcpp (>= 0.9.10), parallel, pcaPP (>= 1.8-1)\cr Imports: \tab RcppArmadillo (>= 0.2.36)\cr LinkingTo: \tab Rcpp, RcppArmadillo\cr diff --git a/src/cca.cpp b/src/cca.cpp index 4ed4059..33bec07 100644 --- a/src/cca.cpp +++ b/src/cca.cpp @@ -248,6 +248,8 @@ class GridControl { vec&, vec&, bool&); // get equispaced grid of angles in a plane vec getGrid(const uword&); + // compute updated weighting vector in grid search + vec getVector(const vec&, const double&, const uword&); // grid search to update one weighting vector template void gridSearch(const mat&, const uvec&, const vec&, CorControl, vec&, @@ -418,6 +420,14 @@ vec GridControl::getGrid(const uword& i) { return grid; } +// compute updated weighting vector in grid search +vec GridControl::getVector(const vec& a, const double& angle, const uword& j) { + // fast computation of cos(angle) * a + sin(angle) * ej + vec b = cos(angle) * a; + b(j) += sin(angle); + return b; +} + // grid search to update one weighting vector // x ............ data matrix for which to update the weighting vector // orderX ......... order in which to browse through the variables @@ -434,14 +444,13 @@ void GridControl::gridSearch(const mat& x, const uvec& orderX, const vec& y, const uword p = x.n_cols, nGrid = grid.n_elem; // perform grid searches for each canonical basis vector for(uword j = 0; j < p; j++) { - // define current canonical basis vector according to order of columns - vec ej = zeros(p); - ej(orderX(j)) = 1; + // current coordinate to be updated (in the order of columns) + uword orderJ = orderX(j); // perform grid search for the current canonical basis vector vec corY(nGrid); for(uword k = 0; k < nGrid; k++) { double angle = grid(k); - vec currentA = cos(angle) * a + sin(angle) * ej; + vec currentA = getVector(a, grid(k), orderJ); corY(k) = abs(corControl.cor(x * currentA, y)); } // find grid point that maximizes the correlation functional and keep @@ -453,8 +462,7 @@ void GridControl::gridSearch(const mat& x, const uvec& orderX, const vec& y, // of the current grid search may be smaller than the previous maximum if(currentMaxCor > maxCor) { maxCor = currentMaxCor; - double optAngle = grid(whichMax); - a = cos(optAngle) * a + sin(optAngle) * ej; + a = getVector(a, grid(whichMax), orderJ); } } } @@ -614,6 +622,8 @@ class SparseGridControl: public GridControl { // constructors SparseGridControl(); SparseGridControl(List&); + // compute updated weighting vector in grid search + vec getVector(const vec&, const double&, const uword&); // grid search to update one weighting vector template void gridSearch(const mat&, const uvec&, const double&, const vec&, @@ -640,6 +650,17 @@ inline SparseGridControl::SparseGridControl(List& control) lambdaY = lambda[1]; } +// compute updated weighting vector in grid search +vec SparseGridControl::getVector(const vec& a, + const double& angle, const uword& j) { + // fast computation of b / ||b|| with b = cos(angle) * a + sin(angle) * ej + double tanAngle = tan(angle); + double denominator = sqrt(1 + 2*tanAngle*a(j) + tanAngle*tanAngle); + vec b = a / denominator; + b(j) += tanAngle / denominator; + return b; +} + // sparse grid search to update one weighting vector // x .............. data matrix for which to update the weighting vector // orderX ......... order in which to browse through the variables @@ -658,15 +679,12 @@ void SparseGridControl::gridSearch(const mat& x, const uvec& orderX, const uword p = x.n_cols, nGrid = grid.n_elem; // perform grid searches for each canonical basis vector for(uword j = 0; j < p; j++) { - // define current canonical basis vector according to order of columns - vec ej = zeros(p); - ej(orderX(j)) = 1; + // current coordinate to be updated (in the order of columns) + uword orderJ = orderX(j); // perform grid search for the current canonical basis vector vec objective(nGrid); for(uword k = 0; k < nGrid; k++) { - double angle = grid(k); - vec currentA = cos(angle) * a + sin(angle) * ej; - currentA = currentA / norm(currentA, 2); + vec currentA = getVector(a, grid(k), orderJ); double currentCorY = abs(corControl.cor(x * currentA, y)); objective(k) = currentCorY - lambda * norm(currentA, 1); } @@ -680,8 +698,7 @@ void SparseGridControl::gridSearch(const mat& x, const uvec& orderX, // previous maximum if(currentMaxObjective > maxObjective) { maxObjective = currentMaxObjective; - double optAngle = grid(whichMax); - a = cos(optAngle) * a + sin(optAngle) * ej; + a = getVector(a, grid(whichMax), orderJ); } } } @@ -708,15 +725,13 @@ void SparseGridControl::gridSearch(const mat& x, const uvec& orderX, const uword p = x.n_cols, nGrid = grid.n_elem; // perform grid searches for each canonical basis vector for(uword j = 0; j < p; j++) { - // define current canonical basis vector according to order of columns - vec ej = zeros(p); - ej(orderX(j)) = 1; + // current coordinate to be updated (in the order of columns) + uword orderJ = orderX(j); // perform grid search for the current canonical basis vector vec corY(nGrid), objective(nGrid), penalties(nGrid); for(uword k = 0; k < nGrid; k++) { double angle = grid(k); - vec currentA = cos(angle) * a + sin(angle) * ej; - currentA = currentA / norm(currentA, 2); + vec currentA = getVector(a, grid(k), orderJ); corY(k) = abs(corControl.cor(x * currentA, y)); penalties(k) = lambda * norm(currentA, 1); objective(k) = corY(k) - penalties(k) - penaltyY; @@ -731,8 +746,7 @@ void SparseGridControl::gridSearch(const mat& x, const uvec& orderX, // previous maximum if(currentMaxObjective > maxObjective) { maxCor = corY(whichMax); - double optAngle = grid(whichMax); - a = cos(optAngle) * a + sin(optAngle) * ej; + a = getVector(a, grid(whichMax), orderJ); penaltyX = penalties(whichMax); maxObjective = currentMaxObjective; } From 903166a77a2c688607362ea487a08f5b1b03f163 Mon Sep 17 00:00:00 2001 From: Andreas Alfons Date: Fri, 23 Nov 2012 19:38:36 +0100 Subject: [PATCH 04/14] added code for multiple values of the penalty parameter for one multivariate and one univariate data set --- .project | 19 --- R/cca.R | 6 +- R/maxCor.R | 127 +++++++++++---- src/cca.cpp | 437 +++++++++++++++++++++++++++++++++------------------- src/cca.h | 4 +- 5 files changed, 379 insertions(+), 214 deletions(-) delete mode 100644 .project diff --git a/.project b/.project deleted file mode 100644 index a807d01..0000000 --- a/.project +++ /dev/null @@ -1,19 +0,0 @@ - - - ccaPP - - - - - - de.walware.statet.r.builders.RSupport - - - - - - de.walware.statet.base.StatetNature - de.walware.statet.r.RNature - de.walware.statet.r.RPkgNature - - diff --git a/R/cca.R b/R/cca.R index ccd2971..ad01bd3 100644 --- a/R/cca.R +++ b/R/cca.R @@ -322,8 +322,7 @@ CCAproj <- function(x, y, k = 1, ## workhorse function ccaPP <- function(x, y, k = 1, method = c("spearman", "kendall", "quadrant", "M", "pearson"), - corControl, forceConsistency = TRUE, - algorithm = c("grid", "sparse", "proj"), + corControl, forceConsistency = TRUE, algorithm = c("grid", "proj"), ppControl, fallback = FALSE, seed = NULL) { ## initializations x <- as.matrix(x) @@ -348,7 +347,7 @@ ccaPP <- function(x, y, k = 1, method <- match.arg(method) corControl <- getCorControl(method, corControl, forceConsistency) # additional checks for grid search algorithm - if(algorithm != "proj") { + if(algorithm == "grid") { # check subset of variables to be used for determining the order of # the variables from the respective other data set select <- ppControl$select @@ -389,7 +388,6 @@ ccaPP <- function(x, y, k = 1, R_corControl=corControl, R_algorithm=algorithm, R_ppControl=ppControl, R_fallback=isTRUE(fallback), PACKAGE="ccaPP") - cca$cor <- drop(cca$cor) } ## assign class and return results class(cca) <- "cca" diff --git a/R/maxCor.R b/R/maxCor.R index b2eb2fe..cd99fc7 100644 --- a/R/maxCor.R +++ b/R/maxCor.R @@ -137,36 +137,32 @@ maxCorGrid <- function(x, y, #' @rdname maxCorGrid #' @export -sMaxCorGrid <- function(x, y, lambda = 0, - method = c("spearman", "kendall", "quadrant", "M", "pearson"), - control = list(...), nIterations = 10, nAlternate = 10, nGrid = 25, - select = NULL, tol = 1e-06, fallback = FALSE, seed = NULL, ...) { - ## initializations - matchedCall <- match.call() - ## define list of control arguments for algorithm - lambda <- rep(as.numeric(lambda), length.out=2) - lambda[is.na(lambda)] <- formals()$lambda - nIterations <- as.integer(nIterations) - nAlternate <- as.integer(nAlternate) - nGrid <- as.integer(nGrid) - tol <- as.numeric(tol) - ppControl <- list(lambda=lambda, nIterations=nIterations, - nAlternate=nAlternate, nGrid=nGrid, select=select, tol=tol) - ## call workhorse function - maxCor <- maxCorPP(x, y, method=method, corControl=control, - algorithm="sparse", ppControl=ppControl, fallback=fallback, - seed=seed) - ## recompute penalty - a <- maxCor$a - b <- maxCor$b - if(length(a) < 2) lambda[1] <- 0 - if(length(b) < 2) lambda[2] <- 0 - penalty <- lambda[1] * sum(abs(a)) + lambda[2] * sum(abs(b)) - ## add information on penalty and objective function - maxCor$lambda <- lambda - maxCor$objective <- maxCor$cor - penalty - maxCor$call <- matchedCall - maxCor +sMaxCorGrid <- function(x, y, lambdaX = 0, lambdaY = 0, + method = c("spearman", "kendall", "quadrant", + "M", "pearson"), + control = list(...), nIterations = 10, nAlternate = 10, + nGrid = 25, select = NULL, tol = 1e-06, + fallback = FALSE, seed = NULL, ...) { + ## initializations + matchedCall <- match.call() + ## define list of control arguments for algorithm + lambdaX <- as.numeric(lambdaX) + lambdaX[is.na(lambdaX)] <- formals()$lambdaX + lambdaX <- sort(unique(lambdaX)) + lambdaY <- as.numeric(lambdaY) + lambdaY[is.na(lambdaY)] <- formals()$lambdaY + lambdaY <- sort(unique(lambdaY)) + nIterations <- as.integer(nIterations) + nAlternate <- as.integer(nAlternate) + nGrid <- as.integer(nGrid) + tol <- as.numeric(tol) + ppControl <- list(lambdaX=lambdaX, lambdaY=lambdaY, nIterations=nIterations, + nAlternate=nAlternate, nGrid=nGrid, select=select, tol=tol) + ## call workhorse function + maxCor <- sMaxCorPP(x, y, method=method, corControl=control, algorithm="grid", + ppControl=ppControl, fallback=fallback, seed=seed) + maxCor$call <- matchedCall + maxCor } @@ -260,7 +256,7 @@ maxCorProj <- function(x, y, } -## workhorse function +## workhorse function for maximum correlation maxCorPP <- function(x, y, ...) { ## call workhorse function for canonical correlation analysis maxCor <- ccaPP(x, y, forceConsistency=FALSE, ...) @@ -269,3 +265,72 @@ maxCorPP <- function(x, y, ...) { class(maxCor) <- "maxCor" maxCor } + + +## workhorse function for sparse maximum correlation +sMaxCorPP <- function(x, y, method = c("spearman", "kendall", "quadrant", + "M", "pearson"), + corControl, algorithm = "grid", ppControl, + fallback = FALSE, seed = NULL) { + ## initializations + x <- as.matrix(x) + y <- as.matrix(y) + n <- nrow(x) + if(nrow(y) != n) stop("'x' and 'y' must have the same number of observations") + p <- ncol(x) + q <- ncol(y) + ## prepare the data and call C++ function + if(n == 0 || p == 0 || q == 0) { + # zero dimension + maxCor <- list(cor=NA, a=numeric(), b=numeric(), objective=NA) + } else { + # check method and get list of control arguments + method <- match.arg(method) + corControl <- getCorControl(method, corControl, forceConsistency) + # additional checks for grid search algorithm + if(algorithm == "grid") { + # check subset of variables to be used for determining the order of + # the variables from the respective other data set + select <- ppControl$select + ppControl$select <- NULL + if(!is.null(select)) { + if(is.list(select)) { + # make sure select is a list with two index vectors and + # drop invalid indices from each vector + select <- rep(select, length.out=2) + select <- mapply(function(indices, max) { + indices <- as.integer(indices) + indices[which(indices > 0 & indices <= max)] - 1 + }, select, c(p, q)) + valid <- sapply(select, length) > 0 + # add the two index vectors to control object + if(all(valid)) { + ppControl$selectX <- select[[1]] + ppControl$selectY <- select[[2]] + } else select <- NULL + } else { + # check number of indices to sample + select <- rep(as.integer(select), length.out=2) + valid <- !is.na(select) & select > 0 & select < c(p, q) + if(all(valid)) { + # generate index vectors and add them to control object + if(!is.null(seed)) set.seed(seed) + ppControl$selectX <- sample.int(p, select[1]) - 1 + ppControl$selectY <- sample.int(q, select[2]) - 1 + } else select <- NULL + } + } + if(is.null(select)) { + ppControl$selectX <- ppControl$selectY <- integer() + } + } + # call C++ function + maxCor <- .Call("R_sMaxCorPP", R_x=x, R_y=y, R_method=method, + R_corControl=corControl, R_algorithm=algorithm, + R_ppControl=ppControl, R_fallback=isTRUE(fallback), + PACKAGE="ccaPP") + } + ## assign class and return results + class(maxCor) <- "maxCor" + maxCor +} diff --git a/src/cca.cpp b/src/cca.cpp index 33bec07..43d38d6 100644 --- a/src/cca.cpp +++ b/src/cca.cpp @@ -173,7 +173,7 @@ double CorMControl::cor(const vec& x, const vec& y) { // for testing purposes template -double fastCor(const vec& x, const vec& y, CorControl control) { +double fastCor(const vec& x, const vec& y, CorControl& control) { return control.cor(x, y); } @@ -241,10 +241,10 @@ class GridControl { GridControl(List&); // get the order in which to update the elements of a weighting vector template - void findOrder(const mat&, const vec&, CorControl, uvec&, double&, vec&); + void findOrder(const mat&, const vec&, CorControl&, uvec&, double&, vec&); // get the order in which to update the elements of the weighting vectors template - void findOrder(const mat&, const mat&, CorControl, uvec&, uvec&, double&, + void findOrder(const mat&, const mat&, CorControl&, uvec&, uvec&, double&, vec&, vec&, bool&); // get equispaced grid of angles in a plane vec getGrid(const uword&); @@ -252,13 +252,13 @@ class GridControl { vec getVector(const vec&, const double&, const uword&); // grid search to update one weighting vector template - void gridSearch(const mat&, const uvec&, const vec&, CorControl, vec&, + void gridSearch(const mat&, const uvec&, const vec&, CorControl&, vec&, double&, vec&); // set counter of consecutive iterations without improvement void setCounter(uword&, const double&, const double&); // maximum correlation between multivariate data sets x and y template - double maxCor(const mat&, const mat&, CorControl, vec&, vec&); + double maxCor(const mat&, const mat&, CorControl&, vec&, vec&); }; // constructors @@ -297,7 +297,7 @@ inline GridControl::GridControl(List& control) { // maxCor ....... maximum correlation to get initial value // a ............ weighting vector to get initial value template -void GridControl::findOrder(const mat& x, const vec& y, CorControl corControl, +void GridControl::findOrder(const mat& x, const vec& y, CorControl& corControl, uvec& orderX, double& maxCor, vec& a) { // compute columnwise absolute correlations of x with y const uword p = x.n_cols; @@ -327,7 +327,7 @@ void GridControl::findOrder(const mat& x, const vec& y, CorControl corControl, // startWithX ... logical to be computed that indicates whether to start with // the first data set in the alternate grid searches template -void GridControl::findOrder(const mat& x, const mat& y, CorControl corControl, +void GridControl::findOrder(const mat& x, const mat& y, CorControl& corControl, uvec& orderX, uvec& orderY, double& maxCor, vec& a, vec& b, bool& startWithX) { const uword p = x.n_cols, q = y.n_cols; @@ -439,7 +439,7 @@ vec GridControl::getVector(const vec& a, const double& angle, const uword& j) { // a ............ weighting vector to be updated template void GridControl::gridSearch(const mat& x, const uvec& orderX, const vec& y, - CorControl corControl, vec& grid, double& maxCor, vec& a) { + CorControl& corControl, vec& grid, double& maxCor, vec& a) { // initializations const uword p = x.n_cols, nGrid = grid.n_elem; // perform grid searches for each canonical basis vector @@ -488,7 +488,7 @@ void GridControl::setCounter(uword& counter, const double& current, // a ............. first weighting vector to be updated // b ............. second weighting vector to be updated template -double GridControl::maxCor(const mat& x, const mat& y, CorControl corControl, +double GridControl::maxCor(const mat& x, const mat& y, CorControl& corControl, vec& a, vec& b) { // initializations uword p = x.n_cols, q = y.n_cols; @@ -618,36 +618,45 @@ double GridControl::maxCor(const mat& x, const mat& y, CorControl corControl, class SparseGridControl: public GridControl { public: - double lambdaX, lambdaY; // penalty parameters + vec lambdaX, lambdaY; // penalty parameters // constructors SparseGridControl(); SparseGridControl(List&); // compute updated weighting vector in grid search vec getVector(const vec&, const double&, const uword&); - // grid search to update one weighting vector + // grid search to update weighting vector for multivariate x and univariate y template void gridSearch(const mat&, const uvec&, const double&, const vec&, - CorControl, vec&, double&, vec&); + CorControl&, vec&, vec&, double&); // grid search to update one weighting vector for multivariate x and y template void gridSearch(const mat&, const uvec&, const double&, const vec&, - const double&, CorControl, vec&, double&, vec&, double&, double&); - // maximum correlation between multivariate data sets x and y + const double&, CorControl&, vec&, double&, vec&, double&, double&); + // workhorse function for maximum correlation between multivariate x and + // univariate y + template + void maxCorFit(const mat&, const uvec&, const double&, const vec&, + CorControl&, double&, vec&, double&); + // workhorse function for maximum correlation between multivariate x and y + template + void maxCorFit(const mat&, const uvec&, const double&, const mat&, + const uvec&, const double&, CorControl&, double&, vec&, vec&, + double&, double&, double&); + // maximum correlation between multivariate data sets x and y template - double maxCor(const mat&, const mat&, CorControl, vec&, vec&); + vec maxCor(const mat&, const mat&, CorControl&, mat&, mat&, vec&); }; // constructors inline SparseGridControl::SparseGridControl() { - lambdaX = 0; - lambdaY = 0; + lambdaX = zeros(1); + lambdaY = zeros(1); } inline SparseGridControl::SparseGridControl(List& control) : GridControl(control) { - SEXP R_lambda = control["lambda"]; - NumericVector lambda(R_lambda); - lambdaX = lambda[0]; - lambdaY = lambda[1]; + SEXP R_lambdaX = control["lambdaX"], R_lambdaY = control["lambdaY"]; + lambdaX = as(R_lambdaX); + lambdaY = as(R_lambdaY); } // compute updated weighting vector in grid search @@ -661,20 +670,20 @@ vec SparseGridControl::getVector(const vec& a, return b; } -// sparse grid search to update one weighting vector +// sparse grid search to update weighting vector for multivariate x and +// univariate y // x .............. data matrix for which to update the weighting vector // orderX ......... order in which to browse through the variables // lambda ......... penalty parameter for the weighting vector to be updated -// y .............. linear combination of the other data matrix according to -// the other weighting vector, which is kept fixed +// y .............. other data vector // corControl ..... control object to compute correlation // grid ........... grid points to be used for grid search -// maxObjective ... maximum correlation to be updated // a .............. weighting vector to be updated +// maxObjective ... maximum value of the objective function to be updated template void SparseGridControl::gridSearch(const mat& x, const uvec& orderX, - const double& lambda, const vec& y, CorControl corControl, - vec& grid, double& maxObjective, vec& a) { + const double& lambda, const vec& y, CorControl& corControl, + vec& grid, vec& a, double& maxObjective) { // initializations const uword p = x.n_cols, nGrid = grid.n_elem; // perform grid searches for each canonical basis vector @@ -713,13 +722,13 @@ void SparseGridControl::gridSearch(const mat& x, const uvec& orderX, // corControl ..... control object to compute correlation // grid ........... grid points to be used for grid search // maxCor ......... maximum correlation to be updated -// maxObjective ... maximum value of the objective function to be updated // a .............. weighting vector to be updated // penaltyX ....... penalty for the weighting vector to be updated +// maxObjective ... maximum value of the objective function to be updated template void SparseGridControl::gridSearch(const mat& x, const uvec& orderX, const double& lambda, const vec& y, const double& penaltyY, - CorControl corControl, vec& grid, double& maxCor, vec& a, + CorControl& corControl, vec& grid, double& maxCor, vec& a, double& penaltyX, double& maxObjective) { // initializations const uword p = x.n_cols, nGrid = grid.n_elem; @@ -753,6 +762,69 @@ void SparseGridControl::gridSearch(const mat& x, const uvec& orderX, } } +// workhorse function for maximum correlation between multivariate x and +// univariate y +template +void SparseGridControl::maxCorFit(const mat& x, const uvec& orderX, + const double& lambda, const vec& y, CorControl& corControl, + double& maxCor, vec& a, double& maxObjective) { + // initializations + double previousMaxObjective = R_NegInf; + // stop if there are two consecutive iterations without improvement + uword i = 0, convCounter = 0; + while((i < nIterations) && (convCounter < 2)) { + previousMaxObjective = maxObjective; + vec grid = getGrid(i+1); // define vector of grid points + gridSearch(x, orderX, lambda, y, corControl, grid, a, maxObjective); + i++; + setCounter(convCounter, maxObjective, previousMaxObjective); + } + // check direction + maxCor = corControl.cor(x * a, y); + if(maxCor < 0) { + maxCor = -maxCor; + a = -a; + } +} + +// workhorse function for maximum correlation between multivariate x and y +template +void SparseGridControl::maxCorFit(const mat& x, const uvec& orderX, + const double& lambda1, const mat& y, const uvec& orderY, + const double& lambda2, CorControl& corControl, double& maxCor, vec& a, + vec& b, double& penaltyX, double& penaltyY, double& maxObjective) { + // initializations + double previousMaxObjective = R_NegInf; + // stop if there are two consecutive iterations without improvement + uword i = 0, convCounter = 0; + while((i < nIterations) && (convCounter < 2)) { + previousMaxObjective = maxObjective; + vec grid = getGrid(i+1); // define vector of grid points + uword j = 0; + double altMaxObjective = R_NegInf; + while((j < nAlternate) && ((maxObjective - altMaxObjective) > tol)) { + altMaxObjective = maxObjective; + // maximize correlation functional over a keeping b fixed + vec yb = y * b; // linear combination of columns of y + gridSearch(x, orderX, lambda1, yb, penaltyY, corControl, + grid, maxCor, a, penaltyX, maxObjective); + // maximize correlation functional over b keeping a fixed + vec xa = x * a; // linear combination of columns of x + gridSearch(y, orderY, lambda2, xa, penaltyX, corControl, + grid, maxCor, b, penaltyY, maxObjective); + j++; + } + i++; + setCounter(convCounter, maxObjective, previousMaxObjective); + } + // check direction + maxCor = corControl.cor(x * a, y * b); + if(maxCor < 0) { + maxCor = -maxCor; + b = -b; + } +} + // maximum correlation between multivariate data sets x and y based on // sparse alternate grid searches // x ............. first data matrix @@ -761,135 +833,116 @@ void SparseGridControl::gridSearch(const mat& x, const uvec& orderX, // a ............. first weighting vector to be updated // b ............. second weighting vector to be updated template -double SparseGridControl::maxCor(const mat& x, const mat& y, - CorControl corControl, vec& a, vec& b) { +vec SparseGridControl::maxCor(const mat& x, const mat& y, + CorControl& corControl, mat& A, mat& B, vec& objective) { // initializations uword p = x.n_cols, q = y.n_cols; + vec r; // initialize vector of maximum correlations // perform alternate grid searches if both data sets are multivariate // if one data set is univariate, alternate grid searches are not necessary - double maxCor; // initialize maximum correlation if((p == 1) && (q == 1)) { // both data sets are univariate - a.ones(p); b.ones(q); + r.set_size(1); objective.set_size(1); + A.ones(p, 1); B.ones(q, 1); + lambdaX.zeros(p); lambdaY.zeros(q); vec xx = x.unsafe_col(0), yy = y.unsafe_col(0); // reuse memory - maxCor = corControl.cor(xx, yy); // compute correlation + // compute correlation + r(0) = corControl.cor(xx, yy); // check sign of correlation - if(maxCor < 0) { - maxCor = -maxCor; - b = -b; + if(r(0) < 0) { + r(0) = -r(0); + B(0, 0) = -B(0, 0); } + objective(0) = r(0); } else { - double maxObjective, previousMaxObjective = R_NegInf; + double maxCor, maxObjective; if((p > 1) && (q == 1)) { // x is multivariate, y is univariate - vec yy = y.unsafe_col(0); // reuse memory - uvec orderX(p); - a.zeros(p); b.ones(q); + uword nLambdaX = lambdaX.n_elem; + r.set_size(nLambdaX); objective.set_size(nLambdaX); + A.set_size(p, nLambdaX); B.ones(q, nLambdaX); + lambdaY.zeros(q); + vec yy = y.unsafe_col(0); // reuse memory + // find order of x variables + uvec orderX(p); + vec a = zeros(p); findOrder(x, yy, corControl, orderX, maxCor, a); // column order - maxObjective = maxCor - lambdaX; // L1 norm of a is 1 - // stop if there are two consecutive iterations without improvement - uword i = 0, convCounter = 0; - while((i < nIterations) && (convCounter < 2)) { - previousMaxObjective = maxObjective; - vec grid = getGrid(i+1); // define vector of grid points - gridSearch(x, orderX, lambdaX, yy, corControl, - grid, maxObjective, a); // grid search - i++; - setCounter(convCounter, maxObjective, previousMaxObjective); - } + // compute maximum correlation for each penalty parameter + for(int k = nLambdaX-1; k >= 0; k--) { + maxObjective = maxCor - lambdaX(k) * norm(a, 1); + maxCorFit(x, orderX, lambdaX(k), yy, corControl, + maxCor, a, maxObjective); + r(k) = maxCor; + A.col(k) = a; + objective(k) = maxObjective; + } } else if((p == 1) && (q > 1)) { // x is univariate, y is multivariate - vec xx = x.unsafe_col(0); // reuse memory - uvec orderY(q); - a.ones(p); b.zeros(q); - findOrder(y, xx, corControl, orderY, maxCor, b); // column order - maxObjective = maxCor - lambdaY; // L1 norm of b is 1 - // stop if there are two consecutive iterations without improvement - uword i = 0, convCounter = 0; - while((i < nIterations) && (convCounter < 2)) { - previousMaxObjective = maxObjective; - vec grid = getGrid(i+1); // define vector of grid points - gridSearch(y, orderY, lambdaY, xx, corControl, - grid, maxObjective, b); // grid search - i++; - setCounter(convCounter, maxObjective, previousMaxObjective); - } + uword nLambdaY = lambdaY.n_elem; + r.set_size(nLambdaY); objective.set_size(nLambdaY); + A.ones(p, nLambdaY); B.set_size(q, nLambdaY); + lambdaX.zeros(p); + vec xx = x.unsafe_col(0); // reuse memory + // find order of y variables + uvec orderY(q); + vec b = zeros(q); + findOrder(y, xx, corControl, orderY, maxCor, b); // column order + // compute maximum correlation for each penalty parameter + for(int k = nLambdaY-1; k >= 0; k--) { + maxObjective = maxCor - lambdaY(k) * norm(b, 1); + maxCorFit(y, orderY, lambdaY(k), xx, corControl, + maxCor, b, maxObjective); + r(k) = maxCor; + B.col(k) = b; + objective(k) = maxObjective; + } } else if((p > 1) && (q > 1)) { // both data sets are multivariate - // compute correlations between variables in x and y - uvec orderX(p), orderY(q); - a.zeros(p); b.zeros(q); - bool startWithX; - findOrder(x, y, corControl, orderX, orderY, - maxCor, a, b, startWithX); // column orders - double penaltyX = lambdaX, penaltyY = lambdaY; // L1 norms are 1 - maxObjective = maxCor - penaltyX - penaltyY; - // perform alternate grid searches - if(startWithX) { - // start with grid search for x - // stop if there are two consecutive iterations without improvement - uword i = 0, convCounter = 0; - while((i < nIterations) && (convCounter < 2)) { - previousMaxObjective = maxObjective; - vec grid = getGrid(i+1); // define vector of grid points - uword j = 0; - double altMaxObjective = R_NegInf; - while((j < nAlternate) && ((maxObjective - altMaxObjective) > tol)) { - altMaxObjective = maxObjective; - // maximize correlation functional over a keeping b fixed - vec yb = y * b; // linear combination of columns of y - gridSearch(x, orderX, lambdaX, yb, penaltyY, corControl, - grid, maxCor, a, penaltyX, maxObjective); - // maximize correlation functional over b keeping a fixed - vec xa = x * a; // linear combination of columns of x - gridSearch(y, orderY, lambdaY, xa, penaltyX, corControl, - grid, maxCor, b, penaltyY, maxObjective); - j++; - } - i++; - setCounter(convCounter, maxObjective, previousMaxObjective); - } - } else { - // start with grid search for y - // stop if there are two consecutive iterations without improvement - uword i = 0, convCounter = 0; - while((i < nIterations) && (convCounter < 2)) { - previousMaxObjective = maxObjective; - vec grid = getGrid(i+1); // define vector of grid points - uword j = 0; - double altMaxObjective = R_NegInf; - while((j < nAlternate) && ((maxObjective - altMaxObjective) > tol)) { - altMaxObjective = maxObjective; - // maximize correlation functional over b keeping a fixed - vec xa = x * a; // linear combination of columns of x - gridSearch(y, orderY, lambdaY, xa, penaltyX, corControl, - grid, maxCor, b, penaltyY, maxObjective); - // maximize correlation functional over a keeping b fixed - vec yb = y * b; // linear combination of columns of y - gridSearch(x, orderX, lambdaX, yb, penaltyY, corControl, - grid, maxCor, a, penaltyX, maxObjective); - j++; - } - i++; - setCounter(convCounter, maxObjective, previousMaxObjective); - } - } - } else { - return NA_REAL; // should never happen - } - // check direction - maxCor = corControl.cor(x * a, y * b); - if(maxCor < 0) { - maxCor = -maxCor; - if((p > 1) && (q == 1)) { - a = -a; - } else { - b = -b; - } +// uword nLambdaX = lambdaX.n_elem, nLambdaY = lambdaY.n_elem; +// uword nLambda = nLambdaX * nLambdaY; +// r.set_size(nLambda); objective.set_size(nLambda); +// A.ones(p, nLambda); B.set_size(q, nLambda); +// // find order of x and y variables +// uvec orderX(p), orderY(q); +// vec a = zeros(p), b = zeros(q); +// bool startWithX; +// findOrder(x, y, corControl, orderX, orderY, maxCor, a, b, startWithX); +// // compute initial value of objective function (L1 norms are 1) +// // compute maximum correlation for each combination of penalty parameters +// for(int k = nLambdaX-1; k >= 0; k--) { +// for(int l = nLambdaY-1; l >= 0; l--) { +// double penaltyX = lambdaX(k) * norm(a, 1); +// double penaltyY = lambdaY(l) * norm(b, 1); +// maxObjective = maxCor - penaltyX - penaltyY; +// Rprintf("lambdaX = %f, lambdaY=%f\n", lambdaX(k), lambdaY(l)); +// Rprintf(" before: maxCor = %f, maxObjective = %f\n", maxCor, maxObjective); +// // perform alternate grid searches +// if(startWithX) { +// // start with grid search for x +// maxCorFit(x, orderX, lambdaX(k), y, orderY, lambdaY(l), corControl, +// maxCor, a, b, penaltyX, penaltyY, maxObjective); +// } else { +// // start with grid search for y +// maxCorFit(y, orderY, lambdaY(l), x, orderX, lambdaX(k), corControl, +// maxCor, b, a, penaltyY, penaltyX, maxObjective); +// } +// Rprintf(" after: maxCor = %f, maxObjective = %f\n", maxCor, maxObjective); +// uword kl = k * nLambdaX + l; +// r(kl) = maxCor; +// A.col(kl) = a; +// B.col(kl) = b; +// objective(kl) = maxObjective; +// } +// // reset values +// uword kl = k * nLambdaX + nLambdaY - 1; +// maxCor = r(kl); +// a = A.col(kl); +// b = B.col(kl); +// } } } // return maximum correlation - return maxCor; + return r; } @@ -907,7 +960,7 @@ class ProjControl { mat getDirections(const mat&); // maximum correlation between multivariate data sets x and y template - double maxCor(const mat&, const mat&, CorControl, vec&, vec&); + double maxCor(const mat&, const mat&, CorControl&, vec&, vec&); }; // constructors @@ -949,7 +1002,7 @@ mat ProjControl::getDirections(const mat& x) { // a ............ first weighting vector to be updated // b ............ second weighting vector to be updated template -double ProjControl::maxCor(const mat& x, const mat& y, CorControl corControl, +double ProjControl::maxCor(const mat& x, const mat& y, CorControl& corControl, vec& a, vec& b) { // initializations uword n = x.n_rows, p = x.n_cols, q = y.n_cols; @@ -1114,6 +1167,7 @@ mat householder(const vec& a) { return P; } + // canonical correlation analysis via projection pursuit // x ............ first data matrix // y ............ second data matrix @@ -1125,8 +1179,8 @@ mat householder(const vec& a) { // A ............ matrix of canonical vectors for first matrix to be updated // B ............ matrix of canonical vectors for second matrix to be updated template -vec ccaPP(const mat& x, const mat& y, const uword& k, CorControl corControl, - PPControl ppControl, const bool& robust, const bool& fallback, +vec ccaPP(const mat& x, const mat& y, const uword& k, CorControl& corControl, + PPControl& ppControl, const bool& robust, const bool& fallback, mat& A, mat& B) { // initializations uword p = x.n_cols, q = y.n_cols; @@ -1192,7 +1246,7 @@ SEXP R_ccaPP(SEXP R_x, SEXP R_y, SEXP R_k, SEXP R_method, SEXP R_corControl, string method = as(R_method); // convert character string List Rcpp_corControl(R_corControl); // list of control parameters string algorithm = as(R_algorithm); // convert character string - List Rcpp_ppControl(R_ppControl); // list of control parameters + List Rcpp_ppControl(R_ppControl); // list of control parameters bool fallback = as(R_fallback); // convert to boolean // initialize results vec r; @@ -1219,12 +1273,9 @@ SEXP R_ccaPP(SEXP R_x, SEXP R_y, SEXP R_k, SEXP R_method, SEXP R_corControl, } else { error("method not available"); } - } else if(algorithm == "sparse") { - // algorithm currently only works for k = 1; larger values of k require - // backtransformation of weighting vectors before computing the penalty - k = 1; - // define control object for alternate grid searches - SparseGridControl ppControl(Rcpp_ppControl); + } else if(algorithm == "proj") { + // define control object for projections through data points + ProjControl ppControl(Rcpp_ppControl); // define control object for the correlations and call the arma version if(method == "spearman") { CorSpearmanControl corControl(Rcpp_corControl); @@ -1244,35 +1295,103 @@ SEXP R_ccaPP(SEXP R_x, SEXP R_y, SEXP R_k, SEXP R_method, SEXP R_corControl, } else { error("method not available"); } - } else if(algorithm == "proj") { - // define control object for projections through data points - ProjControl ppControl(Rcpp_ppControl); + } else { + error("algorithm not available"); + } + // wrap and return result + return List::create( + Named("cor") = wrap(r.begin(), r.end()), + Named("A") = A, + Named("B") = B + ); +} + + +// sparse maximum correlation via projection pursuit +// x ............ first data matrix +// y ............ second data matrix +// corControl ... control object to compute correlation +// ppControl .... control object for sparse algorithm +// robust ....... should the data be robustly standardized? +// fallback ..... should the fallback mode for robust standardization be used? +// A ............ matrix of weighting vectors for first matrix to be updated +// B ............ matrix of weighting vectors for second matrix to be updated +template +vec sMaxCorPP(const mat& x, const mat& y, CorControl& corControl, + PPControl& ppControl, const bool& robust, const bool& fallback, + mat& A, mat& B, vec& objective) { + // initializations + const uword p = x.n_cols, q = y.n_cols; + // standardize the data + vec scaleX, scaleY; + mat xs = standardize(x, robust, fallback, scaleX); + mat ys = standardize(y, robust, fallback, scaleY); + // compute maximum correlations with standardized data + vec r = ppControl.maxCor(xs, ys, corControl, A, B, objective); + // transform weighting vectors back to original scale + for(uword k = 0; k < A.n_cols; k++) { + vec a = A.unsafe_col(k); backtransform(a, scaleX); + } + for(uword k = 0; k < B.n_cols; k++) { + vec b = B.unsafe_col(k); backtransform(b, scaleY); + } + // return vector of maximum correlations + return r; +} + +// R interface +SEXP R_sMaxCorPP(SEXP R_x, SEXP R_y, SEXP R_method, SEXP R_corControl, + SEXP R_algorithm, SEXP R_ppControl, SEXP R_fallback) { + // initializations + NumericMatrix Rcpp_x(R_x), Rcpp_y(R_y); // convert data to Rcpp types + mat x(Rcpp_x.begin(), Rcpp_x.nrow(), Rcpp_x.ncol(), false); // convert data + mat y(Rcpp_y.begin(), Rcpp_y.nrow(), Rcpp_y.ncol(), false); // to arma types + string method = as(R_method); // convert character string + List Rcpp_corControl(R_corControl); // list of control parameters + string algorithm = as(R_algorithm); // convert character string + List Rcpp_ppControl(R_ppControl); // list of control parameters + bool fallback = as(R_fallback); // convert to boolean + // initialize results + vec r, lambdaX, lambdaY, objective; + mat A, B; + if(algorithm == "grid") { + // define control object for alternate grid searches + SparseGridControl ppControl(Rcpp_ppControl); // define control object for the correlations and call the arma version if(method == "spearman") { CorSpearmanControl corControl(Rcpp_corControl); - r = ccaPP(x, y, k, corControl, ppControl, true, fallback, A, B); + r = sMaxCorPP(x, y, corControl, ppControl, true, fallback, + A, B, objective); } else if(method == "kendall") { CorKendallControl corControl(Rcpp_corControl); - r = ccaPP(x, y, k, corControl, ppControl, true, fallback, A, B); + r = sMaxCorPP(x, y, corControl, ppControl, true, fallback, + A, B, objective); } else if(method == "quadrant") { CorQuadrantControl corControl(Rcpp_corControl); - r = ccaPP(x, y, k, corControl, ppControl, true, fallback, A, B); + r = sMaxCorPP(x, y, corControl, ppControl, true, fallback, + A, B, objective); } else if(method == "M") { CorMControl corControl(Rcpp_corControl); - r = ccaPP(x, y, k, corControl, ppControl, true, fallback, A, B); + r = sMaxCorPP(x, y, corControl, ppControl, true, fallback, + A, B, objective); } else if(method == "pearson") { CorPearsonControl corControl; - r = ccaPP(x, y, k, corControl, ppControl, false, false, A, B); + r = sMaxCorPP(x, y, corControl, ppControl, false, false, + A, B, objective); } else { error("method not available"); } + lambdaX = ppControl.lambdaX; lambdaY = ppControl.lambdaY; } else { error("algorithm not available"); } // wrap and return result return List::create( - Named("cor") = r, + Named("cor") = wrap(r.begin(), r.end()), Named("A") = A, - Named("B") = B + Named("B") = B, + Named("lambdaX") = wrap(lambdaX.begin(), lambdaX.end()), + Named("lambdaY") = wrap(lambdaY.begin(), lambdaY.end()), + Named("objective") = wrap(objective.begin(), objective.end()) ); } diff --git a/src/cca.h b/src/cca.h index 781b0ec..9ce6292 100644 --- a/src/cca.h +++ b/src/cca.h @@ -16,6 +16,8 @@ using namespace Rcpp; // functions to export to R RcppExport SEXP R_fastCor(SEXP R_x, SEXP R_y, SEXP R_method, SEXP R_control); // for testing RcppExport SEXP R_ccaPP(SEXP R_x, SEXP R_y, SEXP R_k, SEXP R_method, SEXP R_corControl, - SEXP R_algorithm, SEXP R_ppControl, SEXP R_fallback); + SEXP R_algorithm, SEXP R_ppControl, SEXP R_fallback); +RcppExport SEXP R_sMaxCorPP(SEXP R_x, SEXP R_y, SEXP R_method, SEXP R_corControl, + SEXP R_algorithm, SEXP R_ppControl, SEXP R_fallback); #endif From d58e04ef3011d2c80eb7087e8946ccac641d23f6 Mon Sep 17 00:00:00 2001 From: Andreas Alfons Date: Mon, 26 Nov 2012 12:00:56 +0100 Subject: [PATCH 05/14] working code for sparse maximum correlation with multivariate x and y --- src/cca.cpp | 77 +++++++++++++++++++++++++---------------------------- 1 file changed, 36 insertions(+), 41 deletions(-) diff --git a/src/cca.cpp b/src/cca.cpp index 43d38d6..289b160 100644 --- a/src/cca.cpp +++ b/src/cca.cpp @@ -898,47 +898,42 @@ vec SparseGridControl::maxCor(const mat& x, const mat& y, } } else if((p > 1) && (q > 1)) { // both data sets are multivariate -// uword nLambdaX = lambdaX.n_elem, nLambdaY = lambdaY.n_elem; -// uword nLambda = nLambdaX * nLambdaY; -// r.set_size(nLambda); objective.set_size(nLambda); -// A.ones(p, nLambda); B.set_size(q, nLambda); -// // find order of x and y variables -// uvec orderX(p), orderY(q); -// vec a = zeros(p), b = zeros(q); -// bool startWithX; -// findOrder(x, y, corControl, orderX, orderY, maxCor, a, b, startWithX); -// // compute initial value of objective function (L1 norms are 1) -// // compute maximum correlation for each combination of penalty parameters -// for(int k = nLambdaX-1; k >= 0; k--) { -// for(int l = nLambdaY-1; l >= 0; l--) { -// double penaltyX = lambdaX(k) * norm(a, 1); -// double penaltyY = lambdaY(l) * norm(b, 1); -// maxObjective = maxCor - penaltyX - penaltyY; -// Rprintf("lambdaX = %f, lambdaY=%f\n", lambdaX(k), lambdaY(l)); -// Rprintf(" before: maxCor = %f, maxObjective = %f\n", maxCor, maxObjective); -// // perform alternate grid searches -// if(startWithX) { -// // start with grid search for x -// maxCorFit(x, orderX, lambdaX(k), y, orderY, lambdaY(l), corControl, -// maxCor, a, b, penaltyX, penaltyY, maxObjective); -// } else { -// // start with grid search for y -// maxCorFit(y, orderY, lambdaY(l), x, orderX, lambdaX(k), corControl, -// maxCor, b, a, penaltyY, penaltyX, maxObjective); -// } -// Rprintf(" after: maxCor = %f, maxObjective = %f\n", maxCor, maxObjective); -// uword kl = k * nLambdaX + l; -// r(kl) = maxCor; -// A.col(kl) = a; -// B.col(kl) = b; -// objective(kl) = maxObjective; -// } -// // reset values -// uword kl = k * nLambdaX + nLambdaY - 1; -// maxCor = r(kl); -// a = A.col(kl); -// b = B.col(kl); -// } + uword nLambdaX = lambdaX.n_elem, nLambdaY = lambdaY.n_elem; + uword nLambda = nLambdaX * nLambdaY; + r.set_size(nLambda); objective.set_size(nLambda); + A.ones(p, nLambda); B.set_size(q, nLambda); + // find order of x and y variables + uvec orderX(p), orderY(q); + vec a = zeros(p), b = zeros(q); + bool startWithX; + findOrder(x, y, corControl, orderX, orderY, maxCor, a, b, startWithX); + // compute maximum correlation for each combination of penalty parameters + double initialMaxCor = maxCor; + vec initialA = a, initialB = b; + int kl = nLambda; + for(int k = nLambdaX-1; k >= 0; k--) { + for(int l = nLambdaY-1; l >= 0; l--) { + maxCor = initialMaxCor; + a = initialA; b = initialB; + double penaltyX = lambdaX(k), penaltyY = lambdaY(l); // L1 norms are 1 + maxObjective = maxCor - penaltyX - penaltyY; + // perform alternate grid searches + if(startWithX) { + // start with grid search for x + maxCorFit(x, orderX, lambdaX(k), y, orderY, lambdaY(l), corControl, + maxCor, a, b, penaltyX, penaltyY, maxObjective); + } else { + // start with grid search for y + maxCorFit(y, orderY, lambdaY(l), x, orderX, lambdaX(k), corControl, + maxCor, b, a, penaltyY, penaltyX, maxObjective); + } + kl--; + r(kl) = maxCor; + A.col(kl) = a; + B.col(kl) = b; + objective(kl) = maxObjective; + } + } } } // return maximum correlation From a0bf762144d3a67eb40954223153edd5bf97d665 Mon Sep 17 00:00:00 2001 From: Andreas Alfons Date: Mon, 26 Nov 2012 12:34:06 +0100 Subject: [PATCH 06/14] using the same starting values for all values of the tuning parameters --- src/cca.cpp | 38 +++++++++++++++++++------------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/src/cca.cpp b/src/cca.cpp index 289b160..3119353 100644 --- a/src/cca.cpp +++ b/src/cca.cpp @@ -855,7 +855,7 @@ vec SparseGridControl::maxCor(const mat& x, const mat& y, } objective(0) = r(0); } else { - double maxCor, maxObjective; + double initialMaxCor, maxCor, maxObjective; if((p > 1) && (q == 1)) { // x is multivariate, y is univariate uword nLambdaX = lambdaX.n_elem; @@ -865,15 +865,16 @@ vec SparseGridControl::maxCor(const mat& x, const mat& y, vec yy = y.unsafe_col(0); // reuse memory // find order of x variables uvec orderX(p); - vec a = zeros(p); - findOrder(x, yy, corControl, orderX, maxCor, a); // column order + vec initialA = zeros(p); + findOrder(x, yy, corControl, orderX, initialMaxCor, initialA); // compute maximum correlation for each penalty parameter for(int k = nLambdaX-1; k >= 0; k--) { - maxObjective = maxCor - lambdaX(k) * norm(a, 1); + maxCor = initialMaxCor; + vec a = A.unsafe_col(k); a = initialA; + maxObjective = maxCor - lambdaX(k); // L1 norm is 1 maxCorFit(x, orderX, lambdaX(k), yy, corControl, maxCor, a, maxObjective); r(k) = maxCor; - A.col(k) = a; objective(k) = maxObjective; } } else if((p == 1) && (q > 1)) { @@ -885,15 +886,16 @@ vec SparseGridControl::maxCor(const mat& x, const mat& y, vec xx = x.unsafe_col(0); // reuse memory // find order of y variables uvec orderY(q); - vec b = zeros(q); - findOrder(y, xx, corControl, orderY, maxCor, b); // column order + vec initialB = zeros(q); + findOrder(y, xx, corControl, orderY, initialMaxCor, initialB); // compute maximum correlation for each penalty parameter for(int k = nLambdaY-1; k >= 0; k--) { - maxObjective = maxCor - lambdaY(k) * norm(b, 1); + maxCor = initialMaxCor; + vec b = B.unsafe_col(k); b = initialB; + maxObjective = maxCor - lambdaY(k); // L1 norm is 1 maxCorFit(y, orderY, lambdaY(k), xx, corControl, maxCor, b, maxObjective); r(k) = maxCor; - B.col(k) = b; objective(k) = maxObjective; } } else if((p > 1) && (q > 1)) { @@ -901,19 +903,19 @@ vec SparseGridControl::maxCor(const mat& x, const mat& y, uword nLambdaX = lambdaX.n_elem, nLambdaY = lambdaY.n_elem; uword nLambda = nLambdaX * nLambdaY; r.set_size(nLambda); objective.set_size(nLambda); - A.ones(p, nLambda); B.set_size(q, nLambda); + A.set_size(p, nLambda), B.set_size(q, nLambda); // find order of x and y variables uvec orderX(p), orderY(q); - vec a = zeros(p), b = zeros(q); + vec initialA = zeros(p), initialB = zeros(q); bool startWithX; - findOrder(x, y, corControl, orderX, orderY, maxCor, a, b, startWithX); - // compute maximum correlation for each combination of penalty parameters - double initialMaxCor = maxCor; - vec initialA = a, initialB = b; - int kl = nLambda; + findOrder(x, y, corControl, orderX, orderY, initialMaxCor, + initialA, initialB, startWithX); + // compute maximum correlation for each combination of penalty parameters + int kl = nLambda-1; for(int k = nLambdaX-1; k >= 0; k--) { for(int l = nLambdaY-1; l >= 0; l--) { maxCor = initialMaxCor; + vec a = A.unsafe_col(kl), b = B.unsafe_col(kl); a = initialA; b = initialB; double penaltyX = lambdaX(k), penaltyY = lambdaY(l); // L1 norms are 1 maxObjective = maxCor - penaltyX - penaltyY; @@ -927,11 +929,9 @@ vec SparseGridControl::maxCor(const mat& x, const mat& y, maxCorFit(y, orderY, lambdaY(l), x, orderX, lambdaX(k), corControl, maxCor, b, a, penaltyY, penaltyX, maxObjective); } - kl--; r(kl) = maxCor; - A.col(kl) = a; - B.col(kl) = b; objective(kl) = maxObjective; + kl--; } } } From ef89634e1206ee408c00781da6f69e17751b6bef Mon Sep 17 00:00:00 2001 From: Andreas Alfons Date: Wed, 28 Nov 2012 14:15:02 +0100 Subject: [PATCH 07/14] finished first version of sparse maximum correlation --- DESCRIPTION | 1 + R/cca.R | 38 +------------- R/cv.R | 120 ++++++++++++++++++++++++++++++++++++++++++ R/maxCor.R | 130 +++++++++++++++++++++++++++++----------------- R/utils.R | 54 +++++++++++++++++++ man/maxCorGrid.Rd | 7 ++- src/cca.cpp | 94 ++++++++++++++------------------- 7 files changed, 302 insertions(+), 142 deletions(-) create mode 100644 R/cv.R diff --git a/DESCRIPTION b/DESCRIPTION index 0a9d11e..cd322ac 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -22,6 +22,7 @@ Authors@R: c(person("Andreas", "Alfons", email = "andreas.alfons@kuleuven.be", Collate: 'cca.R' 'cor.R' + 'cv.R' 'fastMAD.R' 'fastMedian.R' 'maxCor.R' diff --git a/R/cca.R b/R/cca.R index ad01bd3..5218942 100644 --- a/R/cca.R +++ b/R/cca.R @@ -346,43 +346,7 @@ ccaPP <- function(x, y, k = 1, # check method and get list of control arguments method <- match.arg(method) corControl <- getCorControl(method, corControl, forceConsistency) - # additional checks for grid search algorithm - if(algorithm == "grid") { - # check subset of variables to be used for determining the order of - # the variables from the respective other data set - select <- ppControl$select - ppControl$select <- NULL - if(!is.null(select)) { - if(is.list(select)) { - # make sure select is a list with two index vectors and - # drop invalid indices from each vector - select <- rep(select, length.out=2) - select <- mapply(function(indices, max) { - indices <- as.integer(indices) - indices[which(indices > 0 & indices <= max)] - 1 - }, select, c(p, q)) - valid <- sapply(select, length) > 0 - # add the two index vectors to control object - if(all(valid)) { - ppControl$selectX <- select[[1]] - ppControl$selectY <- select[[2]] - } else select <- NULL - } else { - # check number of indices to sample - select <- rep(as.integer(select), length.out=2) - valid <- !is.na(select) & select > 0 & select < c(p, q) - if(all(valid)) { - # generate index vectors and add them to control object - if(!is.null(seed)) set.seed(seed) - ppControl$selectX <- sample.int(p, select[1]) - 1 - ppControl$selectY <- sample.int(q, select[2]) - 1 - } else select <- NULL - } - } - if(is.null(select)) { - ppControl$selectX <- ppControl$selectY <- integer() - } - } + ppControl <- getPPControl(algorithm, ppControl, p=p, q=q, seed=seed) # call C++ function cca <- .Call("R_ccaPP", R_x=x, R_y=y, R_k=k, R_method=method, R_corControl=corControl, R_algorithm=algorithm, diff --git a/R/cv.R b/R/cv.R new file mode 100644 index 0000000..70ca2a2 --- /dev/null +++ b/R/cv.R @@ -0,0 +1,120 @@ +# ---------------------- +# Author: Andreas Alfons +# KU Leuven +# ---------------------- + +# perform cross-validation +cvMaxCor <- function(call, x, y, folds, corFun = corSpearman, corArgs = list(), + envir = parent.frame(), cl = NULL) { + # initializations + R <- folds$R + useParallel <- !is.null(cl) + # obtain list of predictions for all replications + if(R == 1) { + s <- getIndices(folds) + if(useParallel) xy <- pcvXY(s, call=call, x=x, y=y, envir=envir, cl=cl) + else xy <- cvXY(s, call=call, x=x, y=y, envir=envir) + corXY(xy, corFun=corFun, corArgs=corArgs) + } else { + cvFun <- function(r) { + s <- getIndices(folds, r) + xy <- cvXY(s, call=call, x=x, y=y, envir=envir) + corXY(xy, corFun=corFun, corArgs=corArgs) + } + if(useParallel) { + cv <- parLapply(cl, seq_len(R), cvFun) + cv <- do.call(cbind, cv) + } + else cv <- sapply(seq_len(R), cvFun) + apply(cv, 1, mean) + } +} + +# compute the weighting vectors for the training data and obtain linear +# combinations of the test data +rsXY <- function(i, call, x, y, envir) { + # plug training data into function call + call$x <- x[-i, , drop=FALSE] + call$y <- y[-i, , drop=FALSE] + # evaluate function call in supplied environment to make sure + # that other arguments are evaluated correctly + fit <- eval(call, envir) + # obtain linear combinations for test data + cbind(x=x[i, , drop=FALSE] %*% fit$a, y=y[i, , drop=FALSE] %*% fit$b) +} + +# compute linear combinations for each left-out block in one replication of +# cross validation +cvXY <- function(folds, call, x, y, envir) { + tmp <- lapply(folds, rsXY, call=call, x=x, y=y, envir=envir) + # instead of collecting the results from the folds in the original order + # of the observations, they are simply stacked on top of each other + do.call(rbind, tmp) +} + +# compute linear combinations for each left-out block in one replication of +# cross validation via parallel computing +pcvXY <- function(folds, call, x, y, envir, cl) { + # parLapply() already has an argument 'x', so don't use the argument names + # for the data or it will throw an error + tmp <- parLapply(cl, folds, rsXY, call, x, y, envir=envir) + # instead of collecting the results from the folds in the original order + # of the observations, they are simply stacked on top of each other + do.call(rbind, tmp) +} + +# compute the correlation between the linear combinations obtained via +# cross-validation +corXY <- function(xy, corFun = corSpearman, corArgs = list()) { + # number of linear combinations of each data set + p <- ncol(xy) / 2 + # split the data into the linear combinations of x and y, respectively, and + # compute the correlation between each pair of linear combinations + if(p == 1) doCall(corFun, xy[, 1], xy[, 2], args=corArgs) + else { + sapply(seq_len(p), function(i) { + doCall(corFun, xy[, i], xy[, p+i], args=corArgs) + }) + } +} + +# set up cross-validation folds +cvFolds <- function(n, K = 5, R = 1, + type = c("random", "consecutive", "interleaved"), + grouping = NULL) { + # check arguments + n <- if(is.null(grouping)) round(rep(n, length.out=1)) else nlevels(grouping) + if(!isTRUE(n > 0)) stop("'n' must be positive") + if(!isTRUE(K <= n)) stop(sprintf("'K' must be smaller or equal to %d", n)) + if(K == n) type <- "leave-one-out" + else type <- match.arg(type) + # obtain CV folds + if(type == "random") { + # random K-fold splits with R replications + subsets <- replicate(R, sample.int(n)) + } else { + # leave-one-out CV or non-random splits, replication not meaningful + R <- 1 + subsets <- as.matrix(seq_len(n)) + } + which <- as.factor(rep(seq_len(K), length.out=n)) + if(type == "consecutive") which <- rep.int(seq_len(K), tabulate(which)) + # construct and return object + folds <- list(n=n, K=K, R=R, subsets=subsets, which=which) + if(!is.null(grouping)) folds$grouping <- split(seq_along(grouping), grouping) + class(folds) <- "cvFolds" + folds +} + +# retrieve list of indices for r-th replication +getIndices <- function(x, r = 1) { + # split permuted items according to the folds + subsets <- split(x$subsets[, r], x$which) + # in case of grouped data, the list contains the group indices in each CV + # fold, so the indices of the respective observations need to be extracted + if(!is.null(grouping <- x$grouping)) + subsets <- lapply(subsets, function(s) unlist(grouping[s], use.names=FALSE)) + # return list of indices for CV folds + names(subsets) <- NULL + subsets +} diff --git a/R/maxCor.R b/R/maxCor.R index cd99fc7..f92b4ea 100644 --- a/R/maxCor.R +++ b/R/maxCor.R @@ -140,18 +140,19 @@ maxCorGrid <- function(x, y, sMaxCorGrid <- function(x, y, lambdaX = 0, lambdaY = 0, method = c("spearman", "kendall", "quadrant", "M", "pearson"), - control = list(...), nIterations = 10, nAlternate = 10, - nGrid = 25, select = NULL, tol = 1e-06, + control = list(...), nIterations = 10, + nAlternate = 10, nGrid = 25, select = NULL, + tol = 1e-06, folds = NULL, K = 5, R = 1, + type = c("random", "consecutive", "interleaved"), + grouping = NULL, nCores = 1, cl = NULL, fallback = FALSE, seed = NULL, ...) { ## initializations matchedCall <- match.call() ## define list of control arguments for algorithm lambdaX <- as.numeric(lambdaX) lambdaX[is.na(lambdaX)] <- formals()$lambdaX - lambdaX <- sort(unique(lambdaX)) lambdaY <- as.numeric(lambdaY) lambdaY[is.na(lambdaY)] <- formals()$lambdaY - lambdaY <- sort(unique(lambdaY)) nIterations <- as.integer(nIterations) nAlternate <- as.integer(nAlternate) nGrid <- as.integer(nGrid) @@ -160,7 +161,9 @@ sMaxCorGrid <- function(x, y, lambdaX = 0, lambdaY = 0, nAlternate=nAlternate, nGrid=nGrid, select=select, tol=tol) ## call workhorse function maxCor <- sMaxCorPP(x, y, method=method, corControl=control, algorithm="grid", - ppControl=ppControl, fallback=fallback, seed=seed) + ppControl=ppControl, folds=folds, K=K, R=R, type=type, + grouping=grouping, nCores=nCores, cl=cl, + fallback=fallback, seed=seed) maxCor$call <- matchedCall maxCor } @@ -268,10 +271,12 @@ maxCorPP <- function(x, y, ...) { ## workhorse function for sparse maximum correlation +#' @import parallel sMaxCorPP <- function(x, y, method = c("spearman", "kendall", "quadrant", "M", "pearson"), - corControl, algorithm = "grid", ppControl, - fallback = FALSE, seed = NULL) { + corControl, algorithm = "grid", ppControl, folds = NULL, + nCores = 1, cl = NULL, fallback = FALSE, seed = NULL, + ...) { ## initializations x <- as.matrix(x) y <- as.matrix(y) @@ -282,55 +287,82 @@ sMaxCorPP <- function(x, y, method = c("spearman", "kendall", "quadrant", ## prepare the data and call C++ function if(n == 0 || p == 0 || q == 0) { # zero dimension - maxCor <- list(cor=NA, a=numeric(), b=numeric(), objective=NA) + maxCor <- list(cor=NA, a=numeric(), b=numeric(), lambdaX=NA, + lambdaY=NA, objective=NA) } else { - # check method and get list of control arguments + # get list of control arguments method <- match.arg(method) - corControl <- getCorControl(method, corControl, forceConsistency) - # additional checks for grid search algorithm - if(algorithm == "grid") { - # check subset of variables to be used for determining the order of - # the variables from the respective other data set - select <- ppControl$select - ppControl$select <- NULL - if(!is.null(select)) { - if(is.list(select)) { - # make sure select is a list with two index vectors and - # drop invalid indices from each vector - select <- rep(select, length.out=2) - select <- mapply(function(indices, max) { - indices <- as.integer(indices) - indices[which(indices > 0 & indices <= max)] - 1 - }, select, c(p, q)) - valid <- sapply(select, length) > 0 - # add the two index vectors to control object - if(all(valid)) { - ppControl$selectX <- select[[1]] - ppControl$selectY <- select[[2]] - } else select <- NULL - } else { - # check number of indices to sample - select <- rep(as.integer(select), length.out=2) - valid <- !is.na(select) & select > 0 & select < c(p, q) - if(all(valid)) { - # generate index vectors and add them to control object - if(!is.null(seed)) set.seed(seed) - ppControl$selectX <- sample.int(p, select[1]) - 1 - ppControl$selectY <- sample.int(q, select[2]) - 1 - } else select <- NULL - } + corControl <- getCorControl(method, corControl, forceConsistency=FALSE) + fallback <- isTRUE(fallback) + if(!is.null(seed)) set.seed(seed) + ppControl <- getPPControl(algorithm, ppControl, p=p, q=q) + # check grid of values for the penalty parameters + lambdaX <- ppControl$lambdaX + lambdaY <- ppControl$lambdaY + ppControl$lambdaX <- NULL + ppControl$lambdaY <- NULL + lambdaX <- if(p == 1) 0 else sort(unique(lambdaX)) + lambdaY <- if(q == 1) 0 else sort(unique(lambdaY)) + lambda <- as.matrix(expand.grid(lambdaX=lambdaX, lambdaY=lambdaY)) + ppControl$lambda <- lambda + # check whether cross-validation over grid of tuning parameters is necessary + useCV <- nrow(lambda) > 1 + if(useCV) { + if(is.null(folds)) folds <- cvFolds(n, ...) + # set up parallel computing if requested + haveCl <- inherits(cl, "cluster") + if(haveCl) haveNCores <- FALSE + else { + if(is.na(nCores)) nCores <- detectCores() # use all available cores + if(!is.numeric(nCores) || is.infinite(nCores) || nCores < 1) { + nCores <- 1 # use default value + warning("invalid value of 'nCores'; using default value") + } else nCores <- as.integer(nCores) + R <- folds$R + nCores <- if(R == 1) min(nCores, folds$K) else min(nCores, R) + haveNCores <- nCores > 1 } - if(is.null(select)) { - ppControl$selectX <- ppControl$selectY <- integer() + # check whether parallel computing should be used + useParallel <- haveNCores || haveCl + # set up multicore or snow cluster if not supplied + if(haveNCores) { + if(.Platform$OS.type == "windows") { + cl <- makePSOCKcluster(rep.int("localhost", nCores)) + } else cl <- makeForkCluster(nCores) + on.exit(stopCluster(cl)) } + # perform cross-validation over grid of tuning parameters + call <- call("sMaxCorPPFit", method=method, corControl=corControl, + algorithm=algorithm, ppControl=ppControl, fallback=fallback) + call[[1]] <- sMaxCorPPFit # necessary to work with parallel computing + corFun <- switch(method, spearman=corSpearman, kendall=corKendall, + quadrant=corQuadrant, M=corM, pearson=corPearson) + cv <- cvMaxCor(call, x, y, folds=folds, corFun=corFun, + corArgs=corControl, cl=cl) + # determine the optimal tuning parameters + ppControl$lambda <- lambda[which.max(cv), , drop=FALSE] } - # call C++ function - maxCor <- .Call("R_sMaxCorPP", R_x=x, R_y=y, R_method=method, - R_corControl=corControl, R_algorithm=algorithm, - R_ppControl=ppControl, R_fallback=isTRUE(fallback), - PACKAGE="ccaPP") + # compute the final solution + maxCor <- sMaxCorPPFit(x, y, method=method, corControl=corControl, + algorithm=algorithm, ppControl=ppControl, + fallback=fallback) + maxCor$a <- drop(maxCor$a) + maxCor$b <- drop(maxCor$b) + if(useCV) maxCor$cv <- cbind(lambda, CV=cv) } ## assign class and return results class(maxCor) <- "maxCor" maxCor } + +## simple wrapper for calling the C++ function +sMaxCorPPFit <- function(x, y, method = c("spearman", "kendall", "quadrant", + "M", "pearson"), + corControl, algorithm = "grid", ppControl, + fallback = FALSE) { + # call C++ function + maxCor <- .Call("R_sMaxCorPP", R_x=x, R_y=y, R_method=method, + R_corControl=corControl, R_algorithm=algorithm, + R_ppControl=ppControl, R_fallback=fallback, + PACKAGE="ccaPP") +} diff --git a/R/utils.R b/R/utils.R index b300e20..d603d23 100644 --- a/R/utils.R +++ b/R/utils.R @@ -18,6 +18,17 @@ checkIndices <- function(indices, max) { indices[which(indices > 0 & indices <= max)] } +## call a function by either +# 1) simply evaluating a supplied function for the basic arguments if there are +# no additional arguments in list format +# 2) evaluating a supplied function with 'do.call' if there are additional +# arguments in list format +doCall <- function(fun, ..., args = list()) { + if(length(args) == 0) { + fun(...) + } else do.call(fun, c(list(...), args)) +} + ## call C++ to compute ranks of observations in a vector (for testing) fastRank <- function(x) { x <- as.numeric(x) @@ -61,6 +72,49 @@ getCorControl <- function(method, control, forceConsistency = TRUE) { out } +## get list of control arguments for projection pursuit algorithm +getPPControl <- function(algorithm, control, p, q, seed = NULL) { + # additional checks for grid search algorithm + if(algorithm == "grid") { + # check subset of variables to be used for determining the order of + # the variables from the respective other data set + select <- control$select + control$select <- NULL + if(!is.null(select)) { + if(is.list(select)) { + # make sure select is a list with two index vectors and + # drop invalid indices from each vector + select <- rep(select, length.out=2) + select <- mapply(function(indices, max) { + indices <- as.integer(indices) + indices[which(indices > 0 & indices <= max)] - 1 + }, select, c(p, q)) + valid <- sapply(select, length) > 0 + # add the two index vectors to control object + if(all(valid)) { + control$selectX <- select[[1]] + control$selectY <- select[[2]] + } else select <- NULL + } else { + # check number of indices to sample + select <- rep(as.integer(select), length.out=2) + valid <- !is.na(select) & select > 0 & select < c(p, q) + if(all(valid)) { + # generate index vectors and add them to control object + if(!is.null(seed)) set.seed(seed) + control$selectX <- sample.int(p, select[1]) - 1 + control$selectY <- sample.int(q, select[2]) - 1 + } else select <- NULL + } + } + if(is.null(select)) { + control$selectX <- control$selectY <- integer() + } + } + # return list of control arguments + control +} + ## L1 median (for testing) l1Median <- function(x) { # initializations diff --git a/man/maxCorGrid.Rd b/man/maxCorGrid.Rd index 873fab8..22ab85d 100644 --- a/man/maxCorGrid.Rd +++ b/man/maxCorGrid.Rd @@ -10,10 +10,13 @@ nGrid = 25, select = NULL, tol = 1e-06, fallback = FALSE, seed = NULL, ...) - sMaxCorGrid(x, y, lambda = 0, + sMaxCorGrid(x, y, lambdaX = 0, lambdaY = 0, method = c("spearman", "kendall", "quadrant", "M", "pearson"), control = list(...), nIterations = 10, nAlternate = 10, - nGrid = 25, select = NULL, tol = 1e-06, + nGrid = 25, select = NULL, tol = 1e-06, folds = NULL, + K = 5, R = 1, + type = c("random", "consecutive", "interleaved"), + grouping = NULL, nCores = 1, cl = NULL, fallback = FALSE, seed = NULL, ...) } \arguments{ diff --git a/src/cca.cpp b/src/cca.cpp index 3119353..b1c81ee 100644 --- a/src/cca.cpp +++ b/src/cca.cpp @@ -618,7 +618,7 @@ double GridControl::maxCor(const mat& x, const mat& y, CorControl& corControl, class SparseGridControl: public GridControl { public: - vec lambdaX, lambdaY; // penalty parameters + vec lambda; // grid of values for penalty parameters // constructors SparseGridControl(); SparseGridControl(List&); @@ -649,14 +649,12 @@ class SparseGridControl: public GridControl { // constructors inline SparseGridControl::SparseGridControl() { - lambdaX = zeros(1); - lambdaY = zeros(1); + lambda = zeros(1, 2); } inline SparseGridControl::SparseGridControl(List& control) : GridControl(control) { - SEXP R_lambdaX = control["lambdaX"], R_lambdaY = control["lambdaY"]; - lambdaX = as(R_lambdaX); - lambdaY = as(R_lambdaY); + SEXP R_lambda = control["lambda"]; + lambda = as(R_lambda); } // compute updated weighting vector in grid search @@ -790,8 +788,8 @@ void SparseGridControl::maxCorFit(const mat& x, const uvec& orderX, // workhorse function for maximum correlation between multivariate x and y template void SparseGridControl::maxCorFit(const mat& x, const uvec& orderX, - const double& lambda1, const mat& y, const uvec& orderY, - const double& lambda2, CorControl& corControl, double& maxCor, vec& a, + const double& lambdaX, const mat& y, const uvec& orderY, + const double& lambdaY, CorControl& corControl, double& maxCor, vec& a, vec& b, double& penaltyX, double& penaltyY, double& maxObjective) { // initializations double previousMaxObjective = R_NegInf; @@ -806,11 +804,11 @@ void SparseGridControl::maxCorFit(const mat& x, const uvec& orderX, altMaxObjective = maxObjective; // maximize correlation functional over a keeping b fixed vec yb = y * b; // linear combination of columns of y - gridSearch(x, orderX, lambda1, yb, penaltyY, corControl, + gridSearch(x, orderX, lambdaX, yb, penaltyY, corControl, grid, maxCor, a, penaltyX, maxObjective); // maximize correlation functional over b keeping a fixed vec xa = x * a; // linear combination of columns of x - gridSearch(y, orderY, lambda2, xa, penaltyX, corControl, + gridSearch(y, orderY, lambdaY, xa, penaltyX, corControl, grid, maxCor, b, penaltyY, maxObjective); j++; } @@ -844,7 +842,6 @@ vec SparseGridControl::maxCor(const mat& x, const mat& y, // both data sets are univariate r.set_size(1); objective.set_size(1); A.ones(p, 1); B.ones(q, 1); - lambdaX.zeros(p); lambdaY.zeros(q); vec xx = x.unsafe_col(0), yy = y.unsafe_col(0); // reuse memory // compute correlation r(0) = corControl.cor(xx, yy); @@ -856,53 +853,46 @@ vec SparseGridControl::maxCor(const mat& x, const mat& y, objective(0) = r(0); } else { double initialMaxCor, maxCor, maxObjective; - if((p > 1) && (q == 1)) { + uword nLambda = lambda.n_rows; + r.set_size(nLambda); objective.set_size(nLambda); + if((p > 1) && (q == 1)) { // x is multivariate, y is univariate - uword nLambdaX = lambdaX.n_elem; - r.set_size(nLambdaX); objective.set_size(nLambdaX); - A.set_size(p, nLambdaX); B.ones(q, nLambdaX); - lambdaY.zeros(q); + A.set_size(p, nLambda); B.ones(q, nLambda); vec yy = y.unsafe_col(0); // reuse memory // find order of x variables uvec orderX(p); vec initialA = zeros(p); findOrder(x, yy, corControl, orderX, initialMaxCor, initialA); // compute maximum correlation for each penalty parameter - for(int k = nLambdaX-1; k >= 0; k--) { + for(uword k = 0; k < nLambda; k++) { maxCor = initialMaxCor; vec a = A.unsafe_col(k); a = initialA; - maxObjective = maxCor - lambdaX(k); // L1 norm is 1 - maxCorFit(x, orderX, lambdaX(k), yy, corControl, + maxObjective = maxCor - lambda(k,0); // L1 norm is 1 + maxCorFit(x, orderX, lambda(k,0), yy, corControl, maxCor, a, maxObjective); r(k) = maxCor; objective(k) = maxObjective; } } else if((p == 1) && (q > 1)) { // x is univariate, y is multivariate - uword nLambdaY = lambdaY.n_elem; - r.set_size(nLambdaY); objective.set_size(nLambdaY); - A.ones(p, nLambdaY); B.set_size(q, nLambdaY); - lambdaX.zeros(p); + A.ones(p, nLambda); B.set_size(q, nLambda); vec xx = x.unsafe_col(0); // reuse memory // find order of y variables uvec orderY(q); vec initialB = zeros(q); findOrder(y, xx, corControl, orderY, initialMaxCor, initialB); // compute maximum correlation for each penalty parameter - for(int k = nLambdaY-1; k >= 0; k--) { + for(uword k = 0; k < nLambda; k++) { maxCor = initialMaxCor; vec b = B.unsafe_col(k); b = initialB; - maxObjective = maxCor - lambdaY(k); // L1 norm is 1 - maxCorFit(y, orderY, lambdaY(k), xx, corControl, + maxObjective = maxCor - lambda(k,1); // L1 norm is 1 + maxCorFit(y, orderY, lambda(k,1), xx, corControl, maxCor, b, maxObjective); r(k) = maxCor; objective(k) = maxObjective; } } else if((p > 1) && (q > 1)) { // both data sets are multivariate - uword nLambdaX = lambdaX.n_elem, nLambdaY = lambdaY.n_elem; - uword nLambda = nLambdaX * nLambdaY; - r.set_size(nLambda); objective.set_size(nLambda); A.set_size(p, nLambda), B.set_size(q, nLambda); // find order of x and y variables uvec orderX(p), orderY(q); @@ -911,28 +901,24 @@ vec SparseGridControl::maxCor(const mat& x, const mat& y, findOrder(x, y, corControl, orderX, orderY, initialMaxCor, initialA, initialB, startWithX); // compute maximum correlation for each combination of penalty parameters - int kl = nLambda-1; - for(int k = nLambdaX-1; k >= 0; k--) { - for(int l = nLambdaY-1; l >= 0; l--) { - maxCor = initialMaxCor; - vec a = A.unsafe_col(kl), b = B.unsafe_col(kl); - a = initialA; b = initialB; - double penaltyX = lambdaX(k), penaltyY = lambdaY(l); // L1 norms are 1 - maxObjective = maxCor - penaltyX - penaltyY; - // perform alternate grid searches - if(startWithX) { - // start with grid search for x - maxCorFit(x, orderX, lambdaX(k), y, orderY, lambdaY(l), corControl, - maxCor, a, b, penaltyX, penaltyY, maxObjective); - } else { - // start with grid search for y - maxCorFit(y, orderY, lambdaY(l), x, orderX, lambdaX(k), corControl, - maxCor, b, a, penaltyY, penaltyX, maxObjective); - } - r(kl) = maxCor; - objective(kl) = maxObjective; - kl--; - } + for(uword k = 0; k < nLambda; k++) { + maxCor = initialMaxCor; + vec a = A.unsafe_col(k), b = B.unsafe_col(k); + a = initialA; b = initialB; + double penaltyX = lambda(k,0), penaltyY = lambda(k,1); // L1 norms are 1 + maxObjective = maxCor - penaltyX - penaltyY; + // perform alternate grid searches + if(startWithX) { + // start with grid search for x + maxCorFit(x, orderX, lambda(k,0), y, orderY, lambda(k,1), corControl, + maxCor, a, b, penaltyX, penaltyY, maxObjective); + } else { + // start with grid search for y + maxCorFit(y, orderY, lambda(k,1), x, orderX, lambda(k,0), corControl, + maxCor, b, a, penaltyY, penaltyX, maxObjective); + } + r(k) = maxCor; + objective(k) = maxObjective; } } } @@ -1376,15 +1362,15 @@ SEXP R_sMaxCorPP(SEXP R_x, SEXP R_y, SEXP R_method, SEXP R_corControl, } else { error("method not available"); } - lambdaX = ppControl.lambdaX; lambdaY = ppControl.lambdaY; + lambdaX = ppControl.lambda.col(0); lambdaY = ppControl.lambda.col(1); } else { error("algorithm not available"); } // wrap and return result return List::create( Named("cor") = wrap(r.begin(), r.end()), - Named("A") = A, - Named("B") = B, + Named("a") = A, + Named("b") = B, Named("lambdaX") = wrap(lambdaX.begin(), lambdaX.end()), Named("lambdaY") = wrap(lambdaY.begin(), lambdaY.end()), Named("objective") = wrap(objective.begin(), objective.end()) From 8ade4e105b804bdca8bc799fd3ea303a6ac3998e Mon Sep 17 00:00:00 2001 From: Andreas Alfons Date: Wed, 28 Nov 2012 16:03:13 +0100 Subject: [PATCH 08/14] added documentation for sparse maximum correlation --- NAMESPACE | 2 ++ R/cv.R | 54 +++++++++++++++++++++++++++++++-- R/maxCor.R | 56 ++++++++++++++++++++++++++-------- R/print.R | 23 ++++++++++++++ man/cvFolds.Rd | 77 +++++++++++++++++++++++++++++++++++++++++++++++ man/maxCorGrid.Rd | 65 ++++++++++++++++++++++++++++++++++----- 6 files changed, 256 insertions(+), 21 deletions(-) create mode 100644 man/cvFolds.Rd diff --git a/NAMESPACE b/NAMESPACE index 38e6bc5..a3361a8 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,4 +1,5 @@ S3method(print,cca) +S3method(print,cvFolds) S3method(print,maxCor) S3method(print,permTest) export(CCAgrid) @@ -10,6 +11,7 @@ export(corM) export(corPearson) export(corQuadrant) export(corSpearman) +export(cvFolds) export(fastMAD) export(fastMedian) export(maxCorGrid) diff --git a/R/cv.R b/R/cv.R index 70ca2a2..39e2d2d 100644 --- a/R/cv.R +++ b/R/cv.R @@ -3,7 +3,7 @@ # KU Leuven # ---------------------- -# perform cross-validation +## perform cross-validation cvMaxCor <- function(call, x, y, folds, corFun = corSpearman, corArgs = list(), envir = parent.frame(), cl = NULL) { # initializations @@ -78,7 +78,57 @@ corXY <- function(xy, corFun = corSpearman, corArgs = list()) { } } -# set up cross-validation folds + +#' Cross-validation folds +#' +#' Split observations or groups of observations into \eqn{K} folds to be used +#' for (repeated) \eqn{K}-fold cross-validation. \eqn{K} should thereby be +#' chosen such that all folds are of approximately equal size. +#' +#' @aliases print.cvFolds +#' +#' @param n an integer giving the number of observations to be split into +#' folds. This is ignored if \code{grouping} is supplied in order to split +#' groups of observations into folds. +#' @param K an integer giving the number of folds into which the observations +#' should be split (the default is five). Setting \code{K} equal to the number +#' of observations or groups yields leave-one-out cross-validation. +#' @param R an integer giving the number of replications for repeated +#' \eqn{K}-fold cross-validation. This is ignored for for leave-one-out +#' cross-validation and other non-random splits of the data. +#' @param type a character string specifying the type of folds to be +#' generated. Possible values are \code{"random"} (the default), +#' \code{"consecutive"} or \code{"interleaved"}. +#' @param grouping a factor specifying groups of observations. If supplied, +#' the data are split according to the groups rather than individual +#' observations such that all observations within a group belong to the same +#' fold. +#' +#' @returnClass cvFolds +#' @returnItem n an integer giving the number of observations or groups. +#' @returnItem K an integer giving the number of folds. +#' @returnItem R an integer giving the number of replications. +#' @returnItem subsets an integer matrix in which each column contains a +#' permutation of the indices of the observations or groups. +#' @returnItem which an integer vector giving the fold for each permuted +#' observation or group. +#' @returnItem grouping a list giving the indices of the observations +#' belonging to each group. This is only returned if a grouping factor +#' has been supplied. +#' +#' @author Andreas Alfons +#' +#' @seealso \code{\link{sMaxCorGrid}} +#' +#' @examples +#' set.seed(1234) # set seed for reproducibility +#' cvFolds(20, K = 5) +#' cvFolds(20, K = 5, R = 10) +#' +#' @keywords utilities +#' +#' @export + cvFolds <- function(n, K = 5, R = 1, type = c("random", "consecutive", "interleaved"), grouping = NULL) { diff --git a/R/maxCor.R b/R/maxCor.R index f92b4ea..05f87bc 100644 --- a/R/maxCor.R +++ b/R/maxCor.R @@ -9,7 +9,7 @@ #' based on alternating series of grid searches in two-dimensional subspaces of #' each data set, with a focus on robust and nonparametric methods. #' -#' The algorithm is based on alternating series of grid searches in +#' \code{maxCorGrid} is based on alternating series of grid searches in #' two-dimensional subspaces of each data set. In each grid search, #' \code{nGrid} grid points on the unit circle in the corresponding plane are #' obtained, and the directions from the center to each of the grid points are @@ -33,9 +33,19 @@ #' selected subsets of variables, or to specify the subsets of variables #' directly. #' +#' \code{sMaxCorGrid} enforces sparsity by adding \eqn{L_{1}}{L1} penalties on +#' the weighting vectors to the objective function (i.e., the correlation of +#' the resulting projections). The optimal combination of penalty paramters is +#' thereby determined via cross-validation over a grid of values. To increase +#' the computational performance of this cross-validation procedure, parallel +#' computing is implemented via package \pkg{parallel}. +#' #' @aliases print.maxCor #' #' @param x,y each can be a numeric vector, matrix or data frame. +#' @param lambdaX,lambdaY numeric vectors of non-negative values giving the +#' penalty parameters for the weighting vectors of \code{x} and \code{y}, +#' respectively. #' @param method a character string specifying the correlation functional to #' maximize. Possible values are \code{"spearman"} for the Spearman #' correlation, \code{"kendall"} for the Kendall correlation, \code{"quadrant"} @@ -61,6 +71,19 @@ #' the second integer (see \dQuote{Details}). #' @param tol a small positive numeric value to be used for determining #' convergence. +#' @param K,R,type,grouping additional arguments for generating +#' cross-validation folds (see \code{\link{cvFolds}}). +#' @param folds an object of class \code{"cvFolds"} (as returned by +#' \code{\link{cvFolds}}) defining the folds of the data for cross-validation. +#' If supplied, this is preferred over the arguments for generating +#' cross-validation folds. +#' @param nCores a positive integer giving the number of processor cores to be +#' used for parallel computing in cross-validation (the default is 1 for no +#' parallelization). If this is set to \code{NA}, all available processor +#' cores are used. +#' @param cl a \pkg{parallel} cluster for parallel computing to be used in +#' cross-validation, as generated by \code{\link[parallel]{makeCluster}}. If +#' supplied, this is preferred over \code{nCores}. #' @param fallback logical; if a correlation functional other than the #' Pearson correlation is maximized, the data are first robustly standardized #' via median and MAD. This indicates whether standardization via mean and @@ -79,6 +102,16 @@ #' @returnItem cor a numeric giving the maximum correlation estimate. #' @returnItem a numeric; the weighting vector for \code{x}. #' @returnItem b numeric; the weighting vector for \code{y}. +#' @returnItem lambdaX a numeric giving the optimal penalty parameter for the +#' weighting vector of \code{x} (only \code{sMaxCorGrid}). +#' @returnItem lambdaY a numeric giving the optimal penalty parameter for the +#' weighting vector of \code{y} (only \code{sMaxCorGrid}). +#' @returnItem objective a numeric giving the value of the objective function +#' for the optimal combination of penalty parameters (only \code{sMaxCorGrid}). +#' @returnItem cv a numeric matrix containing the (average) results from +#' cross-validation for each combination of penalty parameters (only +#' \code{sMaxCorGrid} and if a grid of values for the penalty parameters has +#' been supplied). #' @returnItem call the matched function call. #' #' @author Andreas Alfons @@ -99,11 +132,11 @@ #' y <- xy[, (p+1):m] #' #' ## Spearman correlation -#' maxCorGrid(x, y, method = "spearman") -#' maxCorGrid(x, y, method = "spearman", consistent = TRUE) +#' maxCorProj(x, y, method = "spearman") +#' maxCorProj(x, y, method = "spearman", consistent = TRUE) #' #' ## Pearson correlation -#' maxCorGrid(x, y, method = "pearson") +#' maxCorProj(x, y, method = "pearson") #' #' @keywords multivariate robust #' @@ -140,11 +173,10 @@ maxCorGrid <- function(x, y, sMaxCorGrid <- function(x, y, lambdaX = 0, lambdaY = 0, method = c("spearman", "kendall", "quadrant", "M", "pearson"), - control = list(...), nIterations = 10, - nAlternate = 10, nGrid = 25, select = NULL, - tol = 1e-06, folds = NULL, K = 5, R = 1, + control = list(...), nIterations = 10, nAlternate = 10, + nGrid = 25, select = NULL, tol = 1e-06, K = 5, R = 1, type = c("random", "consecutive", "interleaved"), - grouping = NULL, nCores = 1, cl = NULL, + grouping = NULL, folds = NULL, nCores = 1, cl = NULL, fallback = FALSE, seed = NULL, ...) { ## initializations matchedCall <- match.call() @@ -160,10 +192,10 @@ sMaxCorGrid <- function(x, y, lambdaX = 0, lambdaY = 0, ppControl <- list(lambdaX=lambdaX, lambdaY=lambdaY, nIterations=nIterations, nAlternate=nAlternate, nGrid=nGrid, select=select, tol=tol) ## call workhorse function - maxCor <- sMaxCorPP(x, y, method=method, corControl=control, algorithm="grid", - ppControl=ppControl, folds=folds, K=K, R=R, type=type, - grouping=grouping, nCores=nCores, cl=cl, - fallback=fallback, seed=seed) + maxCor <- sMaxCorPP(x, y, method=method, corControl=control, + algorithm="grid", ppControl=ppControl, K=K, + R=R, type=type, grouping=grouping, folds=folds, + nCores=nCores, cl=cl, fallback=fallback, seed=seed) maxCor$call <- matchedCall maxCor } diff --git a/R/print.R b/R/print.R index 52aa3a3..df61bd5 100644 --- a/R/print.R +++ b/R/print.R @@ -17,6 +17,29 @@ print.cca <- function(x, ...) { invisible(x) } +#' @S3method print cvFolds +print.cvFolds <- function(x, ...) { + # print general information + cvText <- if(x$n == x$K) "Leave-one-out CV" else sprintf("%d-fold CV", x$K) + if(x$R > 1) cvText <- paste(cvText, "with", x$R, "replications") + cat(paste("\n", cvText, ":", sep="")) + # print information on folds (add space between folds and subsets) + subsets <- x$subsets + if(x$R == 1) { + cn <- if("grouping" %in% names(x)) "Group index" else "Index" + nblanks <- 2 + } else { + cn <- as.character(seq_len(x$R)) + nblanks <- 3 + } + nblanks <- max(nchar(as.character(subsets[, 1]))-nchar(cn[1]), 0) + nblanks + cn[1] <- paste(c(rep.int(" ", nblanks), cn[1]), collapse="") + dimnames(subsets) <- list(Fold=x$which, cn) + print(subsets, ...) + # return object invisibly + invisible(x) +} + #' @S3method print maxCor print.maxCor <- function(x, ...) { # print function call diff --git a/man/cvFolds.Rd b/man/cvFolds.Rd new file mode 100644 index 0000000..842bfe0 --- /dev/null +++ b/man/cvFolds.Rd @@ -0,0 +1,77 @@ +\name{cvFolds} +\alias{cvFolds} +\alias{print.cvFolds} +\title{Cross-validation folds} +\usage{ + cvFolds(n, K = 5, R = 1, + type = c("random", "consecutive", "interleaved"), + grouping = NULL) +} +\arguments{ + \item{n}{an integer giving the number of observations to + be split into folds. This is ignored if \code{grouping} + is supplied in order to split groups of observations into + folds.} + + \item{K}{an integer giving the number of folds into which + the observations should be split (the default is five). + Setting \code{K} equal to the number of observations or + groups yields leave-one-out cross-validation.} + + \item{R}{an integer giving the number of replications for + repeated \eqn{K}-fold cross-validation. This is ignored + for for leave-one-out cross-validation and other + non-random splits of the data.} + + \item{type}{a character string specifying the type of + folds to be generated. Possible values are + \code{"random"} (the default), \code{"consecutive"} or + \code{"interleaved"}.} + + \item{grouping}{a factor specifying groups of + observations. If supplied, the data are split according + to the groups rather than individual observations such + that all observations within a group belong to the same + fold.} +} +\value{ + An object of class \code{"cvFolds"} with the following + components: + + \item{n}{an integer giving the number of observations or + groups.} + + \item{K}{an integer giving the number of folds.} + + \item{R}{an integer giving the number of replications.} + + \item{subsets}{an integer matrix in which each column + contains a permutation of the indices of the observations + or groups.} + + \item{which}{an integer vector giving the fold for each + permuted observation or group.} + + \item{grouping}{a list giving the indices of the + observations belonging to each group. This is only + returned if a grouping factor has been supplied.} +} +\description{ + Split observations or groups of observations into \eqn{K} + folds to be used for (repeated) \eqn{K}-fold + cross-validation. \eqn{K} should thereby be chosen such + that all folds are of approximately equal size. +} +\examples{ +set.seed(1234) # set seed for reproducibility +cvFolds(20, K = 5) +cvFolds(20, K = 5, R = 10) +} +\author{ + Andreas Alfons +} +\seealso{ + \code{\link{sMaxCorGrid}} +} +\keyword{utilities} + diff --git a/man/maxCorGrid.Rd b/man/maxCorGrid.Rd index 22ab85d..815e90d 100644 --- a/man/maxCorGrid.Rd +++ b/man/maxCorGrid.Rd @@ -13,16 +13,19 @@ sMaxCorGrid(x, y, lambdaX = 0, lambdaY = 0, method = c("spearman", "kendall", "quadrant", "M", "pearson"), control = list(...), nIterations = 10, nAlternate = 10, - nGrid = 25, select = NULL, tol = 1e-06, folds = NULL, - K = 5, R = 1, + nGrid = 25, select = NULL, tol = 1e-06, K = 5, R = 1, type = c("random", "consecutive", "interleaved"), - grouping = NULL, nCores = 1, cl = NULL, + grouping = NULL, folds = NULL, nCores = 1, cl = NULL, fallback = FALSE, seed = NULL, ...) } \arguments{ \item{x,y}{each can be a numeric vector, matrix or data frame.} + \item{lambdaX,lambdaY}{numeric vectors of non-negative + values giving the penalty parameters for the weighting + vectors of \code{x} and \code{y}, respectively.} + \item{method}{a character string specifying the correlation functional to maximize. Possible values are \code{"spearman"} for the Spearman correlation, @@ -63,6 +66,27 @@ \item{tol}{a small positive numeric value to be used for determining convergence.} + \item{K,R,type,grouping}{additional arguments for + generating cross-validation folds (see + \code{\link{cvFolds}}).} + + \item{folds}{an object of class \code{"cvFolds"} (as + returned by \code{\link{cvFolds}}) defining the folds of + the data for cross-validation. If supplied, this is + preferred over the arguments for generating + cross-validation folds.} + + \item{nCores}{a positive integer giving the number of + processor cores to be used for parallel computing in + cross-validation (the default is 1 for no + parallelization). If this is set to \code{NA}, all + available processor cores are used.} + + \item{cl}{a \pkg{parallel} cluster for parallel computing + to be used in cross-validation, as generated by + \code{\link[parallel]{makeCluster}}. If supplied, this + is preferred over \code{nCores}.} + \item{fallback}{logical; if a correlation functional other than the Pearson correlation is maximized, the data are first robustly standardized via median and MAD. This @@ -94,6 +118,23 @@ \item{b}{numeric; the weighting vector for \code{y}.} + \item{lambdaX}{a numeric giving the optimal penalty + parameter for the weighting vector of \code{x} (only + \code{sMaxCorGrid}).} + + \item{lambdaY}{a numeric giving the optimal penalty + parameter for the weighting vector of \code{y} (only + \code{sMaxCorGrid}).} + + \item{objective}{a numeric giving the value of the + objective function for the optimal combination of penalty + parameters (only \code{sMaxCorGrid}).} + + \item{cv}{a numeric matrix containing the (average) + results from cross-validation for each combination of + penalty parameters (only \code{sMaxCorGrid} and if a grid + of values for the penalty parameters has been supplied).} + \item{call}{the matched function call.} } \description{ @@ -103,7 +144,7 @@ with a focus on robust and nonparametric methods. } \details{ - The algorithm is based on alternating series of grid + \code{maxCorGrid} is based on alternating series of grid searches in two-dimensional subspaces of each data set. In each grid search, \code{nGrid} grid points on the unit circle in the corresponding plane are obtained, and the @@ -131,6 +172,16 @@ thereby possible to use randomly selected subsets of variables, or to specify the subsets of variables directly. + + \code{sMaxCorGrid} enforces sparsity by adding + \eqn{L_{1}}{L1} penalties on the weighting vectors to the + objective function (i.e., the correlation of the + resulting projections). The optimal combination of + penalty paramters is thereby determined via + cross-validation over a grid of values. To increase the + computational performance of this cross-validation + procedure, parallel computing is implemented via package + \pkg{parallel}. } \examples{ ## generate data @@ -145,11 +196,11 @@ x <- xy[, 1:p] y <- xy[, (p+1):m] ## Spearman correlation -maxCorGrid(x, y, method = "spearman") -maxCorGrid(x, y, method = "spearman", consistent = TRUE) +maxCorProj(x, y, method = "spearman") +maxCorProj(x, y, method = "spearman", consistent = TRUE) ## Pearson correlation -maxCorGrid(x, y, method = "pearson") +maxCorProj(x, y, method = "pearson") } \author{ Andreas Alfons From 1b828fbc59b59cd98a4c09d91fc1786478e746e4 Mon Sep 17 00:00:00 2001 From: Andreas Alfons Date: Wed, 13 Feb 2013 16:04:19 +0100 Subject: [PATCH 09/14] minimal cleanup --- R/cv.R | 4 ++-- src/cca.cpp | 1 - 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/R/cv.R b/R/cv.R index 39e2d2d..16a1987 100644 --- a/R/cv.R +++ b/R/cv.R @@ -72,8 +72,8 @@ corXY <- function(xy, corFun = corSpearman, corArgs = list()) { # compute the correlation between each pair of linear combinations if(p == 1) doCall(corFun, xy[, 1], xy[, 2], args=corArgs) else { - sapply(seq_len(p), function(i) { - doCall(corFun, xy[, i], xy[, p+i], args=corArgs) + sapply(seq_len(p), function(j) { + doCall(corFun, xy[, j], xy[, p+j], args=corArgs) }) } } diff --git a/src/cca.cpp b/src/cca.cpp index b1c81ee..14e8c30 100644 --- a/src/cca.cpp +++ b/src/cca.cpp @@ -737,7 +737,6 @@ void SparseGridControl::gridSearch(const mat& x, const uvec& orderX, // perform grid search for the current canonical basis vector vec corY(nGrid), objective(nGrid), penalties(nGrid); for(uword k = 0; k < nGrid; k++) { - double angle = grid(k); vec currentA = getVector(a, grid(k), orderJ); corY(k) = abs(corControl.cor(x * currentA, y)); penalties(k) = lambda * norm(currentA, 1); From d9b1221084ac6e6c10ec553c9b696760dd6d5450 Mon Sep 17 00:00:00 2001 From: Andreas Alfons Date: Thu, 14 Feb 2013 16:39:21 +0100 Subject: [PATCH 10/14] sparse maximum correlation estimators now have their own class --- NAMESPACE | 1 + R/maxCor.R | 2 +- R/print.R | 16 ++++++++++++++++ src/cca.cpp | 3 --- 4 files changed, 18 insertions(+), 4 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index a3361a8..264d97d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,6 +2,7 @@ S3method(print,cca) S3method(print,cvFolds) S3method(print,maxCor) S3method(print,permTest) +S3method(print,sMaxCor) export(CCAgrid) export(CCAproj) export(ccaGrid) diff --git a/R/maxCor.R b/R/maxCor.R index 05f87bc..1e6e092 100644 --- a/R/maxCor.R +++ b/R/maxCor.R @@ -383,7 +383,7 @@ sMaxCorPP <- function(x, y, method = c("spearman", "kendall", "quadrant", if(useCV) maxCor$cv <- cbind(lambda, CV=cv) } ## assign class and return results - class(maxCor) <- "maxCor" + class(maxCor) <- c("sMaxCor", "maxCor") maxCor } diff --git a/R/print.R b/R/print.R index df61bd5..370a8c8 100644 --- a/R/print.R +++ b/R/print.R @@ -67,3 +67,19 @@ print.permTest <- function(x, ...) { # return object invisibly invisible(x) } + +#' @S3method print sMaxCor +print.sMaxCor <- function(x, ...) { + print.maxCor(x, ...) + # print optimal tuning parameters + if(!is.null(cv <- x$cv)) { + haveGridX <- any(cv[, "lambdaX"] != 0) + haveGridY <- any(cv[, "lambdaY"] != 0) + if(haveGridX && haveGridY) cat("\nOptimal tuning parameters:\n") + else cat("\nOptimal tuning parameter:\n") + if(haveGridX) cat(sprintf("lambdaX: %f\n", x$lambdaX)) + if(haveGridY) cat(sprintf("lambdaY: %f\n", x$lambdaY)) + } + # return object invisibly + invisible(x) +} diff --git a/src/cca.cpp b/src/cca.cpp index 14e8c30..bcc9645 100644 --- a/src/cca.cpp +++ b/src/cca.cpp @@ -449,7 +449,6 @@ void GridControl::gridSearch(const mat& x, const uvec& orderX, const vec& y, // perform grid search for the current canonical basis vector vec corY(nGrid); for(uword k = 0; k < nGrid; k++) { - double angle = grid(k); vec currentA = getVector(a, grid(k), orderJ); corY(k) = abs(corControl.cor(x * currentA, y)); } @@ -1300,8 +1299,6 @@ template vec sMaxCorPP(const mat& x, const mat& y, CorControl& corControl, PPControl& ppControl, const bool& robust, const bool& fallback, mat& A, mat& B, vec& objective) { - // initializations - const uword p = x.n_cols, q = y.n_cols; // standardize the data vec scaleX, scaleY; mat xs = standardize(x, robust, fallback, scaleX); From c1c2fb1ff61f7a9fa5225c7094f31583c745a13f Mon Sep 17 00:00:00 2001 From: Andreas Alfons Date: Wed, 3 Apr 2013 17:24:11 +0200 Subject: [PATCH 11/14] updated maintainer e-mail --- DESCRIPTION | 14 ++------------ man/ccaPP-package.Rd | 4 ++-- 2 files changed, 4 insertions(+), 14 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index cd322ac..d9aaf02 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,7 +2,7 @@ Package: ccaPP Type: Package Title: (Robust) canonical correlation analysis via projection pursuit Version: 0.3.0 -Date: 2012-10-18 +Date: 2013-04-03 Depends: Rcpp (>= 0.9.10), parallel, @@ -17,15 +17,5 @@ Description: Canonical correlation analysis and maximum correlation via estimators, with a focus on robust and nonparametric methods. License: GPL (>= 2) LazyLoad: yes -Authors@R: c(person("Andreas", "Alfons", email = "andreas.alfons@kuleuven.be", +Authors@R: c(person("Andreas", "Alfons", email = "alfons@ese.eur.nl", role = c("aut", "cre")), person("David", "Simcha", role = "ctb")) -Collate: - 'cca.R' - 'cor.R' - 'cv.R' - 'fastMAD.R' - 'fastMedian.R' - 'maxCor.R' - 'permTest.R' - 'print.R' - 'utils.R' diff --git a/man/ccaPP-package.Rd b/man/ccaPP-package.Rd index 80f64e7..72f33d2 100644 --- a/man/ccaPP-package.Rd +++ b/man/ccaPP-package.Rd @@ -15,7 +15,7 @@ Canonical correlation analysis and maximum correlation via Package: \tab ccaPP\cr Type: \tab Package\cr Version: \tab 0.3.0\cr -Date: \tab 2012-10-18\cr +Date: \tab 2013-04-03\cr Depends: \tab Rcpp (>= 0.9.10), parallel, pcaPP (>= 1.8-1)\cr Imports: \tab RcppArmadillo (>= 0.2.36)\cr LinkingTo: \tab Rcpp, RcppArmadillo\cr @@ -48,6 +48,6 @@ permTest (Robust) permutation test for independence Andreas Alfons [aut, cre], David Simcha [ctb] -Maintainer: Andreas Alfons +Maintainer: Andreas Alfons } \keyword{package} From 927e9f67523c0a17433c87dde150724858d05e4c Mon Sep 17 00:00:00 2001 From: Andreas Alfons Date: Sun, 23 Feb 2014 10:46:27 +0100 Subject: [PATCH 12/14] optional standardization in sMaxCorGrid() --- R/maxCor.R | 52 +++++++++-------- man/maxCorGrid.Rd | 11 ++-- src/cca.cpp | 143 ++++++++++++++++++++++++++-------------------- src/cca.h | 2 +- 4 files changed, 114 insertions(+), 94 deletions(-) diff --git a/R/maxCor.R b/R/maxCor.R index dfb3c3f..4f8708a 100644 --- a/R/maxCor.R +++ b/R/maxCor.R @@ -146,11 +146,11 @@ #' y <- xy[, (p+1):m] #' #' ## Spearman correlation -#' maxCorProj(x, y, method = "spearman") -#' maxCorProj(x, y, method = "spearman", consistent = TRUE) +#' maxCorGrid(x, y, method = "spearman") +#' maxCorGrid(x, y, method = "spearman", consistent = TRUE) #' #' ## Pearson correlation -#' maxCorProj(x, y, method = "pearson") +#' maxCorGrid(x, y, method = "pearson") #' #' @keywords multivariate robust #' @@ -161,10 +161,9 @@ maxCorGrid <- function(x, y, method = c("spearman", "kendall", "quadrant", "M", "pearson"), - control = list(...), nIterations = 10, - nAlternate = 10, nGrid = 25, select = NULL, - tol = 1e-06, standardize = TRUE, fallback = FALSE, - seed = NULL, ...) { + control = list(...), nIterations = 10, nAlternate = 10, + nGrid = 25, select = NULL, tol = 1e-06, + standardize = TRUE, fallback = FALSE, seed = NULL, ...) { ## initializations matchedCall <- match.call() ## define list of control arguments for algorithm @@ -188,13 +187,13 @@ maxCorGrid <- function(x, y, #' @export sMaxCorGrid <- function(x, y, lambdaX = 0, lambdaY = 0, - method = c("spearman", "kendall", "quadrant", - "M", "pearson"), + method = c("spearman", "kendall", "quadrant", "M", "pearson"), control = list(...), nIterations = 10, nAlternate = 10, - nGrid = 25, select = NULL, tol = 1e-06, K = 5, R = 1, + nGrid = 25, select = NULL, tol = 1e-06, + standardize = TRUE, fallback = FALSE, K = 5, R = 1, type = c("random", "consecutive", "interleaved"), grouping = NULL, folds = NULL, nCores = 1, cl = NULL, - fallback = FALSE, seed = NULL, ...) { + seed = NULL, ...) { ## initializations matchedCall <- match.call() ## define list of control arguments for algorithm @@ -210,9 +209,10 @@ sMaxCorGrid <- function(x, y, lambdaX = 0, lambdaY = 0, nAlternate=nAlternate, nGrid=nGrid, select=select, tol=tol) ## call workhorse function maxCor <- sMaxCorPP(x, y, method=method, corControl=control, - algorithm="grid", ppControl=ppControl, K=K, + algorithm="grid", ppControl=ppControl, + standardize=standardize, fallback=fallback, K=K, R=R, type=type, grouping=grouping, folds=folds, - nCores=nCores, cl=cl, fallback=fallback, seed=seed) + nCores=nCores, cl=cl, seed=seed) maxCor$call <- matchedCall maxCor } @@ -335,11 +335,11 @@ maxCorPP <- function(x, y, ...) { ## workhorse function for sparse maximum correlation #' @import parallel -sMaxCorPP <- function(x, y, method = c("spearman", "kendall", "quadrant", - "M", "pearson"), - corControl, algorithm = "grid", ppControl, folds = NULL, - nCores = 1, cl = NULL, fallback = FALSE, seed = NULL, - ...) { +sMaxCorPP <- function(x, y, + method = c("spearman", "kendall", "quadrant", "M", "pearson"), + corControl, algorithm = "grid", ppControl, + standardize = TRUE, fallback = FALSE, folds = NULL, + nCores = 1, cl = NULL, seed = NULL, ...) { ## initializations x <- as.matrix(x) y <- as.matrix(y) @@ -356,6 +356,7 @@ sMaxCorPP <- function(x, y, method = c("spearman", "kendall", "quadrant", # get list of control arguments method <- match.arg(method) corControl <- getCorControl(method, corControl, forceConsistency=FALSE) + standardize <- isTRUE(standardize) fallback <- isTRUE(fallback) if(!is.null(seed)) set.seed(seed) ppControl <- getPPControl(algorithm, ppControl, p=p, q=q) @@ -396,7 +397,8 @@ sMaxCorPP <- function(x, y, method = c("spearman", "kendall", "quadrant", } # perform cross-validation over grid of tuning parameters call <- call("sMaxCorPPFit", method=method, corControl=corControl, - algorithm=algorithm, ppControl=ppControl, fallback=fallback) + algorithm=algorithm, ppControl=ppControl, + standardize=standardize, fallback=fallback) call[[1]] <- sMaxCorPPFit # necessary to work with parallel computing corFun <- switch(method, spearman=corSpearman, kendall=corKendall, quadrant=corQuadrant, M=corM, pearson=corPearson) @@ -408,7 +410,7 @@ sMaxCorPP <- function(x, y, method = c("spearman", "kendall", "quadrant", # compute the final solution maxCor <- sMaxCorPPFit(x, y, method=method, corControl=corControl, algorithm=algorithm, ppControl=ppControl, - fallback=fallback) + standardize=standardize, fallback=fallback) maxCor$a <- drop(maxCor$a) maxCor$b <- drop(maxCor$b) if(useCV) maxCor$cv <- cbind(lambda, CV=cv) @@ -419,13 +421,13 @@ sMaxCorPP <- function(x, y, method = c("spearman", "kendall", "quadrant", } ## simple wrapper for calling the C++ function -sMaxCorPPFit <- function(x, y, method = c("spearman", "kendall", "quadrant", - "M", "pearson"), +sMaxCorPPFit <- function(x, y, + method = c("spearman", "kendall", "quadrant", "M", "pearson"), corControl, algorithm = "grid", ppControl, - fallback = FALSE) { + standardize = TRUE, fallback = FALSE) { # call C++ function maxCor <- .Call("R_sMaxCorPP", R_x=x, R_y=y, R_method=method, R_corControl=corControl, R_algorithm=algorithm, - R_ppControl=ppControl, R_fallback=fallback, - PACKAGE="ccaPP") + R_ppControl=ppControl, R_standardize=standardize, + R_fallback=fallback, PACKAGE="ccaPP") } diff --git a/man/maxCorGrid.Rd b/man/maxCorGrid.Rd index c2e7b89..5f9a115 100644 --- a/man/maxCorGrid.Rd +++ b/man/maxCorGrid.Rd @@ -13,10 +13,11 @@ sMaxCorGrid(x, y, lambdaX = 0, lambdaY = 0, method = c("spearman", "kendall", "quadrant", "M", "pearson"), control = list(...), nIterations = 10, nAlternate = 10, - nGrid = 25, select = NULL, tol = 1e-06, K = 5, R = 1, + nGrid = 25, select = NULL, tol = 1e-06, + standardize = TRUE, fallback = FALSE, K = 5, R = 1, type = c("random", "consecutive", "interleaved"), grouping = NULL, folds = NULL, nCores = 1, cl = NULL, - fallback = FALSE, seed = NULL, ...) + seed = NULL, ...) } \arguments{ \item{x,y}{each can be a numeric vector, matrix or data @@ -215,11 +216,11 @@ x <- xy[, 1:p] y <- xy[, (p+1):m] ## Spearman correlation -maxCorProj(x, y, method = "spearman") -maxCorProj(x, y, method = "spearman", consistent = TRUE) +maxCorGrid(x, y, method = "spearman") +maxCorGrid(x, y, method = "spearman", consistent = TRUE) ## Pearson correlation -maxCorProj(x, y, method = "pearson") +maxCorGrid(x, y, method = "pearson") } \author{ Andreas Alfons diff --git a/src/cca.cpp b/src/cca.cpp index c3112f1..4eff150 100644 --- a/src/cca.cpp +++ b/src/cca.cpp @@ -1353,84 +1353,101 @@ SEXP R_ccaPP(SEXP R_x, SEXP R_y, SEXP R_k, SEXP R_method, SEXP R_corControl, // y ............ second data matrix // corControl ... control object to compute correlation // ppControl .... control object for sparse algorithm -// robust ....... should the data be robustly standardized? +// standard ..... should the data be standardized? +// robust ....... should robust standardization be used? // fallback ..... should the fallback mode for robust standardization be used? // A ............ matrix of weighting vectors for first matrix to be updated // B ............ matrix of weighting vectors for second matrix to be updated template vec sMaxCorPP(const mat& x, const mat& y, CorControl& corControl, - PPControl& ppControl, const bool& robust, const bool& fallback, - mat& A, mat& B, vec& objective) { - // standardize the data - vec scaleX, scaleY; - mat xs = standardize(x, robust, fallback, scaleX); - mat ys = standardize(y, robust, fallback, scaleY); - // compute maximum correlations with standardized data - vec r = ppControl.maxCor(xs, ys, corControl, A, B, objective); - // transform weighting vectors back to original scale - for(uword k = 0; k < A.n_cols; k++) { - vec a = A.unsafe_col(k); backtransform(a, scaleX); - } - for(uword k = 0; k < B.n_cols; k++) { - vec b = B.unsafe_col(k); backtransform(b, scaleY); - } + PPControl& ppControl, const bool& standard, const bool& robust, + const bool& fallback, mat& A, mat& B, vec& centerX, vec& centerY, + vec& scaleX, vec& scaleY, vec& objective) { + // initializations + vec r; + // standardize the data if requested + mat xs, ys; + if(standard) { + xs = standardize(x, robust, fallback, centerX, scaleX); + ys = standardize(y, robust, fallback, centerY, scaleY); + // compute maximum correlations with standardized data + r = ppControl.maxCor(xs, ys, corControl, A, B, objective); + // transform weighting vectors back to original scale + for(uword k = 0; k < A.n_cols; k++) { + vec a = A.unsafe_col(k); backtransform(a, scaleX); + } + for(uword k = 0; k < B.n_cols; k++) { + vec b = B.unsafe_col(k); backtransform(b, scaleY); + } + } else { + uword p = x.n_cols, q = y.n_cols; + centerX = zeros(p); centerY = zeros(q); + scaleX = ones(p); scaleY = ones(q); + // compute maximum correlations with original data + r = ppControl.maxCor(x, y, corControl, A, B, objective); + } // return vector of maximum correlations - return r; + return r; } // R interface SEXP R_sMaxCorPP(SEXP R_x, SEXP R_y, SEXP R_method, SEXP R_corControl, - SEXP R_algorithm, SEXP R_ppControl, SEXP R_fallback) { + SEXP R_algorithm, SEXP R_ppControl, SEXP R_standardize, SEXP R_fallback) { // initializations NumericMatrix Rcpp_x(R_x), Rcpp_y(R_y); // convert data to Rcpp types mat x(Rcpp_x.begin(), Rcpp_x.nrow(), Rcpp_x.ncol(), false); // convert data mat y(Rcpp_y.begin(), Rcpp_y.nrow(), Rcpp_y.ncol(), false); // to arma types - string method = as(R_method); // convert character string - List Rcpp_corControl(R_corControl); // list of control parameters - string algorithm = as(R_algorithm); // convert character string - List Rcpp_ppControl(R_ppControl); // list of control parameters - bool fallback = as(R_fallback); // convert to boolean - // initialize results - vec r, lambdaX, lambdaY, objective; - mat A, B; + string method = as(R_method); // convert character string + List Rcpp_corControl(R_corControl); // list of control parameters + string algorithm = as(R_algorithm); // convert character string + List Rcpp_ppControl(R_ppControl); // list of control parameters + bool standard = as(R_standardize); // convert to boolean + bool fallback = as(R_fallback); // convert to boolean + // initialize results + vec r, centerX, centerY, scaleX, scaleY, lambdaX, lambdaY, objective; + mat A, B; if(algorithm == "grid") { - // define control object for alternate grid searches - SparseGridControl ppControl(Rcpp_ppControl); - // define control object for the correlations and call the arma version - if(method == "spearman") { - CorSpearmanControl corControl(Rcpp_corControl); - r = sMaxCorPP(x, y, corControl, ppControl, true, fallback, - A, B, objective); - } else if(method == "kendall") { - CorKendallControl corControl(Rcpp_corControl); - r = sMaxCorPP(x, y, corControl, ppControl, true, fallback, - A, B, objective); - } else if(method == "quadrant") { - CorQuadrantControl corControl(Rcpp_corControl); - r = sMaxCorPP(x, y, corControl, ppControl, true, fallback, - A, B, objective); - } else if(method == "M") { - CorMControl corControl(Rcpp_corControl); - r = sMaxCorPP(x, y, corControl, ppControl, true, fallback, - A, B, objective); - } else if(method == "pearson") { - CorPearsonControl corControl; - r = sMaxCorPP(x, y, corControl, ppControl, false, false, - A, B, objective); - } else { - error("method not available"); - } + // define control object for alternate grid searches + SparseGridControl ppControl(Rcpp_ppControl); + // define control object for the correlations and call the arma version + if(method == "spearman") { + CorSpearmanControl corControl(Rcpp_corControl); + r = sMaxCorPP(x, y, corControl, ppControl, standard, true, fallback, + A, B, centerX, centerY, scaleX, scaleY, objective); + } else if(method == "kendall") { + CorKendallControl corControl(Rcpp_corControl); + r = sMaxCorPP(x, y, corControl, ppControl, standard, true, fallback, + A, B, centerX, centerY, scaleX, scaleY, objective); + } else if(method == "quadrant") { + CorQuadrantControl corControl(Rcpp_corControl); + r = sMaxCorPP(x, y, corControl, ppControl, standard, true, fallback, + A, B, centerX, centerY, scaleX, scaleY, objective); + } else if(method == "M") { + CorMControl corControl(Rcpp_corControl); + r = sMaxCorPP(x, y, corControl, ppControl, standard, true, fallback, + A, B, centerX, centerY, scaleX, scaleY, objective); + } else if(method == "pearson") { + CorPearsonControl corControl; + r = sMaxCorPP(x, y, corControl, ppControl, standard, false, false, + A, B, centerX, centerY, scaleX, scaleY, objective); + } else { + error("method not available"); + } lambdaX = ppControl.lambda.col(0); lambdaY = ppControl.lambda.col(1); - } else { - error("algorithm not available"); - } - // wrap and return result - return List::create( - Named("cor") = wrap(r.begin(), r.end()), - Named("a") = A, - Named("b") = B, - Named("lambdaX") = wrap(lambdaX.begin(), lambdaX.end()), - Named("lambdaY") = wrap(lambdaY.begin(), lambdaY.end()), - Named("objective") = wrap(objective.begin(), objective.end()) + } else { + error("algorithm not available"); + } + // wrap and return result + return List::create( + Named("cor") = wrap(r.begin(), r.end()), + Named("a") = A, + Named("b") = B, + Named("centerX") = centerX, + Named("centerY") = centerY, + Named("scaleX") = scaleX, + Named("scaleY") = scaleY, + Named("lambdaX") = wrap(lambdaX.begin(), lambdaX.end()), + Named("lambdaY") = wrap(lambdaY.begin(), lambdaY.end()), + Named("objective") = wrap(objective.begin(), objective.end()) ); } diff --git a/src/cca.h b/src/cca.h index 51a13b1..53f65fb 100644 --- a/src/cca.h +++ b/src/cca.h @@ -18,6 +18,6 @@ RcppExport SEXP R_fastCor(SEXP R_x, SEXP R_y, SEXP R_method, SEXP R_control); // RcppExport SEXP R_ccaPP(SEXP R_x, SEXP R_y, SEXP R_k, SEXP R_method, SEXP R_corControl, SEXP R_algorithm, SEXP R_ppControl, SEXP R_standardize, SEXP R_fallback); RcppExport SEXP R_sMaxCorPP(SEXP R_x, SEXP R_y, SEXP R_method, SEXP R_corControl, - SEXP R_algorithm, SEXP R_ppControl, SEXP R_fallback); + SEXP R_algorithm, SEXP R_ppControl, SEXP R_standardize, SEXP R_fallback); #endif From c44fe041252ff8819a259ac32cad581289ffbc09 Mon Sep 17 00:00:00 2001 From: Andreas Alfons Date: Sun, 23 Feb 2014 10:53:52 +0100 Subject: [PATCH 13/14] updated author information --- R/cca.R | 6 +++--- R/cor.R | 6 +++--- R/cv.R | 6 +++--- R/fastMAD.R | 6 +++--- R/fastMedian.R | 6 +++--- R/maxCor.R | 6 +++--- R/permTest.R | 6 +++--- R/print.R | 6 +++--- R/utils.R | 6 +++--- src/cca.cpp | 2 +- src/cca.h | 2 +- src/cor.cpp | 2 +- src/cor.h | 2 +- src/fastCorKendall.cpp | 2 +- src/fastCorKendall.h | 2 +- src/utils.cpp | 2 +- src/utils.h | 2 +- 17 files changed, 35 insertions(+), 35 deletions(-) diff --git a/R/cca.R b/R/cca.R index a0a690a..88df262 100644 --- a/R/cca.R +++ b/R/cca.R @@ -1,7 +1,7 @@ -# ------------------------------------ +# -------------------------------------- # Author: Andreas Alfons -# Erasmus University Rotterdam -# ------------------------------------ +# Erasmus Universiteit Rotterdam +# -------------------------------------- #' (Robust) CCA via alternating series of grid searches #' diff --git a/R/cor.R b/R/cor.R index 93109dc..f0dce9b 100644 --- a/R/cor.R +++ b/R/cor.R @@ -1,7 +1,7 @@ -# ---------------------- +# -------------------------------------- # Author: Andreas Alfons -# KU Leuven -# ---------------------- +# Erasmus Universiteit Rotterdam +# -------------------------------------- #' Fast implementations of (robust) correlation estimators #' diff --git a/R/cv.R b/R/cv.R index 16a1987..03ee1d9 100644 --- a/R/cv.R +++ b/R/cv.R @@ -1,7 +1,7 @@ -# ---------------------- +# -------------------------------------- # Author: Andreas Alfons -# KU Leuven -# ---------------------- +# Erasmus Universiteit Rotterdam +# -------------------------------------- ## perform cross-validation cvMaxCor <- function(call, x, y, folds, corFun = corSpearman, corArgs = list(), diff --git a/R/fastMAD.R b/R/fastMAD.R index 1ebe5ad..2a120e2 100644 --- a/R/fastMAD.R +++ b/R/fastMAD.R @@ -1,7 +1,7 @@ -# ---------------------- +# -------------------------------------- # Author: Andreas Alfons -# KU Leuven -# ---------------------- +# Erasmus Universiteit Rotterdam +# -------------------------------------- #' Fast implementation of the median absolute deviation #' diff --git a/R/fastMedian.R b/R/fastMedian.R index ba94caa..6015448 100644 --- a/R/fastMedian.R +++ b/R/fastMedian.R @@ -1,7 +1,7 @@ -# ---------------------- +# -------------------------------------- # Author: Andreas Alfons -# KU Leuven -# ---------------------- +# Erasmus Universiteit Rotterdam +# -------------------------------------- #' Fast implementation of the median #' diff --git a/R/maxCor.R b/R/maxCor.R index 4f8708a..c16e623 100644 --- a/R/maxCor.R +++ b/R/maxCor.R @@ -1,7 +1,7 @@ -# ------------------------------------ +# -------------------------------------- # Author: Andreas Alfons -# Erasmus University Rotterdam -# ------------------------------------ +# Erasmus Universiteit Rotterdam +# -------------------------------------- #' (Robust) maximum correlation via alternating series of grid searches #' diff --git a/R/permTest.R b/R/permTest.R index 87da117..553b70a 100644 --- a/R/permTest.R +++ b/R/permTest.R @@ -1,7 +1,7 @@ -# ---------------------- +# -------------------------------------- # Author: Andreas Alfons -# KU Leuven -# ---------------------- +# Erasmus Universiteit Rotterdam +# -------------------------------------- #' (Robust) permutation test for independence #' diff --git a/R/print.R b/R/print.R index 370a8c8..ee7013f 100644 --- a/R/print.R +++ b/R/print.R @@ -1,7 +1,7 @@ -# ---------------------- +# -------------------------------------- # Author: Andreas Alfons -# KU Leuven -# ---------------------- +# Erasmus Universiteit Rotterdam +# -------------------------------------- #' @S3method print cca print.cca <- function(x, ...) { diff --git a/R/utils.R b/R/utils.R index d603d23..a65f277 100644 --- a/R/utils.R +++ b/R/utils.R @@ -1,7 +1,7 @@ -# ---------------------- +# -------------------------------------- # Author: Andreas Alfons -# KU Leuven -# ---------------------- +# Erasmus Universiteit Rotterdam +# -------------------------------------- ### transform canonical vectors back to original scale #backtransform <- function(A, scale) { diff --git a/src/cca.cpp b/src/cca.cpp index 4eff150..8aae800 100644 --- a/src/cca.cpp +++ b/src/cca.cpp @@ -1,6 +1,6 @@ /* * Author: Andreas Alfons - * Erasmus University Rotterdam + * Erasmus Universiteit Rotterdam */ #include diff --git a/src/cca.h b/src/cca.h index 53f65fb..06a4c4a 100644 --- a/src/cca.h +++ b/src/cca.h @@ -1,6 +1,6 @@ /* * Author: Andreas Alfons - * Erasmus University Rotterdam + * Erasmus Universiteit Rotterdam */ #ifndef _ccaPP_CCA_H diff --git a/src/cor.cpp b/src/cor.cpp index 27eff5d..c144901 100644 --- a/src/cor.cpp +++ b/src/cor.cpp @@ -1,6 +1,6 @@ /* * Author: Andreas Alfons - * KU Leuven + * Erasmus Universiteit Rotterdam */ #include diff --git a/src/cor.h b/src/cor.h index 8d646fc..22179a4 100644 --- a/src/cor.h +++ b/src/cor.h @@ -1,6 +1,6 @@ /* * Author: Andreas Alfons - * KU Leuven + * Erasmus Universiteit Rotterdam */ #ifndef _ccaPP_COR_H diff --git a/src/fastCorKendall.cpp b/src/fastCorKendall.cpp index 3be9dc4..537f9aa 100755 --- a/src/fastCorKendall.cpp +++ b/src/fastCorKendall.cpp @@ -36,7 +36,7 @@ /* * Modifications: Andreas Alfons - * KU Leuven + * Erasmus Universiteit Rotterdam */ diff --git a/src/fastCorKendall.h b/src/fastCorKendall.h index 1dd033a..405dd9a 100644 --- a/src/fastCorKendall.h +++ b/src/fastCorKendall.h @@ -1,6 +1,6 @@ /* * Author: Andreas Alfons - * KU Leuven + * Erasmus Universiteit Rotterdam */ #ifndef _ccaPP_FASTCORKENDALL_H diff --git a/src/utils.cpp b/src/utils.cpp index 83470b6..0255f8f 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -1,6 +1,6 @@ /* * Author: Andreas Alfons - * Erasmus University Rotterdam + * Erasmus Universiteit Rotterdam */ #include diff --git a/src/utils.h b/src/utils.h index 2186d60..53c84e8 100644 --- a/src/utils.h +++ b/src/utils.h @@ -1,6 +1,6 @@ /* * Author: Andreas Alfons - * Erasmus University Rotterdam + * Erasmus Universiteit Rotterdam */ #ifndef _ccaPP_UTILS_H From 609722eda7728318a019e13793be07d6586c089c Mon Sep 17 00:00:00 2001 From: Andreas Alfons Date: Tue, 19 Jun 2018 12:58:15 +0200 Subject: [PATCH 14/14] minor clean up in references --- inst/CITATION | 2 +- vignettes/maxCor.bib | 21 +++++++++++---------- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/inst/CITATION b/inst/CITATION index e4f84fe..13e23e2 100644 --- a/inst/CITATION +++ b/inst/CITATION @@ -32,5 +32,5 @@ citEntry(entry = "Article", paste("Andreas Alfons, Christophe Croux, Peter Filzmoser (2017).", "Robust maximum association estimators.", "Journal of the American Statistical Association,", - "112(517), 436--445.") + "112(517), 436-445.") ) diff --git a/vignettes/maxCor.bib b/vignettes/maxCor.bib index 7af734d..44306bd 100644 --- a/vignettes/maxCor.bib +++ b/vignettes/maxCor.bib @@ -36,7 +36,7 @@ @article{CroR05 @ARTICLE{DevG75, AUTHOR = {Devlin, S.J. and Gnanadesikan, R. and Kettenring, J.R.}, - TITLE = {Robust Estimation and Outlier Detection with Correlation + TITLE = {Robust Estimation and Outlier Detection with Correlation Coefficients}, JOURNAL = {Biometrika}, YEAR = {1975}, @@ -57,7 +57,7 @@ @ARTICLE{CroH00 @ARTICLE{CroD10, AUTHOR = {Croux, C. and Dehon, C.}, - TITLE = {Influence functions of the {Spearman} and {Kendall} correlation + TITLE = {Influence functions of the {Spearman} and {Kendall} correlation measures}, JOURNAL = {Statistical Methods \& Applications}, YEAR = {2010}, @@ -68,7 +68,7 @@ @ARTICLE{CroD10 @ARTICLE{CuiH03, AUTHOR = {Cui, H. and He, X. and Ng, K.W.}, YEAR = 2003, - TITLE = {Asymptotic distributions of principal components based on + TITLE = {Asymptotic distributions of principal components based on robust dispersions}, JOURNAL = {Biometrika}, VOLUME = {90}, @@ -188,7 +188,7 @@ @article{TasK03 journal = {Statistics and Probability Letters}, volume = {62}, number = {1}, - pages = {9--21}, + pages = {9--21}, year = {2003} } @@ -250,7 +250,7 @@ @manual{RDev } @article{AlfT10, - title = {An object-oriented framework for statistical simulation: The + title = {An object-oriented framework for statistical simulation: The \proglang{R} package \pkg{simFrame}}, author = {Alfons, A. and Templ, M. and Filzmoser, P.}, journal = {Journal of Statistical Software}, @@ -272,9 +272,10 @@ @article{AlfCXX author = {Alfons, A. and Croux, C. and Filzmoser, P.}, title = {Robust Maximum Association Estimators}, journal = {Journal of the American Statistical Association}, - year = {2016}, - note = {In press}, - doi = {10.1080/01621459.2016.1148609}, + volume = {112}, + number = {517}, + pages = {436--445}, + year = {2017} } @article{alfons16b, @@ -427,7 +428,7 @@ @article{ParT09 @article{WaaV08, author = {Waaijenborg, S. and {Verselewel~de~Witt~Hamer}, P.C. and Zwinderman, A.H.}, title = {Quantifying the association between gene expressions and {DNA} markers by penalized canonical correlation analysis}, - journal = {Statistical Applications in Genetics and Molecular Biology}, + journal = {Statistical Applications in Genetics and Molecular Biology}, volume = {7}, year = {2008} } @@ -453,7 +454,7 @@ @article{WitT09 @techreport{WilC13, author = {Wilms, I. and Croux, C.}, title = {Sparse Canonical Correlation Analysis from a Predictive Point of View}, - type = {FEB Research Report}, + type = {FEB Research Report}, number = {KBI\_1320}, institution = {Faculty of Economics and Business, KU Leuven}, year = {2013}