[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Help-gsl] Question on Cholesky decomposition error - bspline penalty ma
From: |
Foivos Diakogiannis |
Subject: |
[Help-gsl] Question on Cholesky decomposition error - bspline penalty matrix |
Date: |
Mon, 17 Nov 2014 14:06:16 +0800 |
Dear all,
I am getting an error in the Cholesky decomposition of a matrix which by
construction is positive definite, however I do not get this error always,
so I suspect that may be a round off error or something similar.
This is what I do, I define a set of bspline basis (well tested) with some
random breakpoints in the region [0, ... , 1.], then I construct a bspline
penalty matrix Omega_ij, according to the definition
\Omega_{ij} = \int_0^1 { d^2B_i(x)/dx^2 * d^2 B_j(x)/dx^2 } dx
The accuracy of this integration is 1.e-10 (routine GSL_CQUAD), this matrix
is by definition positive definite, so it must have a Cholesky
decomposition. Most frequently though the code exits with: gsl:
cholesky.c:157: ERROR: matrix must be positive definite. When I calculate
the eigenvalues of the matrix, any negative one is smaller than the
accuracy of the integration.
Is the problem, a problem of numerical accuracy of the integration or am I
missing something? If it is, any ideas of how can I overcome this problem?
This is an example of a non successful run:
//
********************************************************************************************
Breakpoints vector: 0 0.0787706 0.255106 0.361625 1
The matrix we are going to decompose:
58925.1 -73663.5 12444.1 2242.28 52.0309
0 0 0
-73663.5 93441.6 -17415.7 -2453.75 84.6702
6.6039 0 0
12444.1 -17415.7 5553.91 -478.595 -131.754 27.7066
0.37945 0
2242.28 -2453.75 -478.595 802.183 -79.4845 -47.1953
7.56762 6.99486
52.0309 84.6702 -131.754 -79.4845 108.321 -13.2056
-40.4281 19.8509
0 6.6039 27.7066 -47.1953 -13.2056 78.5354
-77.8817 25.4367
0 0 0.37945 7.56762 -40.4281
-77.8817 273.35 -162.987
0 0 0 6.99486 19.8509
25.4367 -162.987 110.704
Eigen values of Omega_ij: < 154969, 3316.71, 492.1, 391.685, 90.9888,
33.4476, 1.71923e-12, *-6.31159e-13* >
gsl: cholesky.c:157: ERROR: matrix must be positive definite
Default GSL error handler invoked.
Aborted (core dumped)
//
********************************************************************************************
The negative eigenvalue is below the accuracy of the numerical integration.
Thank you very much for your time.
Foivos
- [Help-gsl] Question on Cholesky decomposition error - bspline penalty matrix,
Foivos Diakogiannis <=