I'm attempting to calculate the SVD of a matrix using RcppArmadillo.svd() is supposed to return 0 if the SVD failed. However, I have encountered a matrix for which I get the error code BLAS/LAPACK routine 'DLASCL' gave error code -4.
I think the bug is fixed in updated packages (I'm not getting it when I run the code on my local machine) but I can't update the version of R on the server I need to use. As such, I would like a manual work-around. However I'm not sure how to capture the LAPACK error in my RCPP code. In particular, this code:
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace arma;
// [[Rcpp::export]]
void test_svd(mat X) {
mat U;
vec s;
mat V;
Rprintf("Program starting \n");
try {
Rprintf("... attempting to find the svd with the dc method \n");
svd(U, s, V, X, "dc");
}
catch (...) {
Rprintf("... that failed. \n");
}
}
throws an error on the problematic matrix, rather than running the catch block.
Is it possible to write a C++ block that would catch the LAPACK error?
Additional details:
As the arma documentation site shows, while svd(X) throws you can also access alternate methods such as (member function) svd(s, X) should not be throwing an exception. We may have a configuration default that prefers a throw: for R that is actually useful. You could check the config file of #define values, maybe there is something to override.
Besides that you can of course inquire about the 'condition number' using cond() before you call svd() (as well as the documented faster rcond()).
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
double call_cond(arma::mat X) {
double val = arma::cond(X);
return val;
}
/*** R
load("bad_svd.RData")
call_cond(draw_cov)
*/
> Rcpp::sourceCpp("answer.cpp")
> load("bad_svd.RData")
> call_cond(draw_cov)
[1] 34.3382
>
The value is fairly far from the ideal values at or near one so it give you a strong hint the matrix may not be singular and SVD could fail.