cholupdates.rank_1.update

cholupdates.rank_1.update(L, v, check_diag=True, overwrite_L=False, overwrite_v=False, method='seeger', **method_kwargs)

Update a Cholesky factorization after addition of a positive-semidefinite symmetric rank-1 matrix.

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.\]

Parameters

L ((N, N) numpy.ndarray) – Lower-triangular Cholesky factor of \(A\). Must have a non-zero diagonal. 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) – The vector \(v\) with shape (N,) defining the symmetric rank-1 update matrix \(v v^T\).
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 may 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 may reuse the array v as an internal computation buffer, which will modify v. Passing False here ensures that the array v is not modified.
method (str) –
Algorithm to be used to compute the updated Cholesky factor. Must be one of
- ”cho_factor”
  Directly uses scipy.linalg.cho_factor() on \(L L^T + v v^T\). This is just here for convenience and should be slower than all other methods.
- ”seeger”
  Calls cholupdates.rank_1.update_seeger().
Defaults to “seeger”.
method_kwargs – Additional keyword arguments which will be passed to the function selected by method.

Returns

Lower triangular Cholesky factor \(L^+\) of \(A + v v^T\). 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.

Return type

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

Raises

ValueError – If L does not have shape (N, N) for some N.
numpy.linalg.LinAlgError – If the diagonal of L contains zeros and check_diag is set to True.
ValueError – If v does not have shape (N,), while L has shape (N, N).
Exception – Any exception raised by the function specified by method.

See also

cholupdates.rank_1.downdate: A similar function which performs a symmetric rank 1 downdate instead of an update.

Examples

Consider the following matrix-vector pair

>>> A = np.diag([1.0, 2.0, 3.0]) + 0.1
>>> A
array([[1.1, 0.1, 0.1],
       [0.1, 2.1, 0.1],
       [0.1, 0.1, 3.1]])
>>> v = np.array([1.0, 25.0, 10.0])
>>> v
array([ 1., 25., 10.])

We want to compute the lower triangular Cholesky factor L_ud of

>>> A_ud = A + np.outer(v, v)
>>> A_ud
array([[  2.1,  25.1,  10.1],
       [ 25.1, 627.1, 250.1],
       [ 10.1, 250.1, 103.1]])

We assume that the lower triangular Cholesky factor of A is given

>>> import scipy.linalg
>>> L = scipy.linalg.cho_factor(A, lower=True)[0]
>>> np.tril(L)
array([[1.04880885, 0.        , 0.        ],
       [0.09534626, 1.44599761, 0.        ],
       [0.09534626, 0.06286946, 1.75697368]])

The function cholupdates.rank_1.update() can compute L_ud from L efficiently

>>> import cholupdates
>>> L_ud = cholupdates.rank_1.update(L, v, method="seeger")
>>> np.tril(L_ud)
array([[ 1.44913767,  0.        ,  0.        ],
       [17.32064554, 18.08577447,  0.        ],
       [ 6.96966215,  7.15374133,  1.82969791]])

Did it work?

>>> np.allclose(A_ud, np.tril(L_ud) @ np.tril(L_ud).T)
True

We could also compute L_ud by directly computing the Cholesky factorization of A_ud (which is however less efficient)

>>> L_ud_cho = cholupdates.rank_1.update(L, v, method="cho_factor")
>>> np.tril(L_ud_cho)
array([[ 1.44913767,  0.        ,  0.        ],
       [17.32064554, 18.08577447,  0.        ],
       [ 6.96966215,  7.15374133,  1.82969791]])
>>> np.allclose(np.tril(L_ud), np.tril(L_ud_cho))
True