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