[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Help-gsl] Question on Cholesky decomposition error - bspline penalt
From: |
Foivos Diakogiannis |
Subject: |
Re: [Help-gsl] Question on Cholesky decomposition error - bspline penalty matrix |
Date: |
Tue, 18 Nov 2014 16:42:02 +0800 |
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
>