Computational complexity and implementation of UDU Modified Cholesky Rank 1 Update

Question

I am attempting to increase the performance of a legacy Kalman Filter implementation. The state covariance is factored in terms of UDU, i.e. $\mathbf{P} = \mathbf{U}\mathbf{D}\mathbf{U}^T$.

Many years ago, Bierman showed how scalar measurements may be used to update the covariance. Helpfully, it includes an implementation in FORTRAN that I have reproduced in MATLAB (though it combines certain operations in a loop and thus the algorithm is not particularly transparent).

Ultimately, I would like to replace the routine to take advantage of modern matrix libraries (BLAS, Eigen, etc) to increase performance.

Ultimately, the core of the operation is to perform a rank-1 update (downdate) of the UD-factored covariance in the form of:

$\mathbf{U_+}\mathbf{D_+}\mathbf{U_+}^T = \mathbf{U_-}\mathbf{D_-}\mathbf{U_-}^T - c\mathbf{v}\mathbf{v}^T$

where $\mathbf{U_-}$ and $\mathbf{D_-}$ is my prior known factorisation, c is a non-negative scalar and $\mathbf{v}$ is (known) a column vector. The goal is to determine $\mathbf{U_+}$ and $\mathbf{D_+}$.

So, my questions are:

• For the modified Cholesky update, what is the cost in terms of FLOPS - I understand it is $O(n^2$)?
• Is the a known efficient implementation in a commonly available library, or so I need to try to roll my own?

(For interest, $n$ is about 40 , but is being implemented on an embedded platform).

Answer 1

Yes, rank-one updates of a Cholesky factorization take $O(n^{2})$ time.

Unfortunately, there isn't a low rank update for the Cholesky factorization in LAPACK. There is an update routine for the Cholesky factorization in Linpack, but that routine works with the convention $LL^{T}$ factorization, not the square root free Cholesky factorization.

Unfortunately, there isn't enough data reuse in the process of rank-one updating a Cholesky factorization to get the kind of performance boost from blocknig and cache reuse that you get in various $O(n^3)$ operations that operate on $O(n^{2})$ data. .

Answer 2

Since a UDU rank 1 update is hard to find, below is a MATLAB implementation (without warranties of any kind)

Downdate

function [Ubar, Dbar] = ThorntonRank1Downdate(U, D, c, a)

Ubar = U;
Dbar = D;

n = length(a);

for j = n : -1 : 2

    s = a(j);
    d = D(j) - c * s^2;
    b = c / d;
    beta = s * b;
    c = b * D(j);
    Dbar(j) = d;

    a( 1:(j-1) ) = a( 1:(j-1) ) - s * U( 1:(j-1) ,j);
    Ubar( 1:(j-1) ,j) = U( 1:(j-1) ,j) - beta * a( 1:(j-1) ) ;

end

Dbar(1) = D(1) - c * a(1)^2;

Update

function [Ubar, Dbar] = ThorntonRank1Update(U, D, c, a)

Ubar = U;
Dbar = D;

n = length(a);

for j = n : -1 : 2

    s = a(j);
    d = D(j) + c * s^2;
    b = c / d;
    beta = s * b;
    c = b * D(j);
    Dbar(j) = d;

    for i = j-1 : -1 : 1

        a(i) = a(i) - s * U(i,j);
        Ubar(i,j) = U(i,j) + beta * a(i);

    end

end

Dbar(1) = D(1) + c * a(1)^2;

Test Program

% Requires symbolic toolbox
syms u12 u13 u14 u23 u24 u34 real
U = [1 u12 u13 u14; 0 1 u23 u24; 0 0 1 u24; 0 0 0 1]
% U = [1 u12 u13; 0 1 u23; 0 0 1]

syms v1 v2 v3 v4 real
v = [v1;v2;v3; v4]
% v = [v1;v2;v3]

syms d1 d2 d3 d4 real
D = diag([d1,d2,d3,d4]);
% D = diag([d1, d2, d3]);

syms c positive

%--------------------------------------------------------------------------

P = U * D * U.' + c * (v * v.');
[Udash, Ddash] = UD(P)

[Ubar, Dbar] = ThorntonRank1Update(U, diag(D), c, v)

update_TestU = simplify(Udash - Ubar, 'seconds', 10)
update_TestD = simplify(Ddash - diag(Dbar), 'seconds', 10)

update_TestPbar  = simplify((Ubar * diag(Dbar) * Ubar.') - P, 'seconds', 10)
update_TestPdash = simplify((Udash * Ddash * Udash.') - P, 'seconds', 10)

%--------------------------------------------------------------------------

P = U * D * U.' - c * (v * v.');
[Udash, Ddash] = UD(P)

[Ubar, Dbar] = ThorntonRank1Downdate(U, diag(D), c, v)

downdate_TestU = simplify(Udash - Ubar, 'seconds', 10)
downdate_TestD = simplify(Ddash - diag(Dbar), 'seconds', 10)

downdate_TestPbar  = simplify((Ubar * diag(Dbar) * Ubar.') - P, 'seconds', 10)
downdate_TestPdash = simplify((Udash * Ddash * Udash.') - P, 'seconds', 10)



%--------------------------------------------------------------------------

Reference

Page 55 of: Bierman, G. J., "Factorization Methods for Discrete Sequential Estimation", Academic Press, 1977

Also useful was the following paper:

Fletcher, R. & Powell, M. J. D. On the modification of LDLT factorizations Mathematics of Computation, 1974, 28, 1067-1087 (PDF)

Complexity

Without doing a proper FLOPS count, it looks to be in the order of $n^2 + O(n)$.