cholupdates.rank_1.update_seeger

cholupdates.rank_1.update_seeger(L, v, check_diag=True, overwrite_L=False, overwrite_v=False, impl=None)

Update a Cholesky factorization after addition of a positive semidefinite symmetric rank-1 matrix using the algorithm from section 2 of [Seeger, 2008].

In other words, given \(A = L L^T \in \mathbb{R}^{N \times N}\) and \(v \in \mathbb{R}^N\), compute \(L^+\) such that

\[A^+ := A + v v^T = L^+ (L^+)^T.\]

This function computes the Cholesky factorization of \(A^+\) from \(L\) in \(O(N^2)\) time, which is faster than the \(O(N^3)\) time complexity of naively applying a Cholesky algorithm to \(A^+\) directly.

Parameters

  • L ((N, N) numpy.ndarray, dtype=numpy.single or numpy.double) – Lower-triangular Cholesky factor of \(A\). Must have a non-zero diagonal. Must have the same dtype as v. Must be a contiguous array (either Fortran- or C-contiguous). The algorithm is most efficient if this array is given in column-major layout, a.k.a. Fortran-contiguous or f-contiguous memory order. Hint: Lower-triangular Cholesky factors in column-major memory layout can be obtained efficiently (i.e. without requiring an additional copy) from scipy.linalg.cho_factor(). The entries in the strict upper-triangular part of L can contain arbitrary values, since the algorithm neither reads from nor writes to this part of the matrix. This behavior is useful when using the Cholesky factors returned by scipy.linalg.cho_factor() which contain arbitrary values on the irrelevant triangular part of the matrix.

  • v ((N,) numpy.ndarray, dtype=numpy.single or numpy.double) – The vector \(v\) which defines the symmetric rank-1 update matrix \(v v^T\). Must have the same dtype as L. Must be a contiguous array.

  • check_diag (bool) – If set to True, the function will check whether the diagonal of the given Cholesky factor L is non-zero and raise a ValueError if this is not the case. Setting check_diag to False can be used to speed up computations if it is clear that the Cholesky factor can not have zeros on its diagonal. Caution: If this argument is set to False and the Cholesky factor does contain zeros on its diagonal, the behavior of the function will be undefined.

  • overwrite_L (bool) – If set to True, the function will overwrite the array L with the upper Cholesky factor \(L^+\) of \(A^+\), i.e. the result is computed in-place. Passing False here ensures that the array L is not modified.

  • overwrite_v (bool) – If set to True, the function will reuse the array v as an internal computation buffer, which will modify v. Passing False here ensures that the array v is not modified. In this case, an additional array of shape (N,) and dtype numpy.double must be allocated.

  • impl (Optional[str]) –

    Defines which implementation of the algorithm to use. Must be one of

    • None

      Choose the Cython implementation if it is available, otherwise use the Python implementation.

    • ”cython”

      Use the Cython implementation. Throws a ValueError if the Cython implementation is not available.

    • ”python”

      Use the Python implementation.

    Defaults to None.

Returns

Lower-triangular Cholesky factor \(L^+\) of \(A + v v^T\) with shape (N, N). The diagonal entries of this matrix are guaranteed to be positive. The strict upper-triangular part of this matrix will contain the values from the upper-triangular part of L. The matrix will inherit dtype and memory order from L.

Return type

(N, N) numpy.ndarray, dtype=L.dtype

Raises

See also

cholupdates.rank_1.update

Interface function for all rank-1 update methods. Can be used to call this function conveniently.

Notes

This function implements the algorithm from 1, section 2. In the following, we will briefly summarize the theory behind the algorithm. Let \(A \in \mathbb{R}^{n \times n}\) be symmetric and positive-definite, and \(v \in \mathbb{R}^n\). The vector \(v\) defines a symmetric rank-1 update \(v v^T\) to \(A\), i.e.

\[A^+ := A + v v^T.\]

Assume that the Cholesky factorization of \(A\), i.e. a lower-triangular matrix \(L \in \mathbb{R}^{n \times n}\) with \(A = L L^T\), is given. We want to find a fast and stable method to update the Cholesky factorization \(L L^T\) of \(A\) to the Cholesky factorization of the updated matrix \(A^+\).

To this end, we rewrite the update equation as

\[\begin{split}A^+ = \underbrace{L L^T}_{= A} + v v^T = \begin{pmatrix} L & v \end{pmatrix} \begin{pmatrix} L^T \\ v^T \end{pmatrix}.\end{split}\]

We can now construct a sequence of orthogonal transformations \(Q := Q_1 \dotsb Q_n \in \mathbb{R}^{(n + 1) \times (n + 1)}\), in this case Givens rotations, such that

\[\begin{pmatrix} L & v \end{pmatrix} Q = \begin{pmatrix} L^+ & 0 \end{pmatrix},\]

where \(L^+ \in \mathbb{R}^{n \times n}\) is lower-triangular. In other words, \(Q\) eliminates the last column, i.e. \(v\), from the augmented matrix \((L \mid v)\), while preserving the lower-triangular structure of the left block. But now we have

\[\begin{split}L^+ (L^+)^T & = \begin{pmatrix} L^+ & 0 \end{pmatrix} \begin{pmatrix} L^+ & 0 \end{pmatrix}^T \\ & = \begin{pmatrix} L & v \end{pmatrix} \underbrace{Q Q^T}_{= I} \begin{pmatrix} L^T \\ v^T \end{pmatrix} \\ & = \underbrace{L L^T}_{= A} + v v^T \\ & = A^+,\end{split}\]

so \(L^+\) is already the lower-triangular Cholesky factor of \(A^+\).

As mentioned above, we need to multiply \(n\) Givens rotation matrices to \((L \mid v)\) which takes \(O(n)\) arithmetic operations for each rotation matrix. Hence, we end up with a total time complexity of \(O(n^2)\).

Note that this algorithm is a minor modification of the LINPACK 2 routine dchud. The exact modifications are described in 1.

References

1(1,2)

  1. Seeger, “Low Rank Updates for the Cholesky Decomposition”, 2008.

2

J. Dongarra, C. Moler, J. Bunch, and G. Steward. “LINPACK User’s Guide”. Society for Industrial and Applied Mathematics, 1979.