Cholesky with singular matrices

[therneau]

I'm the author/maintainer of the survival package in R, which uses generalized Cholesky of the Hessian matrix (H= LDL' with D diagonal) for computation and am interested in switching over to LAPACK routines. The current code has been rock solid for 2 decades, this is motivated by the possible speedup, and even more by an in-development companion library for interval censored models that needs all the speedup it can get.

It is common in statistics for a user to inadvertently create a prediction matrix X that has a redundant column (sometimes more than one!). My routines simply set that element of D to zero and that row of L to a 1 on the diagonal. This looks to be similar to DSYTRF. The part I am not sure of is creating the Newton-Raphson coefficient update and the final variance matrix H-inverse. In both of these my routines simply skip those columns: that element of the solution b will be 0, and likewise that row/col of H-inverse. We don't want to treat these as "errors" in the computation, but rather handle them at the very end when packaging results for the user.

Questions:
Are there versions of the routines which take this path?
If not, is my best bet to create a modified versions of DSYTRS and DSYTRI that replace the "fail if D[j,j]=0" step with "keep going"?
Is there prior discussion of this issue that is pertinent? I am new to this forum so could easily have missed one.

Thanks in advance for any pointers
Terry Therneau. [email protected]

[langou]

I'm the author/maintainer of the survival package in R

Thanks for your work!

Are there versions of the routines which take this path?

LAPACK does not have a Cholesky factorization routine that "skips" a diagonal entry when it is smaller than a threshold. This is a good request though, there could be other applications that would benefit from it.

I would not recommend to use DSYTRF for this application.

Something that could be used is DPSTRF but DPSTRF does more than what you want and so might be unnecessarily slow.

If not, is my best bet to create a modified versions of DSYTRS and DSYTRI that replace the "fail if D[j,j]=0" step with "keep going"?

Yes, that's a good idea but start from DPOTRF for the factorization and adapt it to treat some diagonal entries. (Then yes you will also need to adapt the associated TRS and TRI.)

Is there prior discussion of this issue that is pertinent?

It is pertinent.