Posted by
Dirk Eddelbuettel on
Apr 12, 2017; 1:33am
URL: http://quantlib.414.s1.nabble.com/matrix-inverse-boost-exception-tp18181p18183.html
Peter,
Klaus already nailed it, but here is a little illustration:
--snip---------------------------------------------------------------------------
#include <Rcpp.h>
typedef double Real;
// [[plugins(cpp11)]]
// [[Rcpp::export]]
Rcpp::NumericMatrix getMatrix() {
Real data[8][8] = {{0.0006474874439615692, 0.002486659100900001, 0.0006628267415236041, 0.04673067750857564,
0.1685162462417822, 0.1516341778778243, 0.1219910757564679, 0.507330849328965},
{5.165711089877979e-06, 0.001690175534926361, 2.136045858343317e-05, 0.04475379133340539,
0.1394666326109254, 0.1636645892860261, 0.12811590642367, 0.5222823786413735},
{4.3468717221017e-07, 3.324472127119397e-05, 1.632530899871639e-05, 0.0473437145861729,
0.09087902924968709, 0.1932661669150169, 0.1257053481714577, 0.5427557363602231},
{3.883601299340593e-07, 5.96026473031791e-06, 1.108224823791975e-05, 0.05269312648895734,
0.1170973364539729, 0.2225313741567021, 0.08459583680146567, 0.5230648952258039},
{1.311427219168951e-07, 1.451479950650993e-06, 9.870251521216817e-06, 0.002430611065235592,
0.1217894913776252, 0.2748389234300169, 0.07704266870962403, 0.5238868525433042},
{3.858458712871179e-08, 6.299157109034205e-07, 9.478963012990911e-06, 0.0006947826112014794,
0.007430976289762706, 0.2725490552732179, 0.07467942912312821, 0.6446356092393787},
{8.037084743179829e-09, 1.06623621797325e-07, 7.740443849605495e-06, 0.0005175216683466955,
0.004174347052491365, 0.02170492331427787, 0.06116825771382062, 0.9124270951465073},
{0, 0, 0, 0, 0, 0, 0, 1}};
Rcpp::NumericMatrix m(8, 8);
for (auto i = 0; i < 8; ++i)
for (auto j = 0; j < 8; ++j)
m(i, j) = data[i][j];
return m;
}
/*** R
m <- getMatrix()
kappa(m)
rcond(m)
m
*/
--snip---------------------------------------------------------------------------
If you have R and Rcpp, you can just drop the file somewhere --
/tmp/peter.cpp for me -- and soureCpp() does the rest:
R> sourceCpp("/tmp/peter.cpp")
R> m <- getMatrix()
R> kappa(m)
[1] 2.41674e+23
R> rcond(m)
[1] 2.96131e-22
R> m
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 6.47487e-04 2.48666e-03 6.62827e-04 0.046730678 0.16851625 0.1516342 0.1219911 0.507331
[2,] 5.16571e-06 1.69018e-03 2.13605e-05 0.044753791 0.13946663 0.1636646 0.1281159 0.522282
[3,] 4.34687e-07 3.32447e-05 1.63253e-05 0.047343715 0.09087903 0.1932662 0.1257053 0.542756
[4,] 3.88360e-07 5.96026e-06 1.10822e-05 0.052693126 0.11709734 0.2225314 0.0845958 0.523065
[5,] 1.31143e-07 1.45148e-06 9.87025e-06 0.002430611 0.12178949 0.2748389 0.0770427 0.523887
[6,] 3.85846e-08 6.29916e-07 9.47896e-06 0.000694783 0.00743098 0.2725491 0.0746794 0.644636
[7,] 8.03708e-09 1.06624e-07 7.74044e-06 0.000517522 0.00417435 0.0217049 0.0611683 0.912427
[8,] 0.00000e+00 0.00000e+00 0.00000e+00 0.000000000 0.00000000 0.0000000 0.0000000 1.000000
R> eigen(m)
$values
[1] 1.00000e+00 2.95561e-01 1.08353e-01 5.67581e-02 4.75499e-02 1.69628e-03 6.35697e-04 -1.94874e-18
$vectors
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 0.353553 -0.4648349 0.5813418 -0.5457531 -0.47984632 -9.21607e-01 9.99990e-01 -7.07062e-01
[2,] 0.353553 -0.4280282 0.4971925 -0.5063981 -0.44627394 -3.88072e-01 -4.58103e-03 -4.68906e-15
[3,] 0.353553 -0.3849208 0.3854263 -0.4808008 -0.48972225 -6.41879e-03 -5.34806e-05 7.07151e-01
[4,] 0.353553 -0.4517516 0.4544627 -0.4545852 -0.57391253 4.20525e-05 -4.92962e-06 -2.29916e-18
[5,] 0.353553 -0.4262395 0.2429683 -0.0141278 0.02041553 2.91491e-06 -7.02968e-07 -1.96387e-19
[6,] 0.353553 -0.2564187 -0.0206878 0.0304491 -0.00878029 8.57713e-07 -9.99792e-08 -1.79023e-19
[7,] 0.353553 -0.0323459 0.0170277 -0.0822840 0.02982799 7.72153e-07 8.60803e-09 -8.93925e-05
[8,] 0.353553 0.0000000 0.0000000 0.0000000 0.00000000 0.00000e+00 0.00000e+00 0.00000e+00
R>
The last eigenvalue is below 'epsilon' indicating that something somewhat
shady is going on. Explicitly asking for rank gets me 7:
R> Matrix::rankMatrix(m)
[1] 7
attr(,"method")
[1] "tolNorm2"
attr(,"useGrad")
[1] FALSE
attr(,"tol")
[1] 1.77636e-15
R>
Hth, Dirk
--
http://dirk.eddelbuettel.com | @eddelbuettel |
[hidden email]
------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org!
http://sdm.link/slashdot_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev