R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

>  # tryCchol.r
>  #
>  # demonstration of tryCatch with Cholesky factorization
>  #
>  # three test matrices
>  # first is 3x3, positive definite
>   A1 <- matrix(c(9,2,-2,2,1,0,-2,0,4),3,3)
>   A1
     [,1] [,2] [,3]
[1,]    9    2   -2
[2,]    2    1    0
[3,]   -2    0    4
>  # second is 4x4, cov mx of multinomial, positive semi-definite
>   p <- c(.1,.2,.3,.4)
>   A2 <- diag(p) - outer(p,p)
>   A2
      [,1]  [,2]  [,3]  [,4]
[1,]  0.09 -0.02 -0.03 -0.04
[2,] -0.02  0.16 -0.06 -0.08
[3,] -0.03 -0.06  0.21 -0.12
[4,] -0.04 -0.08 -0.12  0.24
>  # third is 3x3, pos semi-def, but only rank 1
>   A3 <- matrix(1,3,3)
>   A3
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    1    1    1
[3,]    1    1    1
>  #
>  # usual function
>   chol(A1)
     [,1]      [,2]       [,3]
[1,]    3 0.6666667 -0.6666667
[2,]    0 0.7453560  0.5962848
[3,]    0 0.0000000  1.7888544
>  # expressions below not executed since they will stop execution
>  #  of batch job -- try these out in interactive session
>  #
>  #
>   #v2 <- chol(A2)
>   #v2
>   #v3 <- chol(A3)
>   #v3
>  #
>  #
>  # first fix -- just throw a warning flag
>   cholw <- function(A) {
+     tryCatch(chol(A),error=function(e){print("not good")})  }
>  #           expression    what to do with error
>  #
>  # notice response to the problem cases
>  # 
>  #
>  # no problem
>   cholw(A1)
     [,1]      [,2]       [,3]
[1,]    3 0.6666667 -0.6666667
[2,]    0 0.7453560  0.5962848
[3,]    0 0.0000000  1.7888544
>  # usually gives error and stops
>   cholw(A2)
[1] "not good"
>  # also gives error and stops
>   cholw(A3)
[1] "not good"
>  # fancy fix -- try to fit one rank smaller
>   chol1 <- function(A) {
+     tryCatch(chol(A),error=function(e){ # fix if last row/col dependent
+                       print(e)          # usual error message
+                       print("trying one smaller")
+                       lA1 <- nrow(A)-1            # one smaller
+                       R1 <- chol(A[1:lA1,1:lA1])  # Cholesky factor
+                       cbind(R1,forwardsolve(t(R1),A[1:lA1,lA1+1]))})  }
>  # use it on the three cases
>  #
>  # first one should give no change
>    chol1(A1)
     [,1]      [,2]       [,3]
[1,]    3 0.6666667 -0.6666667
[2,]    0 0.7453560  0.5962848
[3,]    0 0.0000000  1.7888544
>  # second one is the interesting one
>    v2 <- chol1(A2)
<simpleError in chol.default(A): the leading minor of order 4 is not positive definite>
[1] "trying one smaller"
>    v2
     [,1]        [,2]       [,3]       [,4]
[1,]  0.3 -0.06666667 -0.1000000 -0.1333333
[2,]  0.0  0.39440532 -0.1690309 -0.2253745
[3,]  0.0  0.00000000  0.4140393 -0.4140393
>    t(v2)%*%v2
      [,1]  [,2]  [,3]  [,4]
[1,]  0.09 -0.02 -0.03 -0.04
[2,] -0.02  0.16 -0.06 -0.08
[3,] -0.03 -0.06  0.21 -0.12
[4,] -0.04 -0.08 -0.12  0.24
>  # third one will still fail -- deficient in rank by 2
>    chol1(A3)
<simpleError in chol.default(A): the leading minor of order 2 is not positive definite>
[1] "trying one smaller"