QuantLib › quantlib-dev

matrix inverse boost exception

Classic

List

Threaded

5 messages Options

Peter Caspers-4

matrix inverse boost exception

Hi,

the following Code

#include <ql/math/matrix.hpp>
#include <iostream>
using namespace QuantLib;
int main() {
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}};
Matrix m(8, 8);
for (Size i = 0; i < 8; ++i)
for (Size j = 0; j < 8; ++j)
m(i, j) = data[i][j];
std::clog << inverse(m) << std::endl;
return 0;
}

results in

Check failed in file /usr/local/include/boost/numeric/ublas/lu.hpp at line 298:
detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, upper> (m), e), cm2)
libc++abi.dylib: terminating with uncaught exception of type boost::numeric::ublas::internal_logic: internal logic

for me. Does anyone have an idea what might go wrong here?

I am using Boost 1.59 + a recent QL version on OSX with clang 4.0.0. I believe QuantLib only wraps ublas functions here, so probably I should create a ticket in the boost project?

Thanks a lot
Peter
------------------------------------------------------------------------------
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

Klaus Spanderen-2

Re: matrix inverse boost exception

Hi Peter,

IMO the condition number of your matrix is very large (matrix is singular or computationally singular) and the ublas check for numerical stability is raising an alert here.

http://stackoverflow.com/questions/17240624/matrix-inversion-in-boost

BOOST_UBLAS_TYPE_CHECK_EPSILON default: sqrt(epsilon), controls how large the difference between the expected result and the computed result may become. Increase this value if you are going to use near singular or badly scaled matrices.

Tuning this epsilon doesn't help. If I'm making epsilon big enough so that the exception disappears then the outcome is flat wrong due to the very large condition number.

regards

Klaus

On Dienstag, 11. April 2017 20:59:11 CEST Peter Caspers wrote:

> Hi,

> the following Code

> #include <ql/math/matrix.hpp>

> #include <iostream>

> using namespace QuantLib;

> int main() {

> 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}};

> Matrix m(8, 8);

> for (Size i = 0; i < 8; ++i)

> for (Size j = 0; j < 8; ++j)

> m(i, j) = data[i][j];

> std::clog << inverse(m) << std::endl;

> return 0;

> }

> results in

> Check failed in file /usr/local/include/boost/numeric/ublas/lu.hpp at line

> 298: detail::expression_type_check (prod

> (triangular_adaptor<const_matrix_type, upper> (m), e), cm2)

> libc++abi.dylib: terminating with uncaught exception of type

> boost::numeric::ublas::internal_logic: internal logic

> for me. Does anyone have an idea what might go wrong here?

> I am using Boost 1.59 + a recent QL version on OSX with clang 4.0.0. I

> believe QuantLib only wraps ublas functions here, so probably I should

> create a ticket in the boost project?

> Thanks a lot

> Peter

> ----------------------------------------------------------------------------

> -- 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

------------------------------------------------------------------------------
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

Dirk Eddelbuettel

Re: matrix inverse boost exception

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

Peter Caspers-4

Re: matrix inverse boost exception

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
>
> --
> http://dirk.eddelbuettel.com | @eddelbuettel | [hidden email]

Klaus Spanderen-2

Re: matrix inverse boost exception

Hi Peter,

Yes, we should improve the error message. The ublas exception message is
really misleading.

regards
Klaus

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