aen*_*ima 4 c++ r matrix armadillo rcpp
我有一个病态的4x4矩阵,使expmatArmadillo 的功能挂起.病理矩阵是:
a<-matrix(c(-2.5654e+060,4.6979e-018,2.5654e+060,7.2765e-035
,2.8913e+000, -3.6633e+001,3.3731e+001,1.0003e-002
,1.0656e-009,1.9037e-002, -1.9732e-001,1.7828e-001
, 0e+000, 0e+000, 0e+000, 0e+000), nrow=4, byrow=T)
.cpp文件如下所示:
# include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;
// [[Rcpp::export]]
mat exp_mat(mat x) {
return(expmat(x));
}
将病理矩阵提供给此函数会使其挂起并显示一条消息:
warning: solve(): system seems singular; attempting approx solution
我知道这个矩阵条件很差,但expmR包"expm"中的函数可以使用它的默认算法处理它没有问题.无论如何在RcppArmadillo解决这个问题？至少我想通过处理警告信息来避免挂起.
还有一个类似的问题在这里,但我不认为我的问题是重复的,因为我正好在发布前更新都RCPP和RcppArmadillo.其他线程的问题应该由Armadillo修复,所以看来还有其他东西在这里发生.
几点说明:
将expmat()在所提供的算法expm包是从不同的arma::expmat()算法.
据说,从链接的帖子中,这已在4.550.4中修复.这个变化似乎已被包括在内,因为文件源(#3)与armadillo源相同,即使RcppArmadillo似乎已经跳过了点发布,例如4.550.0和4.550.1被合并了.
arma::expmat()命令有了这个说法,让我们用一些调试语句快速浏览幕后.注:我选择来选择T为double每个API文档.
在简短的调试程序上:
# include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
void run_exp_mat_routine(const arma::mat& x) {
const double norm_val = norm(x, "inf");
Rcpp::Rcout << "norm:" << norm_val << std::endl;
Rcpp::Rcout << "Value of T(0):" << (double(0)) << std::endl;
Rcpp::Rcout << "Inequality:" << (norm_val > double(0)) << std::endl;
Rcpp::Rcout << "log2: " << std::log2(norm_val) << std::endl;
int exponent = int(0); std::frexp(std::log2(norm_val), &exponent);
Rcpp::Rcout << "exponent: " << exponent << std::endl;
const arma::uword s = arma::uword( (std::max)(int(0), exponent + int(1)) );
Rcpp::Rcout << "s: " << s << std::endl;
const arma::mat AA = x / std::pow(2.0,s);
Rcpp::Rcout << "AA: " << std::endl << AA << std::endl;
double c = 0.5;
arma::mat E(AA.n_rows, AA.n_rows, arma::fill::eye);
Rcpp::Rcout << "Init E:" << std::endl << E << std::endl;
E += c * AA;
Rcpp::Rcout << "Mod E:" << std::endl << E << std::endl;
arma::mat D(AA.n_rows, AA.n_rows, arma::fill::eye);
Rcpp::Rcout << "Init D:" << std::endl << D << std::endl;
D -= c * AA;
Rcpp::Rcout << "Mod D:" << std::endl << D << std::endl;
arma::mat X = AA;
bool positive = true;
const arma::uword N = 6;
for(arma::uword i = 2; i<=N; ++i){
c = c * double(N - i + 1) / double(i * (2*N - i + 1));
X = AA * X;
E += c * X;
if(positive) { D += c * X; } else { D -= c * X; }
positive = (positive) ? false : true;
Rcpp::Rcout << "Loop: " << i << ", c: " << c << ", positive:" << positive << std::endl;
Rcpp::Rcout << "X: " << std::endl << X << std::endl << "E: " << std::endl << E << std::endl;
}
//arma::mat out = solve(D,E);
// Rcpp::Rcout << "out:" << std::endl << out << std::endl;
//
// for(arma::uword i = 0; i < s; ++i){
// out = out*out;
// }
// Rcpp::Rcout << "out: " << out <<std::endl;
}
/*** R
a <- matrix(c(-2.5654e+060,4.6979e-018,2.5654e+060,7.2765e-035
,2.8913e+000, -3.6633e+001,3.3731e+001,1.0003e-002
,1.0656e-009,1.9037e-002, -1.9732e-001,1.7828e-001
, 0e+000, 0e+000, 0e+000, 0e+000), nrow=4, byrow=T)
run_exp_mat_routine(a)
*/
结果分为初始化阶段,然后是问题所在的循环阶段.
norm: 5.1308e+60
Value of double(0): 0
Inequality: 1
log2: 201.675
exponent: 8
s: 9
AA:
-5.0105e+57 9.1756e-21 5.0105e+57 1.4212e-37
5.6471e-03 -7.1549e-02 6.5881e-02 1.9537e-05
2.0812e-12 3.7182e-05 -3.8539e-04 3.4820e-04
0 0 0 0
Init E:
1.0000 0 0 0
0 1.0000 0 0
0 0 1.0000 0
0 0 0 1.0000
Mod E:
-2.5053e+57 4.5878e-21 2.5053e+57 7.1060e-38
2.8235e-03 9.6423e-01 3.2940e-02 9.7686e-06
1.0406e-12 1.8591e-05 9.9981e-01 1.7410e-04
0 0 0 1.0000e+00
Init D:
1.0000 0 0 0
0 1.0000 0 0
0 0 1.0000 0
0 0 0 1.0000
Mod D:
2.5053e+57 -4.5878e-21 -2.5053e+57 -7.1060e-38
-2.8235e-03 1.0358e+00 -3.2940e-02 -9.7686e-06
-1.0406e-12 -1.8591e-05 1.0002e+00 -1.7410e-04
0 0 0 1.0000e+00
现在,循环部分似乎在最后一次迭代中触发错误(例如i = 6),因为该数字变得太大而无法在双重结构内表示.
Loop: 2, c: 0.113636, positive:0
X:
2.5106e+115 1.8630e+53 -2.5106e+115 1.7447e+54
-2.8295e+55 5.1217e-03 2.8295e+55 2.1542e-05
-1.0428e+46 -2.6746e-06 1.0428e+46 -1.3347e-07
0 0 0 0
E:
2.8529e+114 2.1170e+52 -2.8529e+114 1.9826e+53
-3.2153e+54 9.6481e-01 3.2153e+54 1.2217e-05
-1.1850e+45 1.8287e-05 1.1850e+45 1.7409e-04
0 0 0 1.0000e+00
Loop: 3, c: 0.0151515, positive:1
X:
-1.2579e+173 -9.3347e+110 1.2579e+173 -8.7418e+111
1.4177e+113 1.0521e+51 -1.4177e+113 9.8524e+51
5.2251e+103 3.8774e+41 -5.2251e+103 3.6311e+42
0 0 0 0
E:
-1.9059e+171 -1.4143e+109 1.9059e+171 -1.3245e+110
2.1481e+111 1.5940e+49 -2.1481e+111 1.4928e+50
7.9168e+101 5.8748e+39 -7.9168e+101 5.5017e+40
0 0 0 1.0000e+00
Loop: 4, c: 0.00126263, positive:0
X:
6.3029e+230 4.6772e+168 -6.3029e+230 4.3801e+169
-7.1036e+170 -5.2714e+108 7.1036e+170 -4.9366e+109
-2.6181e+161 -1.9428e+99 2.6181e+161 -1.8194e+100
0 0 0 0
E:
7.9582e+227 5.9055e+165 -7.9582e+227 5.5305e+166
-8.9692e+167 -6.6557e+105 8.9692e+167 -6.2331e+106
-3.3056e+158 -2.4530e+96 3.3056e+158 -2.2972e+97
0 0 0 1.0000e+00
Loop: 5, c: 6.31313e-05, positive:1
X:
-3.1581e+288 -2.3435e+226 3.1581e+288 -2.1947e+227
3.5593e+228 2.6412e+166 -3.5593e+228 2.4735e+167
1.3118e+219 9.7344e+156 -1.3118e+219 9.1162e+157
0 0 0 0
E:
-1.9937e+284 -1.4795e+222 1.9937e+284 -1.3855e+223
2.2470e+224 1.6674e+162 -2.2470e+224 1.5616e+163
8.2815e+214 6.1454e+152 -8.2815e+214 5.7552e+153
0 0 0 1.0000e+00
Loop: 6, c: 1.50313e-06, positive:0
X:
inf 1.1742e+284 -inf 1.0997e+285
-1.7834e+286 -1.3234e+224 1.7834e+286 -1.2394e+225
-6.5728e+276 -4.8775e+214 6.5728e+276 -4.5677e+215
0 0 0 0
E:
inf 1.7650e+278 -inf 1.6529e+279
-2.6807e+280 -1.9892e+218 2.6807e+280 -1.8629e+219
-9.8797e+270 -7.3314e+208 9.8797e+270 -6.8658e+209
0 0 0 1.0000e+00
因此,无穷大符号被传递给solve参数,这将破坏程序.
除了运行一个单独的函数来检查矩阵是否是无限的之外,我不确定还有另一种方法,因为与http://www.cs.cornell.edu/cv/researchpdf/19ways+.pdf相比,算法似乎是合理的.虽然,我会让更多有经验的民间评论这方面.
编辑
R中例程的快速示例:
# install.packages("expm")
library("expm")
a <- matrix(c(-2.5654e+060,4.6979e-018,2.5654e+060,7.2765e-035
,2.8913e+000, -3.6633e+001,3.3731e+001,1.0003e-002
,1.0656e-009,1.9037e-002, -1.9732e-001,1.7828e-001
, 0e+000, 0e+000, 0e+000, 0e+000), nrow=4, byrow=T)
结果:
[,1] [,2] [,3] [,4]
[1,] 2.403680e-62 0.02132743 1.369318 0.1998306
[2,] 1.543272e-60 1.36931834 41.028506 3.4698436
[3,] 2.403680e-62 0.02132743 1.369318 0.1998306
[4,] 0.000000e+00 0.00000000 0.000000 1.0000000
由于此函数应该向MATLAB提供类似的输出(根据链接帖子中的作者),让我们快速运行.
A = [-2.5654e+060,4.6979e-018,2.5654e+060,7.2765e-035;
2.8913e+000, -3.6633e+001,3.3731e+001,1.0003e-002;
1.0656e-009,1.9037e-002, -1.9732e-001,1.7828e-001;
0e+000, 0e+000, 0e+000, 0e+000]
A在MATLAB中的表示是:A =
1.0e+60 *
-2.5654 0.0000 2.5654 0.0000
0.0000 -0.0000 0.0000 0.0000
0.0000 0.0000 -0.0000 0.0000
0 0 0 0
要获得指数矩阵(不是要素的指数),我们使用MATLAB的expm(A):
ans =
0.0000 0.0213 1.3693 0.1998
0.0000 1.3693 41.0285 3.4698
0.0000 0.0213 1.3693 0.1998
0 0 0 1.0000
因此,R和MATLAB版本一致.因此,为犰狳中的矩阵分解选择的实现可能不是理想的.