**Rcpp Gallery**, and kindly contributed to R-bloggers)

The ubiquitous LAPACK library provides several implementations for the singular-value decomposition (SVD). We will illustrate possible speed gains from using the divide-and-conquer method by comparing it to the base case.

`#include `
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
arma::vec baseSVD(const arma::mat & X) {
arma::mat U, V;
arma::vec S;
arma::svd(U, S, V, X, "standard");
return S;
}
// [[Rcpp::export]]
arma::vec dcSVD(const arma::mat & X) {
arma::mat U, V;
arma::vec S;
arma::svd(U, S, V, X, "dc");
return S;
}

Having the two implementations, which differ only in the `method`

argument (added recently in Armadillo 3.930), we are ready to do a simple timing comparison.

```
library(microbenchmark)
set.seed(42)
X <- matrix(rnorm(16e4), 4e2, 4e2)
microbenchmark(baseSVD(X), dcSVD(X))
```

Unit: milliseconds expr min lq median uq max neval baseSVD(X) 421.2 422.6 424.2 426.2 442.1 100 dcSVD(X) 111.0 111.5 111.9 113.6 126.1 100

The speed gain is noticeable which a ratio of about 3.9 at the median. However, we can also look at complex-valued matrices.

```
// [[Rcpp::export]]
arma::vec cxBaseSVD(const arma::cx_mat & X) {
arma::cx_mat U, V;
arma::vec S;
arma::svd(U, S, V, X, "standard");
return S;
}
// [[Rcpp::export]]
arma::vec cxDcSVD(const arma::cx_mat & X) {
arma::cx_mat U, V;
arma::vec S;
arma::svd(U, S, V, X, "dc");
return S;
}
```

```
A <- matrix(rnorm(16e4), 4e2, 4e2)
B <- matrix(rnorm(16e4), 4e2, 4e2)
X <- A + 1i * B
microbenchmark(cxBaseSVD(X), cxDcSVD(X))
```

Unit: milliseconds expr min lq median uq max neval cxBaseSVD(X) 1248.7 1253.7 1257.5 1262.3 1311.7 100 cxDcSVD(X) 259.2 259.8 260.5 263.2 327.9 100

Here the difference is even more pronounced at about 4.8. However, it is precisely this complex-value divide-and-conquer method which is missing in R’s own Lapack. So R built with the default configuration can not currently take advantage of the complex-valued divide-and-conquer algorithm. Only builds which use an external Lapack library (as for example the Debian and Ubuntu builds) can. Let’s hope that R will add this functionality to its next release R 3.1.0. *Update: And the underlying `zgesdd`
routine has now been added to the upcoming R 3.1.0 release. Nice.*

**leave a comment**for the author, please follow the link and comment on their blog:

**Rcpp Gallery**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...