Chapter 20: Mathematical Support

MCHOL

Computes an upper triangular factorization of a real symmetric matrix A plus a diagonal matrix D, where D is determined sequentially during the Cholesky factorization in order to make A + D nonnegative definite.

Required Arguments

A N by N symmetric matrix for which a Cholesky factorization is attempted.   (Input)
Only elements in the upper triangle and diagonal of A are referenced.

IRANK Rank of A + D.   (Output)

R N by N upper triangular matrix containing the R matrix from a Cholesky decomposition RTR of A + D.   (Output)
The lower triangle of R is not referenced. If A is not needed, then R and A can share the same storage locations.

DMAX Largest diagonal element of D.   (Output)
If DMAX equals 0.0, then A is nonnegative definite, and R is a Cholesky factorization of A. If DMAX is positive, then A is indefinite, i.e., A has at least one negative eigenvalue. In this case, DMAX is an upper bound on the absolute value of the minimum eigenvalue.

IND Index for subsequent computation of a descent direction in the case of a saddle point.   (Output)
If IND = 0, then A is nonnegative definite. For positive IND, let e be a unit vector with IND-th element 1 and remaining elements 0. The solution s of Rs = e is a direction of negative curvature, i.e., sT As is negative.

Optional Arguments

N Order of the matrix.   (Input)
Default: N = size (A,2).

LDA Leading dimension of A exactly as specified in the dimension statement in the calling program.   (Input)
Default: LDA = size (A,1).

TOL Tolerance used in determining linear dependence.   (Input)
TOL = 100 * AMACH(4) is a common choice. See documentation for routine AMACH.
Default: TOL = 1.e-5 for single precision and 2.d –14 for double precision.

LDR Leading dimension of R exactly as specified in the dimension statement in the calling program.   (Input)
Default: LDR = size (R,1).

FORTRAN 90 Interface

Generic:         CALL MCHOL (A, IRANK, R, DMAX, IND [,…])

Specific:                             The specific interface names are S_MCHOL and D_MCHOL.

FORTRAN 77 Interface

Single:                                CALL MCHOL (N, A, LDA, TOL, IRANK, R, LDR, DMAX, IND)

Double:                              The double precision name is DMCHOL.

Description

Routine MCHOL computes a Cholesky factorization, RTR, of A + D where A is symmetric and D is a diagonal matrix with sufficiently large diagonal elements such that A + D is nonnegative definite. The routine is similar to one described by Gill, Murray, and Wright (1981, pages 108111). Here, though, we allow A + D to be singular.

The algorithm proceeds sequentially by rows. If A + D is singular, the Cholesky factor R is taken to have some rows that are entirely zero. The i-th row of A + D is declared to be linearly dependent on the first i 1 rows if the following two conditions are satisfied:

where ε is the input argument TOL.

The routine MCHOL is often used to find a descent direction in a minimization problem. Let A and g be the current Hessian and gradient, respectively, associated with the minimization problem. The solution s of As = g may not give a descent direction if A is not nonnegative definite. Instead, in order to guarantee a descent direction, a solution s of (A + D)s = g can be found where A + D is nonnegative definite. Routine MCHOL is useful for computing the upper triangular Cholesky factor R of A + D so that routine GIRTS can be invoked to compute the descent direction s by solving successively the two triangular linear systems RTx = g and Rs = x for x and then s. Also if g = 0 and A is not nonnegative definite, i.e., the current solution is a saddle point, GIRTS can be used to compute a descent direction s from the linear system Rs = e where e is a unit vector with

Example 1

A Cholesky factorization of a 5 × 5 symmetric nonnegative definite matrix is computed. Maindonald (1984, pages 8586) discusses the example.

 

!                                 SPECIFICATIONS FOR PARAMETERS

      USE MCHOL_INT

      USE UMACH_INT

      USE WRRRN_INT

 

      IMPLICIT   NONE

 

      INTEGER    LDA, LDR, N, J

      PARAMETER  (N=5, LDA=N, LDR=N)

