http://quantlib.414.s1.nabble.com/matrix-inverse-boost-exception-tp18181p18186.html
thanks a lot, that makes sense. I somehow expected to get a ql error if the input matrix is (numerically) singular: The return value of lu_factorize seems to be zero for an invertible matrix and non-zero for a singluar matrix and this value is already checked in the ql code,
But in the case below singular is zero (so the lu_factorize step still seems to works ok) and the exception is thrown later in the lu_substitute step. I wonder if we should add some error handling to the ql code to have a consistent “answer” from the inverse function, e.g. something like that
inverse = boost::numeric::ublas::identity_matrix<Real>(m.rows());
> On 12 Apr 2017, at 03:33, Dirk Eddelbuettel <
[hidden email]> wrote:
>
>
> 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]