# Performance of the divide-and-conquer SVD algorithm

December 9, 2013
By

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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::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. R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.