Appendix D: The square root of the inverse of the covariance matrix

There are various ways in which the square root of the inverse of the covariance matrix can be obtained. A numerically stable and simple method is by taking the square root of the covariance matrix first, and then inverting this lower triangular matrix. Alternatively, one could invert the covariance matrix, and then take the square root (but this could involve less numerically stable calculations).

The algorithm for the Cholesky factorization of the symmetric matrix CV is as follows:

    M = N*(N+1)/2
    DO I = 1, M
      LT(I) = CV(I)
      END DO
    DO I = N, 2, -1
      LT(M) = SQRT(LT(M))
      K = 0
      M = M - I
      L = M
      DO J1 = 1, I-1
        L = L+1
        LT(L) = LT(L)/LT(M)
        DO J2 = 1, J1
          K = K+1
          LT(K) = LT(K) - LT(M+J2)*LT(L)
          END DO
        END DO
      END DO
    LT(1) = SQRT(LT(1))

from which is obtained the lower (or upper) triangular matrix LT. The determinant of the covariance matrix CV, which provides a measure of the total amount of information available in the observations, can be obtained from the product of the squared diagonal elements of LT. The matrix LT is inverted using:

    K = 0
C
C  Creating unit matrix
C
    DO I = 1, N
      K= K + 1
      LI(K) = 1.0
      DO J = 2, I
        K = K + 1
        LI(K) = 0.0
        END DO
      END DO
C
C   Starting inversion process
C
    M = N*(N+1)/2
    DO I = N, 1, -1
      L = M
      DO J1 = I, 1, -1
        LI(M) = LI(M) / LT(L)
        K1 = M
        DO J2 = J1-1, 1, -1
          L = L-1
          K1 = K1-1
          LI(K1) = LI(K1) - LI(M)*LT(L)
          END DO
        M = M-1
        L = L-1
        END DO
      END DO

giving the square root LI of the inverse of the covariance matrix CV.