Re: [Quantlib-dev] matrix inverse boost exception

[Dirk E. <ed...@de...>]

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 | ed...@de...