【How Kalman Fliters Works，卡尔曼滤波（三）附录1：根据协方差矩阵进行随机抽样】

Y0UNG-code

已于 2024-07-22 15:44:36 修改

阅读量637

点赞数 23

分类专栏：卡尔曼滤波文章标签：矩阵线性代数笔记

于 2024-07-22 11:02:51 首次发布

本文链接：https://blog.csdn.net/m0_46626448/article/details/140602878

版权

卡尔曼滤波专栏收录该内容

4 篇文章 0 订阅

订阅专栏

本文翻译自本文翻译自How Kalman Filters Work，考虑到专业水平有限，如有问题和建议，欢迎在评论区讨论！

在对上一节中的粒子滤波器进行仿真时，我们需要以根据的协方差矩阵 $Q$ 进行随机抽样。
如何做到这一点还不是很清楚，所以让我们从现有的开始：任何相关的编程语言都有生成具有高斯分布和单位方差的随机数/向量的库。例如，在MATLAB我们可以用x1 = randn(3，1) 生成一个 $3\times1$ 的随机向量，也可以用x = randn(3, 1000)一次性生成 1000 个 $3\times1$ 的随机向量。randn函数产生的随机数具有单位方差，即
$E(\bm{x}\bm{x}^{T}) = I$ 其中 $I$ 表示单位矩阵。（在不面向矩阵的语言中，我们可以用一维高斯随机数发生器的独立抽样数来填充二维数组，这样也能得到相同的结果。）

Cholesky分解

现在我们希望抽样得到的随机数的协方差矩阵是 $Q$ ，该怎么做？根据Cholesky分解，我们可以解得一个矩阵 $C$ ，满足：
$Q = C C^{T}$ 所以 $C x$ 的协方差矩阵是：
$E(C\bm{x}\bm{x}^{T}C^{T}) = C E(\bm{x}\bm{x}^{T}) C^{T} = CIC^{T} = CC^{T} = 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}Pz \geq 0$ ），但Cholesky只适用于严格的正定矩阵（ $z^{T}Pz > 0$ ）。因此，实际上有更好的方法来计算 $C$ ，即使在半正定的情况下也行得通。

奇异值分解

使用奇异值分解（SVD）来计算，它将矩阵 $Q$ 分解成三个部分，即 $U$ 、 $S$ 和 $V$ ，满足：
$Q = USV^{T}$ 其中 $S$ 是对角矩阵，且对角线上的值大于等于0。对协方差矩阵进行奇异值分解时，可以发现 $U = V$ ，所以 $Q = USU^{T}$ 。

由于 $S$ 是对角矩阵，因此计算其平方根 $S^{1/2}$ 非常简单：只需要分别求出对角线上每个值的平方根。然后我们就可以构造 $C$ ：
$C = US^{1/2}$
检验一下：
$CC^{T} = US^{1/2}S^{1/2^{T}}U^{T} = USU^{T} = 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