r"""In-tree rank-1 Cholesky update and downdate.
Replaces the external ``choldate`` package, which has no PyPI release, must be
installed from git with ``--no-build-isolation``, and therefore made the whole
package impossible to ``pip install`` cleanly.
Convention (identical to ``choldate``): ``R`` is **upper triangular** with
.. math::
A = R^\top R
``cholupdate(R, x)`` overwrites ``R`` in place with the factor of
:math:`A + x x^\top`; ``choldowndate(R, x)`` with the factor of
:math:`A - x x^\top`. Both use ``x`` as scratch space, so a copy is taken
internally and the caller's array is left untouched (``choldate`` clobbers it).
The algorithm is the standard LINPACK sequence of Givens rotations: one
rotation per row, each annihilating a component of ``x`` into the factor. It is
:math:`O(n^2)` and allocation-free.
**One deliberate difference from choldate.** On a downdate whose result is not
positive definite, ``choldate`` neither raises nor returns NaN -- it silently
computes :math:`\sqrt{|r^2|}` and returns a plausible-looking but wrong factor.
This implementation writes NaN instead, which is appropriate for the existing
``np.any(np.isnan(...))`` guard.
"""
__all__ = ["cholupdate", "choldowndate"]
import numpy as np
from numba import njit
@njit(cache=True)
def _cholupdate_impl(R: np.ndarray, x: np.ndarray) -> None:
n = R.shape[0]
for k in range(n):
Rkk = R[k, k]
xk = x[k]
r = np.sqrt(Rkk * Rkk + xk * xk)
if Rkk == 0.0:
# A structurally zero row. With x_k == 0 the Givens rotation is the
# identity; with x_k != 0 it is the swap (c = 0, s = sign(x_k)).
# Both are well defined, so do not propagate NaN.
if xk == 0.0:
continue
sign = 1.0 if xk > 0.0 else -1.0
R[k, k] = abs(xk)
for i in range(k + 1, n):
previous = R[k, i]
R[k, i] = sign * x[i]
x[i] = -sign * previous
continue
c = r / Rkk
s = xk / Rkk
R[k, k] = r
for i in range(k + 1, n):
R[k, i] = (R[k, i] + s * x[i]) / c
x[i] = c * x[i] - s * R[k, i]
@njit(cache=True)
def _choldowndate_impl(R: np.ndarray, x: np.ndarray) -> None:
n = R.shape[0]
for k in range(n):
Rkk = R[k, k]
xk = x[k]
r2 = Rkk * Rkk - xk * xk
if Rkk == 0.0 and xk == 0.0:
# A structurally zero row that the downdate does not touch.
continue
if Rkk == 0.0 or r2 <= 0.0:
# A - x x^T is not positive definite; there is no real factor.
# Fill with NaN so callers' finite-checks trip instead of
# silently accepting a wrong factor.
for i in range(k, n):
for j in range(n):
R[i, j] = np.nan
return
r = np.sqrt(r2)
c = r / Rkk
s = xk / Rkk
R[k, k] = r
for i in range(k + 1, n):
R[k, i] = (R[k, i] - s * x[i]) / c
x[i] = c * x[i] - s * R[k, i]
def _checked(R: np.ndarray, x: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
"""Validate shapes before entering the bounds-check-free compiled kernels."""
if not isinstance(R, np.ndarray) or R.ndim != 2 or R.shape[0] != R.shape[1]:
raise ValueError("R must be a square 2-D numpy array")
if R.dtype != np.float64:
raise TypeError(f"R must have dtype float64, got {R.dtype}")
x = np.ascontiguousarray(x, dtype=np.float64)
if x.shape != (R.shape[0],):
raise ValueError(f"x must have shape ({R.shape[0]},), got {x.shape}")
return R, x.copy()
[docs]
def cholupdate(R: np.ndarray, x: np.ndarray) -> None:
r"""Rank-1 update: overwrite ``R`` with the factor of :math:`A + x x^\top`.
:param R: Upper-triangular Cholesky factor, modified **in place**.
:type R: numpy.ndarray
:param x: Update vector. Not modified (a working copy is taken).
:type x: numpy.ndarray
:return: None
:rtype: None
"""
_cholupdate_impl(*_checked(R, x))
[docs]
def choldowndate(R: np.ndarray, x: np.ndarray) -> None:
r"""Rank-1 downdate: overwrite ``R`` with the factor of
:math:`A - x x^\top`.
If the downdated matrix is not positive definite, ``R`` is filled with NaN.
:param R: Upper-triangular Cholesky factor, modified **in place**.
:type R: numpy.ndarray
:param x: Downdate vector. Not modified (a working copy is taken).
:type x: numpy.ndarray
:return: None
:rtype: None
"""
_choldowndate_impl(*_checked(R, x))