http://quantlib.414.s1.nabble.com/matrix-inverse-boost-exception-tp18181p18189.html
Yes, we should improve the error message. The ublas exception message is
On Mittwoch, 12. April 2017 10:26:16 CEST Peter Caspers wrote:
> Hi Klaus, Dirk,
>
> 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,
>
> QL_REQUIRE(singular == 0, "singular matrix given");
>
> 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
>
> // lu decomposition
> try {
> const Size singular = lu_factorize(a, pert);
> } catch (const boost::numeric::ublas::internal_logic& e) {
> QL_FAIL("lu_factorize error: " << e.what());
> } catch (const boost::numeric::ublas::external_logic& e) {
> QL_FAIL("lu_factorize error: " << e.what());
> }
> QL_REQUIRE(singular == 0, "singular matrix given");
>
> Matrix retVal(m.rows(), m.columns());
> boost::numeric::ublas::matrix<Real>
> inverse =
> boost::numeric::ublas::identity_matrix<Real>(m.rows());
>
> // backsubstitution
> try {
> boost::numeric::ublas::lu_substitute(a, pert, inverse);
> } catch (const boost::numeric::ublas::internal_logic& e) {
> QL_FAIL("lu_substitute error: " << e.what());
> }
>
> Is it worth adding this?
>
> Best Regards
> Peter
>
> > 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
>
> ----------------------------------------------------------------------------
> -- 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