Hi Foivos, > I guess it must be a matter of numerical accuracy...
I'm not convinced: B-splines are wonderfully stable. Gaussian integration is wonderfully scare on FLOPS during which accumulation could creep in. Cholesky factorization is well-behaved. > 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... Supposing you copy the following C99 source code into penalty.c, compile/link it in the usual fashion into a.out, and run 'GSL_RNG_SEED=$(date +%s) ./a.out' you'll see realizations like the following: GSL_RNG_SEED=1416368157 nbreak: 8 korder: 4 left: 0 right: 1 ndof: 10 nderiv: 2 Breakpoints 0 0.109812 0.238165 0.398182 0.448533 0.605244 0.674791 1 Penalty: 1.3460992e+09 2.6746594e+09 2.066323e+08 3627533.6 0 0 0 0 0 0 2.6746594e+09 5.3587223e+09 4.2768015e+08 6082748.9 180632.43 0 0 0 0 0 2.066323e+08 4.2768015e+08 58933490 4140286.1 574969.35 87924.828 0 0 0 0 3627533.6 6082748.9 4140286.1 7514918.3 6168763.3 1184369.7 130591.26 0 0 0 0 180632.43 574969.35 6168763.3 22428201 10986747 1810423.9 75268.263 0 0 0 0 87924.828 1184369.7 10986747 34866380 13873749 621075.37 31586.637 0 0 0 0 130591.26 1810423.9 13873749 15645715 1629739.1 108965.28 26503.364 0 0 0 0 75268.263 621075.37 1629739.1 1206476.1 1009717.2 286360.17 0 0 0 0 0 31586.637 108965.28 1009717.2 5502316.4 1875860.5 0 0 0 0 0 0 26503.364 286360.17 1875860.5 673752.58 Eigenvalues: 6.72975e+09 4.71519e+07 2.80286e+07 2.04277e+07 3.77445e+06 8.24125e+06 7.10832e+06 6.29724e+06 789209 26786.9 I can run this repeatedly without the included Cholesky factorization complaining, suggesting it is at least producing positive definite matrices. Before you use any of this logic, check it against a few B-spline orders and basis function counts where you know the exact solution. Note the TODOs. I do not claim this code is more than structurally correct. I've not gotten to play with numerics, especially not my contributions to the GSL, in the few months since I graduated. Thanks for the divertissement. - Rhys #include <stdio.h> #include <gsl/gsl_bspline.h> #include <gsl/gsl_eigen.h> #include <gsl/gsl_integration.h> #include <gsl/gsl_linalg.h> #include <gsl/gsl_math.h> #include <gsl/gsl_matrix.h> #include <gsl/gsl_rng.h> #include <gsl/gsl_vector.h> int main(const int argc, const char *argv[]) { // Establish RNG from environment gsl_rng_env_setup(); gsl_rng * const r = gsl_rng_alloc(gsl_rng_default); // Parse command line arguments // I got into an fprintf kick for reasons I don't quite understand const size_t nbreak = argc > 1 ? atoi(argv[1]) : 8; const size_t korder = argc > 2 ? atoi(argv[2]) : 4; const double left = argc > 3 ? atof(argv[3]) : 0.0; const double right = argc > 4 ? atof(argv[4]) : 1.0; const size_t nderiv = argc > 5 ? atof(argv[5]) : 2; const size_t ndof = nbreak + korder - 2; fprintf(stdout, "nbreak:\t%zd\n", nbreak); fprintf(stdout, "korder:\t%zd\n", korder); fprintf(stdout, "left:\t%g\n", left); fprintf(stdout, "right:\t%g\n", right); fprintf(stdout, "ndof:\t%zd\n", ndof); fprintf(stdout, "nderiv:\t%zd\n", nderiv); // Prepare basic storage gsl_bspline_workspace * const bw = gsl_bspline_alloc(korder, nbreak); gsl_bspline_deriv_workspace * const dbw = gsl_bspline_deriv_alloc(korder); gsl_vector * const breakpts = gsl_vector_alloc(nbreak); gsl_matrix * const penalty = gsl_matrix_calloc(ndof, ndof); // Obtain a randomly-spaced set of breakpoints gsl_vector_set(breakpts, 0, left); for (size_t i = 1; i < nbreak - 1; ++i) { const double p = gsl_vector_get(breakpts, i - 1); gsl_vector_set(breakpts, i, p + gsl_rng_uniform(r)/3 * (right - p)); } gsl_vector_set(breakpts, nbreak - 1, right); gsl_bspline_knots(breakpts, bw); fprintf(stdout, "\nBreakpoints\n"); gsl_vector_fprintf(stdout, breakpts, "\t%g"); // Compute penalty d^nderiv B_i * d^nderiv B_j via Gauss-Legendre rule // TODO Check indexing is stride one before using this important places // TODO Confirm that the order of integration is correct per korder, nderiv // TODO Could presumably update only one triangular portion of penalty gsl_integration_glfixed_table * const tbl = gsl_integration_glfixed_table_alloc(4 * (korder - nderiv + 1)/2); gsl_matrix * const dB = gsl_matrix_alloc(korder, nderiv + 1); for (size_t k = 0; k < (nbreak - 1); ++k) { for (size_t l = 0; l < tbl->n; ++l) { // Lookup the l-th Gauss point and weight in breakpoint subinterval double xl, wl; gsl_integration_glfixed_point(gsl_bspline_breakpoint(k, bw), gsl_bspline_breakpoint(k + 1, bw), l, &xl, &wl, tbl); // Evaluate basis functions at point xl size_t start, end; gsl_bspline_deriv_eval_nonzero(xl, nderiv, dB, &start, &end, bw, dbw); // Compute (i, j)-th penalty contribution using known basis support for (size_t i = start; i <= end; ++i) { const double dB_i = gsl_matrix_get(dB, i - start, nderiv); for (size_t j = start; j <= end; ++j) { const double dB_j = gsl_matrix_get(dB, j - start, nderiv); double * const p = gsl_matrix_ptr(penalty, i, j); *p += wl * gsl_pow_2(dB_i * dB_j); } } } } gsl_matrix_free(dB); gsl_integration_glfixed_table_free(tbl); // Display the penalty matrix fprintf(stdout, "\nPenalty:\n"); for (size_t i = 0; i < ndof; ++i) { for (size_t j = 0; j < ndof; ++j) { fprintf(stdout, "%16.8g", gsl_matrix_get(penalty, i, j)); } fputc('\n', stdout); } // Display the eigenvalues gsl_eigen_symm_workspace * const eigs = gsl_eigen_symm_alloc(ndof); { gsl_matrix * const A = gsl_matrix_alloc(ndof, ndof); gsl_vector * const eval = gsl_vector_alloc(ndof); gsl_matrix_memcpy(A, penalty); gsl_eigen_symm (A, eval, eigs); fprintf(stdout, "\nEigenvalues:\n"); gsl_vector_fprintf(stdout, eval, "\t%g"); gsl_vector_free(eval); gsl_matrix_free(A); gsl_eigen_symm_free(eigs); } // Compute the Cholesky decomposition, which should not fail gsl_matrix * const cholesky = gsl_matrix_alloc(ndof, ndof); gsl_matrix_memcpy(cholesky, penalty); gsl_linalg_cholesky_decomp(cholesky); // Clean up after ourselves gsl_matrix_free(cholesky); gsl_matrix_free(penalty); gsl_vector_free(breakpts); gsl_bspline_deriv_free(dbw); gsl_bspline_free(bw); gsl_rng_free(r); return 0; }