Posting link here because it is as much math as it is programming:
(Feel free to delete if inappropriate)
(Apologies for having lapsed into abject mathematical ignorance and incompetency over the last few decades).
Thank you for getting this far!
From the documentation of the function I'm trying to replicate:
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.
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 108−111). 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: (images that I can't really put in here, see the documentation here)
I want to make sure I really understand what is being solved for here. Cholesky factorization wants a positive definite matrix, so we're solving for the smallest modifier that does that? And after that, it's merely a matter of algorithm?