Catching LAPACK errors in Armadillo

Question

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:

OS: Red Hat Enterprise Linux
R version: 4.0.2 (2020-06-22)
Rcpp version: 1.0.5
RcppArmadillo version: 0.12.2.0.0
The problematic matrix: stored in this RData file

score 2 · Accepted Answer · answered May 22 '23 at 12:44

2

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()).

Code

#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)
*/

Demo

> Rcpp::sourceCpp("answer.cpp")

> load("bad_svd.RData")

> call_cond(draw_cov)
[1] 34.3382
>

Comment

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.

answered May 22 '23 at 12:44

Dirk Eddelbuettel

360,940
56
644
725

Thanks Dirk. Can you clarify what you mean by "the config file of `#define` values"? – Wilbur May 23 '23 at 09:10
Your code loads a file `RcppArmadillo.h` as an include. If you look at that file, you will see other files include by it and eg one named `RcppArmadilloConfig.h`. We may have set a preference for exceptions there with a value. You can override this with a _local `#define ...`_ you make before including `RcppArmadillo.h`. There are likely examples here at SO from when we had to, say, turn 64bit indexing on as an option or when we select the updated return type. You should be able to find those via search here. – Dirk Eddelbuettel May 23 '23 at 12:07
Or, of course, do as my answer shows you how to: don't change anything, but check the condition number before you would call a `svd()` operation that would throw. – Dirk Eddelbuettel May 23 '23 at 12:58

Catching LAPACK errors in Armadillo

1 Answers1

Code

Demo

Comment