RcppArmadillo expmat 挂起 4x4 矩阵
RcppArmadillo expmat hangs for a 4x4 matrix
我有一个病态的 4x4 矩阵，它使 Armadillo 中的 expmat 函数挂起。病理矩阵为：
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
我知道这个矩阵的条件很差，但是 R 包 "expm" 中的 expm 函数可以使用其默认算法毫无问题地处理它。无论如何在 RcppArmadillo 中解决这个问题？至少我想通过处理警告信息来避免挂起。
有一个类似的问题 here，但我不认为我的问题是重复的，因为我刚刚在发布之前更新了 Rcpp 和 RcppArmadillo。另一个线程中的问题应该已由 Armadillo 修复，所以这里似乎还有其他问题。
快速观察
几点注意事项：
expm package中提供的expmat()算法与arma::expmat()算法不同。
据推测，从链接 post，这已在 4.550.4 中修复。更改似乎已包含在内，因为文件源 (#3) 与 armadillo 源相同，即使 RcppArmadillo 似乎已跳过点发布，例如 4.550.0 and 4.550.1 were merged in。
实际交付的实现（在撰写本文时）是 here and it definitely looks like the unsigned int issue raised on SO 已修复。
正在调试 arma::expmat() 命令
话虽如此，让我们通过一些调试语句快速浏览一下幕后情况。 注意： 根据 API 文档，我选择 select T 作为 double。
进入简短的调试程序：
# 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），因为数字变得太大而无法在 double 结构中表示。
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
因此，无穷大符号被传递给求解参数，这会破坏程序。
除了 运行 一个单独的函数来检查矩阵是否无限之外，我不确定还有另一种方法，因为与 http://www.cs.cornell.edu/cv/researchpdf/19ways+.pdf 相比，该算法似乎很合理。不过，我会让更多有经验的民间人士对此发表评论。
编辑
用 expm::expm(a)
验证 R 中的输出
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
因为这个函数应该提供与 MATLAB 相似的输出（根据作者在链接 post 中的说法）让我们快速 运行.
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中的表示是：
一个=
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 版本一致。因此，在 armadillo 中为矩阵分解选择的实现可能并不理想。
我有一个病态的 4x4 矩阵，它使 Armadillo 中的 expmat 函数挂起。病理矩阵为：
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
我知道这个矩阵的条件很差，但是 R 包 "expm" 中的 expm 函数可以使用其默认算法毫无问题地处理它。无论如何在 RcppArmadillo 中解决这个问题？至少我想通过处理警告信息来避免挂起。
有一个类似的问题 here，但我不认为我的问题是重复的，因为我刚刚在发布之前更新了 Rcpp 和 RcppArmadillo。另一个线程中的问题应该已由 Armadillo 修复，所以这里似乎还有其他问题。
快速观察
几点注意事项：
expmpackage中提供的expmat()算法与arma::expmat()算法不同。据推测，从链接 post，这已在 4.550.4 中修复。更改似乎已包含在内，因为文件源 (#3) 与
armadillo源相同，即使RcppArmadillo似乎已跳过点发布，例如 4.550.0 and 4.550.1 were merged in。实际交付的实现（在撰写本文时）是 here and it definitely looks like the
unsigned intissue raised on SO 已修复。
正在调试 arma::expmat() 命令
话虽如此，让我们通过一些调试语句快速浏览一下幕后情况。 注意： 根据 API 文档，我选择 select T 作为 double。
进入简短的调试程序：
# 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），因为数字变得太大而无法在 double 结构中表示。
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
因此，无穷大符号被传递给求解参数，这会破坏程序。
除了 运行 一个单独的函数来检查矩阵是否无限之外，我不确定还有另一种方法，因为与 http://www.cs.cornell.edu/cv/researchpdf/19ways+.pdf 相比，该算法似乎很合理。不过，我会让更多有经验的民间人士对此发表评论。
编辑
用 expm::expm(a)
验证 R 中的输出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
因为这个函数应该提供与 MATLAB 相似的输出（根据作者在链接 post 中的说法）让我们快速 运行.
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中的表示是：
一个=
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 版本一致。因此，在 armadillo 中为矩阵分解选择的实现可能并不理想。