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

[Rhys Ulerich]

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;
}