Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Full rank Cholesky updates? #14

Open
pnkraemer opened this issue Feb 26, 2021 · 0 comments
Open

Full rank Cholesky updates? #14

pnkraemer opened this issue Feb 26, 2021 · 0 comments

Comments

@pnkraemer
Copy link

In ProbNum we use full-rank Cholesky updates (see the utility function probnum.statespace.cholesky_update) for square-root forward transitions.

Do you want the utility functions in here? For starters, we could merge the code (https://github.com/probabilistic-numerics/probnum/blob/master/src/probnum/filtsmooth/statespace/discrete_transition_utils.py)

def cholesky_update(
    S1: np.ndarray, S2: typing.Optional[np.ndarray] = None
) -> np.ndarray:
    r"""Compute Cholesky update/factorization :math:`C C^\top = S_1 S_1^\top + S_2 S_2^\top`.
    This can be used in various ways.
    For example, :math:`S_1` and :math:`S_2` do not need to be Cholesky factors; any matrix square-root is sufficient.
    As long as :math:`C C^\top = S_1 S_1^\top + S_2 S_2^\top` is well-defined (and admits a Cholesky-decomposition),
    :math:`S_1` and :math:`S_2` do not even have to be square.
    """
    # doc might need a bit more explanation in the future
    # perhaps some doctest or so?
    if S2 is not None:
        stacked_up = np.vstack((S1.T, S2.T))
    else:
        stacked_up = np.vstack(S1.T)
    upper_sqrtm = np.linalg.qr(stacked_up, mode="r")
    return triu_to_positive_tril(upper_sqrtm)


def triu_to_positive_tril(triu_mat: np.ndarray) -> np.ndarray:
    r"""Change an upper triangular matrix into a valid lower Cholesky factor.
    Transpose, and change the sign of the diagonals to '+' if necessary.
    The name of the function is leaned on `np.triu` and `np.tril`.
    """
    tril_mat = triu_mat.T
    d = np.sign(np.diag(tril_mat))
    d[d == 0] = 1.0
    with_pos_diag = tril_mat @ np.diag(d)
    return with_pos_diag

with its respective tests (subject to interface changes to remain compliant with your current implementation of rank-1 updates) as soon as ProbNum has an explicit dependency on cholupdates anyway.

This could then be used inside random_variables as well as inside statespace, and optimisation of this functionality can be outsourced to optimising cholupdates.
Note that the cholesky_update function above already handles your Issue #13 in a basic way (more optimized code seems possible, but can be a future project?!).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant