英文:
Catching LAPACK errors in Armadillo
问题
I'm here to provide the translation without any additional content. Here is the translation of the provided text:
我试图使用RcppArmadillo计算矩阵的奇异值分解（SVD）。svd（）应该在SVD失败时返回0。然而，我遇到了一个矩阵，对于它，我得到了错误代码BLAS/LAPACK例程' DLASCL'给出了错误代码-4。
我认为这个错误在更新的软件包中已经修复（当我在本地计算机上运行代码时没有出现），但我无法更新需要使用的服务器上的R版本。因此，我想要一个手动的解决方法。但是，我不确定如何在我的RCPP代码中捕获LAPACK错误。特别是，这段代码：
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace arma;
// [[Rcpp::export]]
void test_svd(mat X) {
mat U;
vec s;
mat V;
Rprintf("程序开始 \n");
try {
Rprintf("... 尝试使用dc方法找到svd \n");
svd(U, s, V, X, "dc");
}
catch (...) {
Rprintf("... 失败了。 \n");
}
}
在有问题的矩阵上引发错误，而不是运行catch块。
是否可能编写一个C++块来捕获LAPACK错误？
其他细节：
- 操作系统：Red Hat Enterprise Linux
- R版本：4.0.2（2020-06-22）
- Rcpp版本：1.0.5
- RcppArmadillo版本：0.12.2.0.0
- 有问题的矩阵：存储在此RData文件中
英文:
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
答案1
得分: 2
根据arma文档网站显示，svd(X)可能会抛出异常，但也可以使用替代方法（成员函数）svd(s, X)，不应该抛出异常。我们可能有一个配置默认项，它更倾向于抛出异常：对于R来说，这实际上是有用的。你可以检查#define值的配置文件，也许有一些可以覆盖的内容。
此外，在调用svd()之前，你当然可以使用cond()查询'condition number'（条件数）（以及文档中提到的更快的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
>
评论
这个值与理想值（接近或等于1）相差很远，所以它给出了一个强烈的提示，矩阵可能不是奇异的，SVD可能会失败。
英文:
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.
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论