[Help-gsl] Question on Cholesky decomposition error - bspline penalty matrix

[Foivos Diakogiannis]

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