!

      INTEGER    IND, IRANK, NOUT

      REAL       A(LDA,N), DMAX, R(LDR,N), TOL

!

      DATA (A(1,J),J=1,N)/36.0, 12.0, 30.0, 6.0, 18.0/

      DATA (A(2,J),J=1,N)/12.0, 20.0, 2.0, 10.0, 22.0/

      DATA (A(3,J),J=1,N)/30.0, 2.0, 29.0, 1.0, 7.0/

      DATA (A(4,J),J=1,N)/6.0, 10.0, 1.0, 14.0, 20.0/

      DATA (A(5,J),J=1,N)/8.0, 22.0, 7.0, 20.0, 40.0/

!

      TOL = 0.00001

      CALL MCHOL (A, IRANK, R, DMAX, IND, TOL=TOL)

      CALL UMACH (2, NOUT)

      WRITE (NOUT,99998) ' IRANK = ', IRANK

      WRITE (NOUT,99999) ' DMAX =  ', DMAX

      WRITE (NOUT,99998) ' IND =   ', IND

99998 FORMAT (A, I3)

99999 FORMAT (A, 1PE10.3)

      CALL WRRRN ('R', R)

      END

Output

 

IRANK =   4

DMAX =   0.000E+00

IND =     0

                    R

        1       2       3       4       5

1   6.000   2.000   5.000   1.000   3.000

2   0.000   4.000  -2.000   2.000   4.000

3   0.000   0.000   0.000   0.000   0.000

4   0.000   0.000   0.000   3.000   3.000

5   0.000   0.000   0.000   0.000   2.449

Additional Example

Example 2

A modified Cholesky factorization of a 3 × 3 symmetric indefinite matrix A is computed. A solution of Rs = e3 is also obtained using routine GIRTS. Note that sT As  is negative as veri­fied by using routine BLINF (IMSL MATH/LIBRARY). Gill, Murray, and Wright (1981, page 111) discuss the example.

 

!                                 SPECIFICATIONS FOR PARAMETERS

      USE IMSL_LIBRARIES

 

      IMPLICIT   NONE

 

      INTEGER    LDA, LDR, N, J

      PARAMETER  (N=3, LDA=N, LDR=N)

!

      INTEGER    IND, IRANK, NOUT

      REAL       A(LDA,N), DMAX, E(N,1), R(LDR,N), S(N, 1), SPAS, TOL

!

      DATA (A(1,J),J=1,N)/1, 1, 2/

      DATA (A(2,J),J=1,N)/1, 1, 3/

      DATA (A(3,J),J=1,N)/2, 3, 1/

!

      TOL = 0.00001

      CALL MCHOL (A, IRANK, R, DMAX, IND, TOL=TOL)

      CALL UMACH (2, NOUT)

      WRITE (NOUT,99998) ' IRANK = ', IRANK

      WRITE (NOUT,99999) ' DMAX  = ', DMAX

      WRITE (NOUT,99998) ' IND   = ', IND

      CALL WRRRN ('R', R)

      IF (IND .GT. 0) THEN

         E= 0.0

         E(IND, 1) = 1.0

         CALL GIRTS (R, E, IRANK, S)

         SPAS = BLINF(A, S(:,1), S(:,1))

         WRITE (NOUT,*) ' '

         WRITE (NOUT,99999) ' trans(s)*A*s = ', SPAS

      END IF

99998 FORMAT (A, I3)

99999 FORMAT (A, F10.3)

      END

Output

 

IRANK =   3

DMAX  =      5.016

IND   =   3


            R

        1       2       3

1   1.942   0.515   1.030

2   0.000   2.398   1.030

3   0.000   0.000   1.059


trans(s)*A*s =     -2.254



http://www.vni.com/
PHONE: 713.784.3131
FAX:713.781.9260