我正在尝试提高传统卡尔曼滤波器实现的性能。状态协方差是根据 UDU 考虑的，即.
许多年前，Bierman展示了如何使用标量测量来更新协方差。有用的是，它包括我在 MATLAB 中复制的 FORTRAN 中的实现（尽管它将某些操作组合在一个循环中，因此算法不是特别透明）。
最终，我想替换例程以利用现代矩阵库（BLAS、Eigen 等）来提高性能。
最终，操作的核心是执行 UD 因子协方差的 rank-1 更新（更新），形式为：
其中和是我之前已知的分解，c 是非负标量，是（已知的）列向量。目标是确定和。
所以，我的问题是：
（出于兴趣，约为 40 ，但正在嵌入式平台上实现）。
是的，Cholesky 分解的一级更新需要时间。
不幸的是，LAPACK 中的 Cholesky 分解没有低等级更新。Linpack 中有一个用于 Cholesky 分解的更新例程，但该例程适用于约定分解，而不是无平方根的 Cholesky 分解。
不幸的是，在秩一更新 Cholesky 分解的过程中没有足够的数据重用来从块尼格和缓存重用中获得性能提升，就像在 O(n) 上运行的各种操作中数据。.
由于 UDU 等级 1 更新很难找到，下面是一个 MATLAB 实现（没有任何形式的保证）
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;
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;
% 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)
%--------------------------------------------------------------------------
第 55 页：Bierman, GJ，“离散序列估计的因式分解方法”，学术出版社，1977 年
以下论文也很有用：
Fletcher, R. & Powell, MJD On the modify of LDLT factorisations of Computation, 1974, 28, 1067-1087 ( PDF )
如果没有进行适当的 FLOPS 计数，它看起来应该是.