Bfgs Hessian direction by failure-aware Cholesky solve

Input a real symmetric positive-definite matrix $B\in {\mathbb{R}}^{n\times n}$, $n\ge 1$, a vector $g\in {\mathbb{R}}^n$, and representable finite entries in the operationwise arithmetic model of Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add. Output: either a computed direction $p$ or a failure status.

1. Call Failure-aware rounded Cholesky factorization on $B$. If it returns a failure status, return factorization-failure with that pivot and operation index.

2. If it returns a lower-triangular factor $\widehat{L}$ with positive normal diagonal, implement the recursions of Forward and backward substitution for triangular linear systems as follows, with indices $1,\dots ,n$. For the forward solve, process rows $i=1,\dots ,n$ in increasing order. Initialize $s:=0$; for $j=1,\dots ,i-1$ in increasing order, first compute $t:=\mathrm{fl}({\widehat{l}}_{ij}{z}_j)$, then $s:=\mathrm{fl}(s+t)$. After this inner loop, compute $r:=\mathrm{fl}((-g_i)-s)$, then $z_i:=\mathrm{fl}(r\mathrm{/}{\widehat{l}}_{ii})$. For the backward solve, process rows $i=n,\dots ,1$ in decreasing order. Initialize $s:=0$; for $j=i+1,\dots ,n$ in increasing order, first compute $t:=\mathrm{fl}({\widehat{l}}_{ji}{p}_j)$, then $s:=\mathrm{fl}(s+t)$. After this inner loop, compute $r:=\mathrm{fl}(z_i-s)$, then $p_i:=\mathrm{fl}(r\mathrm{/}{\widehat{l}}_{ii})$. Empty inner loops leave $s=0$, and sign changes are exact. Each displayed $\mathrm{fl}$ is a separate scalar operation in the given arithmetic model; do not reassociate or fuse operations. On the first range-failure, stop immediately and return triangular-solve-failure with the solve (forward or backward), row $i$, inner index $j$ when applicable, and operation (multiply, add, subtract, or divide).

3. Otherwise return $p$. The exact-real reference problem is $Bp=-g$ and has a unique solution by Cholesky solves Hermitian positive-definite systems and has about half the factorisation cost of LU. The returned vector is the computed solve from $\widehat{L}$; this definition asserts no backward-error bound for it.

The procedure performs a finite factorization loop and two finite triangular loops. Its use in the full BFGS algorithm includes a subsequent descent check; if the computed direction is not descending, that algorithm resets to a positive scalar multiple of the negative gradient.