4
$\begingroup$

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).

$\endgroup$
3
$\begingroup$

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

$\endgroup$
  • $\begingroup$ So, it looks like I need to hand-roll a routine. Is there a good reference for this? Any techniques I should watch out for to try to maximise performance? $\endgroup$ – Damien Aug 27 '13 at 5:41
  • $\begingroup$ Also, is there a reference for a FLOPS count (not just $O(n^2)$ ) for the cholupdate - I haven't managed to find one so far. $\endgroup$ – Damien Aug 27 '13 at 5:41
  • $\begingroup$ The references that I have at hand are all for updating the LDLT factorization rather than the UDUT factorization. You could rederive this from the Sherman-Morison-Woodbury formula, but I expect that some library work will help you find a version for UDUT. The flops count should be easy to get from the actual code that you use (and you can easily instrument the code to explicitly count the flops used if you want to check your formula.) $\endgroup$ – Brian Borchers Aug 27 '13 at 20:26
1
$\begingroup$

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)$.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.