卡尔曼滤波器实现中,"预测估计协方差"P(k|k-1)(参见此处的维基)是否有可能是一个奇异矩阵?如果不是,我的代码有问题吗?
这是状态空间模型
% y{t}=Z{t} b{t} + eps{t}, eps{t} ~ N(0,H{t})
% b{t} = Pi{t} b{t-1} + tao{t} tao{t} ~ N(0,Q{t})
% b{1} ~ N(b0,P0)
% t=1,...,T
这是作为卡尔曼滤波算法主要部分的向后递归:
for t=1:T
v{t} = y{t} - Z{t} * b_tt_1{t};
M{t} = P_tt_1{t} * Z{t}';
F{t} = Z{t} * M{t} + H{t};
F_{t}= inv(F{t});
MF_{t}= M{t} * F_{t};
b_tt{t}=b_tt_1{t} + MF_{t} * v{t};
P_tt{t}=P_tt_1{t} - MF_{t} * M{t}';
b_tt_1{t+1} = Pi{t} * b_tt{t};
P_tt_1{t+1} = Pi{t} * P_tt{t} * Pi{t}' + Q{t};
end
当我使用实际数据时,就会发生这种情况。为了查看问题可能出在哪里,我编写了一些代码来生成随机状态空间模型(如果需要,我可以提供代码)。
当 T 很大时,在某个 t0 之后,P_tt_1{t0} 是奇异的,状态 (b{t0}) 发散。
编辑:我使用了协方差更新方程的"约瑟夫形式"(见维基百科)。它有所帮助,但是当状态空间模型很大(在方程或状态的数量意义上)时,结果仍然发散。我认为这意味着问题与数值稳定性有关。有没有办法解决这个问题?
矩阵唯一可以变成单数的地方是在以下行中:
F_{t}= inv(F{t});
你可以改用伪逆的"pinv"。
如果你重写这些行,那就更好了:
F_{t}= inv(F{t});
MF_{t}= M{t} * F{t};
自
MF_{t}= M{t} / F{t};
Matlab 将求解线性方程:MF_{t} * F{t} = M{t}
- 即使 F_{t} 是奇异的,也可能有一个解 - 或者如果它仍然是通过伪逆求解的奇异。
问题出在 inv 或 pinv 函数的反演过程中。 矩阵很大,但它们也是正定的。 所以我使用cholesky分解进行反演。函数为:
function A_=inversePD(A)
%A:positive definite matrix
M=size(A,1);
[R b] = chol(A);
if b~=0
return
end
R_ = R eye(M);
A_ = R_ * R_';
end