Curve fit matrices are “almost always” positive definite

“A real matrix A is positive definite if and only it is symmetric, that is A^T=A, and the quadratic form x^TAx is positive for all nonzero vectors x.”[1]

                Consider

                 (1)

Which is greater than zero for all values of S_j(x_i).  Note that this includes the case

                            (2)

Expanding the squared term out and interchanging the order of the sums

                (3)

The sum of the S’s in 3 with the definitions in (2) appears in FitData.doc .htm in equation 9 as an approximation to the second derivative expansion of the c² to be minimized in curve fitting.   This approximation is very good, and the second derivative correction terms usually do not keep the matrix from being positive definite, but there are pathological cases.  The nlfit codes (nlfit\Welcome.doc .htm) always use the approximation of the first term given in equation 9.  Thus making the scond derivative matrix is always appropriate for Cholesky matrix inversion.

In the event that one goes to the effort of putting in the second term of equation 8, one should also use GaussJ.doc htm to invert the matrix.

               

So that PPCC is seen to be a positive definite matrix.  Then we note that adding to this ål_ic_ic_i leaves it positive definite for all l_I > 0.  Thus when symmetric the matrices that we wish to invert are Positive definite matrices.

                The first relevant fact is that in the decomposition

A=GG^T which precedes the full matrix inversion we end up with

|g_i,j|£Öa_i,i. “In other words, the elements of G are bounded, even when pivoting is not used.”[2]  Begin with two very nice features

1. Only half the storage is required n(n+1)/2  versus n²

2. Only half the operations are required 

Then add to this the fact that when the system fails due to truncation error it gives a very easily interpreted sign (the need to take the square root of a negative number).

                Cholesky.for  -- the comments are ready for Bailey’s multiple precsion ..\..\..\MultiplePrecision\bailey\Welcome.htm

The most tested version of Cholesky minimization is the subroutine sminv in Fittery\nlfit\formp\robmin.for.  Multiple precision versions of this mphsminv and mpsminv also appear.  These do not use Bailey’s version of multiple precision, but rather use the very similar one discussed in MultiplePrecision\Bob\Welcome.htm. 

 

      JJ=0

      DO J=1,N

        S=0

        KJ=JJ

        KK=0

        IF(J.GE.2)THEN

          JM1=J-1

          DO K=1,JM1

            KJ=KJ+1

            T=AP(KJ)

            KM1=K-1

            DO I=1,KM1

              T=T-AP(I+KK)*AP(I+JJ)

            ENDDO

            KK=KK+K

            T=T/AP(KK)

            AP(KJ)=T

            IF(T.GT.1.D75.OR.T.LT.-1.D75)THEN

              PRINT*,' AP(KK)=',KK,AP(KK)

              T=DMAX1(-1.D75,DMIN1(1.D75,T))

            ENDIF

            S=S+T*T

          ENDDO

        ENDIF

        JJ=JJ+J

        PRINT*,' JJ=',JJ

        SMT=DMAX1(S,AP(JJ))

        DIFF=AP(JJ)-S

        IF(DIFF.LT.1.D-10*SMT)THEN

C         PRINT*, 'POS DEF < 0, J, JJ AP(JJ),S',J,JJ,AP(JJ),S

C *** FIX IS MORE SMOOTHING FOR A(J,J)

          IFL=-1

          DIFF=1.D30

        ENDIF

        AP(JJ)=DSQRT(DIFF)

      ENDDO

      do i=1,6

        print*,ap(i)

      enddo

      STOP

C *** NEXT WE CONSTRUCT THE INVERSE MATRIX

CHOLESKY.CPP

void main(void)

    {int i,j,k,n=3,ix=57,iy=89,ifl,itest;

    double a[3][3],at[3][3],ai[3][3],ap[6];

    rseed(ix,iy);

    for (i=0;i<3;++i){

      for (j=i;j<3;++j){

        a[i][j]=rndmf()*.1;;

        a[j][i]=a[i][j];}}

 

    for (j=0;j<3;++j){

      a[j][j]=1;}

    k=0;

    for(i=0;i<3;++i){

      for (j=0;j<=i;++j){

        ap[k]=a[j][i];

        k+=1;}}

    ifl=cholesky(ap,n);

    printf(" ifl from cholesky is %d\n",ifl);

 

 

    k=0;

    for (i=0;i<3;++i){

      for (j=0;j<=i;++j){

        ai[i][j]=ap[k];

        ai[j][i]=ai[i][j];

        k+=1;}}

 

    for(i=0;i<3;i++){

      for (j=0;j<3;j++){

        at[j][i]=0;

        for(k=0;k<3;k++){

          at[j][i]+=ai[i][k]*a[k][j];}}}

 

    printf(" orig matrix \n");

    for (i=0;i<3;i++)

      {for (j=0;j<3;j++)printf(" %d %d %6f ",j,i,a[j][i]);

       printf("\n");}

    printf(" \n\n");

    printf(" unit matrix? \n");

    for (i=0;i<3;i++)

      {for (j=0;j<3;j++)printf(" %d %d %6f ",j,i,at[j][i]);

       printf("\n");}

    scanf("%d",&itest);

    return; }