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

[Foivos Diakogiannis]

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 <address@hidden>
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
>