Hi Rhys and all, Thank you very much for the links, the code and your useful comments.
------------------------------------------------ When you say it works sometimes, do you mean it works for "nicely placed" B-spline breakpoints and then fails when they're bunched up in any fashion? -------------------------------- No it breaks down even for "nice" (relatively evenly spaced) breakpoints. I've also tested the numerical integration result with Mathematica10.0 it gives the same results. I would expect that it may break down for cases when the breakpoints are too close to each other, but this is not always the case plus in the following I use a method independent of the breakpoints. I tried the Gauss Legendre integration from GSL, which is exact, and I must be mixing something pretty badly since I can't make it to work (I have to give n >> order of polynomial, still changing n, changes slightly the result, differs from GSL_CQUAD). Nevertheless I tried a different approach which is exact (machine precision exact, I guess...), and tested with GSL_CQUAD as well --> gives same results. The method I followed is based on Piegl & Tiller 1997, Symbolic operators for NURBS. Essentially they describe an algorithm on how to create a new bspline function from the product of two bspline functions. So I applied this for each pair of derivatives of b-spline basis. That is Construct f(x) = d^n B_i(x)/dx d^n B_j(x)/dx = ... = \sum_r a^r Bnew_r(x) for some new a^r coefficients and new b-spline basis Bnew_r(x). This is a new b-spline function and thus has an exact formula for the evaluation of its integral \int_0^1 f(x) dx (e.g. de Boor ). I get the same results as CQUAD, yet again the Cholesky decomposition breaks down. This is actually dependent on how many digits of accuracy I will keep. I've tested this against an online tool that performs Cholesky decomposition and matrix operations, here: http://www.bluebit.gr/matrix-calculator/ With some digits of accuracy, goes well, if I increase them it breaks down. Another peculiar thing is that sometimes the GSL implementation gives me one negative (but very small) eigenvalue, but it does not exit with error, so I assume there is some tolerance for the calculations. Here is the result of some of the output calculations (GSL) where I have the matrix and the eigenvalues. // ***************************************************************************************** // Run 1: *success* Breakpoints vector: 0 0.22823 0.50927 0.928461 1 The matrix we are going to decompose: 10.015 -7.20715 -2.51511 -0.28011 -0.01259 0 0 0 -7.20715 7.97199 0.349051 -0.895536 -0.214722 -0.00363725 0 0 -2.51511 0.349051 2.73437 0.359693 -0.654201 -0.253061 -0.0207373 0 -0.28011 -0.895536 0.359693 1.38283 0.352124 -0.625593 -0.291864 -0.00154437 -0.01259 -0.214722 -0.654201 0.352124 1.54113 0.367551 -1.25646 -0.122834 0 -0.00363725 -0.253061 -0.625593 0.367551 3.14948 0.535466 -3.17021 0 0 -0.0207373 -0.291864 -1.25646 0.535466 29.6895 -28.6559 0 0 0 -0.00154437 -0.122834 -3.17021 -28.6559 31.9505 Eigen values of Omega_ij: < 59.6342, 16.6121, 4.82787, 3.73401, 2.2145, *-9.59056e-17*, 0.779547, 0.632468, > // *********************************************************************************************** observe that despite the small negative value for the eigenvalue, it exits with success. // ********************************************************************************************** // Run 2: *error * Breakpoints vector: 0 0.101753 0.134297 0.321611 1 The matrix we are going to decompose: 22.4633 -13.0785 -8.05851 -1.29888 -0.0273963 0 0 0 -13.0785 15.2429 0.613246 -2.63086 -0.146771 -9.70897e-06 0 0 -8.05851 0.613246 9.36671 -0.782056 -1.07152 -0.0660831 -0.00178624 0 -1.29888 -2.63086 -0.782056 5.52198 0.27101 -0.679176 -0.334383 -0.067637 -0.0273963 -0.146771 -1.07152 0.27101 1.47206 0.309883 -0.430117 -0.377149 0 -9.70897e-06 -0.0660831 -0.679176 0.309883 0.992239 0.438947 -0.995801 0 0 -0.00178624 -0.334383 -0.430117 0.438947 2.25608 -1.92874 0 0 0 -0.067637 -0.377149 -0.995801 -1.92874 3.36933 Eigen values of Omega_ij: < 34.3333, 11.1503, 6.83063, 5.07682, 1.93399, *-7.40293e-16,* 0.819623, 0.539926, > gsl: cholesky.c:157: ERROR: matrix must be positive definite Default GSL error handler invoked. Aborted (core dumped) // **************************************************************************************************************** Thanks for all the help. I guess it must be a matter of numerical accuracy since according to my knowledge, the matrix Omega_{i,j} = \int_rmin^rmax d^n B_i(x) d^n B_j(x) dx is positive definite by construction (e.g. Wood, Generalized Additive models, an introduction with R). All the best, Foivos On Mon, Nov 17, 2014 at 8:58 PM, Rhys Ulerich <[email protected]> wrote: > 'Morning, > > > 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. > > It may very well be round off in the construction, but Cholesky > doesn't care. That any eigenvalue is coming back negative means you've > fed Cholesky bad input. > > When you say it works sometimes, do you mean it works for "nicely > placed" B-spline breakpoints and then fails when they're bunched up in > any fashion? > > > 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. > > Adaptive integration crossing breakpoints tends to be, frankly, lousy > relative to what one can do with knowledge of the B-spline basis and > the high-order Gauss rules in the GSL. Choose an exact Gauss rule > based on your B-spline piecewise polynomial order and what it implies > about your integrand. Iterate across the breakpoints and apply the > Gauss rule on each subinterval. Sum the subinterval results to get > the full integral. Nothing adaptive required-- you've just done the > integral to O(eps). > > An approach similar to > https://github.com/RhysU/suzerain/blob/master/suzerain/bspline.h#L240 > https://github.com/RhysU/suzerain/blob/master/suzerain/bspline.c#L228 > which uses the helper function > https://github.com/RhysU/suzerain/blob/master/suzerain/omsect.h#L39 > is what you want. That's not precisely your use case, so be sure to > read what it's doing. In particular, I don't think you'll need any of > the omsect(...) goofiness. > > If you want to avoid messing with Gauss rules, you could try the same > breakpoint-by-breakpoint approach with GSL_CQUAD. GIven sufficiently > many points it should agree with the Gauss point computation. > > Hope that helps, > Rhys >
