本文翻译自本文翻译自How Kalman Filters Work,考虑到专业水平有限,如有问题和建议,欢迎在评论区讨论!
在对上一节中的粒子滤波器进行仿真时,我们需要以根据的协方差矩阵
Q
Q
Q进行随机抽样。
如何做到这一点还不是很清楚,所以让我们从现有的开始:任何相关的编程语言都有生成具有高斯分布和单位方差的随机数/向量的库。例如,在MATLAB我们可以用x1 = randn(3,1)
生成一个
3
×
1
3\times1
3×1的随机向量,也可以用x = randn(3, 1000)
一次性生成 1000 个
3
×
1
3\times1
3×1的随机向量。randn
函数产生的随机数具有单位方差,即
E
(
x
x
T
)
=
I
E(\bm{x}\bm{x}^{T}) = I
E(xxT)=I其中
I
I
I表示单位矩阵。(在不面向矩阵的语言中,我们可以用一维高斯随机数发生器的独立抽样数来填充二维数组,这样也能得到相同的结果。)
Cholesky分解
现在我们希望抽样得到的随机数的协方差矩阵是
Q
Q
Q,该怎么做?根据Cholesky分解,我们可以解得一个矩阵
C
C
C,满足:
Q
=
C
C
T
Q = C C^{T}
Q=CCT所以
C
x
Cx
Cx的协方差矩阵是:
E
(
C
x
x
T
C
T
)
=
C
E
(
x
x
T
)
C
T
=
C
I
C
T
=
C
C
T
=
Q
E(C\bm{x}\bm{x}^{T}C^{T}) = C E(\bm{x}\bm{x}^{T}) C^{T} = CIC^{T} = CC^{T} = Q
E(CxxTCT)=CE(xxT)CT=CICT=CCT=Q
因此,要获得具有协方差 𝑄 的随机抽样,我们先从具有单位方差的随机抽样开始,然后乘以 𝑄 的 Cholesky 因子。
这里是 MATLAB 代码:
C = chol(Q, 'lower'); % Make sure to use *lower* Cholesky factor!
x = randn(nx, n); % Make n random nx-dimensioned vectors.
y = C * randn(nx, n); % Give them the proper covariance.
还有一个细节是,我们在前面中提到协方差矩阵必须是半正定矩阵( z T P z ≥ 0 z^{T}Pz \geq 0 zTPz≥0),但Cholesky只适用于严格的正定矩阵( z T P z > 0 z^{T}Pz > 0 zTPz>0)。因此,实际上有更好的方法来计算 C C C,即使在半正定的情况下也行得通。
奇异值分解
使用奇异值分解(SVD)来计算,它将矩阵
Q
Q
Q分解成三个部分,即
U
U
U、
S
S
S和
V
V
V,满足:
Q
=
U
S
V
T
Q = USV^{T}
Q=USVT其中
S
S
S是对角矩阵,且对角线上的值大于等于0。对协方差矩阵进行奇异值分解时,可以发现
U
=
V
U = V
U=V,所以
Q
=
U
S
U
T
Q = USU^{T}
Q=USUT。
由于
S
S
S是对角矩阵,因此计算其平方根
S
1
/
2
S^{1/2}
S1/2非常简单:只需要分别求出对角线上每个值的平方根。然后我们就可以构造
C
C
C:
C
=
U
S
1
/
2
C = US^{1/2}
C=US1/2
检验一下:
C
C
T
=
U
S
1
/
2
S
1
/
2
T
U
T
=
U
S
U
T
=
Q
CC^{T} = US^{1/2}S^{1/2^{T}}U^{T} = USU^{T} = Q
CCT=US1/2S1/2TUT=USUT=Q
这里是 MATLAB 代码:
[U, S] = svd(Q); % Works even when Q is positive-semidefinite.
C = U * sqrt(S); % Construct the matrix square root.
x = randn(nx, n); % Make n random nx-dimensioned vectors.
y = C * randn(nx, n); % Give them the proper covariance.
打包成一个函数(matlab)——用奇异值分解
function c = covdraw(C, n)
% Draw n random samples from a zero-mean Gaussian distribution with covariance C.
% 从协方差为 C 的零均值高斯分布中随机抽取 n 个样本。
% Inputs:
%
% C Covariance matrix (nC-by-nC) 协方差矩阵,nC:变量个数
% n Number of samples to draw (default is 1) 需要抽取的样本数(默认为1)
%
% Outputs:
%
% c Draws from C (nC-by-n) 抽取的样本
%
% Defaults
if nargin < 2, n = 1; end;
[u, s] = svd(C);
s(s < 0) = 0; % Nothing should be negative, but just in case... 应该不会出现负数,但以防万一
s = sqrt(s);
% At this point, C == (u * s) * (u * s).', so (u * s) is the Cholesky
% factor (if C was positive definite).
% Construct the draws. 构造采样
c = u * s * randn(size(C, 1), n);
end % covdraw