使协方差矩阵的平方根正定（Matlab）

机器算法验证 matlab 协方差矩阵 数字
2022-03-29 16:34:24

动机：我正在 MATLAB 中编写状态估计器（无迹卡尔曼滤波器），它要求更新协方差矩阵的（上三角）平方根S在每次迭代中（即，对于协方差矩阵P, 确实P=SST）。为了让我执行必要的计算，我需要使用 MATLAB 的cholupdate函数进行 Rank-1 Cholesky Update 和 Downdate。

问题：不幸的是，在迭代过程中，这个矩阵S有时会失去正定性。Cholesky 更新在非 PD 矩阵上失败。

我的问题是：在 MATLAB 中是否有任何简单可靠的方法来制作S正定？

（或更一般地说，有没有一种很好的方法来制作任何给定的协方差X矩阵正定？)

备注：

  • S满级
  • 我尝试了特征分解方法（不起作用）。这基本上涉及发现S=VDVT, 设置所有负元素V,D=1×10−8，并重建一个新的S′=V′D′V′T在哪里V′,D′是只有正元素的矩阵。
  • 我知道 Higham 方法（在 R as 中实现nearpd），但它似乎只投影到最近的 PSD 矩阵。我需要一个用于 Cholesky 更新的 PD 矩阵。
4个回答

这是我过去使用过的代码（使用 SVD 方法）。我知道你说你已经尝试过了，但它一直对我有用，所以我想我会发布它，看看它是否有帮助。

function [sigma] = validateCovMatrix(sig)

% [sigma] = validateCovMatrix(sig)
%
% -- INPUT --
% sig:      sample covariance matrix
%
% -- OUTPUT --
% sigma:    positive-definite covariance matrix
%

EPS = 10^-6;
ZERO = 10^-10;

sigma = sig;
[r err] = cholcov(sigma, 0);

if (err ~= 0)
    % the covariance matrix is not positive definite!
    [v d] = eig(sigma);

    % set any of the eigenvalues that are <= 0 to some small positive value
    for n = 1:size(d,1)
        if (d(n, n) <= ZERO)
            d(n, n) = EPS;
        end
    end
    % recompose the covariance matrix, now it should be positive definite.
    sigma = v*d*v';

    [r err] = cholcov(sigma, 0);
    if (err ~= 0)
        disp('ERROR!');
    end
end

在 Matlab 中：

help cholupdate

我明白了

CHOLUPDATE Rank 1 update to Cholesky factorization.
    If R = CHOL(A) is the original Cholesky factorization of A, then
    R1 = CHOLUPDATE(R,X) returns the upper triangular Cholesky factor of A + X*X',
    where X is a column vector of appropriate length.  CHOLUPDATE uses only the
    diagonal and upper triangle of R.  The lower triangle of R is ignored.

    R1 = CHOLUPDATE(R,X,'+') is the same as R1 = CHOLUPDATE(R,X).

    R1 = CHOLUPDATE(R,X,'-') returns the Cholesky factor of A - X*X'.  An error
    message reports when R is not a valid Cholesky factor or when the downdated
    matrix is not positive definite and so does not have a Cholesky factorization.

    [R1,p] = CHOLUPDATE(R,X,'-') will not return an error message.  If p is 0
    then R1 is the Cholesky factor of A - X*X'.  If p is greater than 0, then
    R1 is the Cholesky factor of the original A.  If p is 1 then CHOLUPDATE failed
    because the downdated matrix is not positive definite.  If p is 2, CHOLUPDATE
    failed because the upper triangle of R was not a valid Cholesky factor.

    CHOLUPDATE works only for full matrices.

    See also chol.

计算 Cholesky 分解的另一种方法是将 S 的对角元素固定为 1，然后引入具有正元素的对角矩阵 D。

这避免了在进行计算时需要取平方根，这在处理“小”数字时可能会导致问题（即数字足够小，以至于由于浮点运算而发生的舍入很重要）。维基百科页面有这个调整后的算法的样子。

所以而不是P=SST你得到P=RDRT和S=RD12

希望这可以帮助！

当您的矩阵不是“真正”确定时，Cholesky 分解可能会失败。出现两种情况，或者你的本征值是负的，或者你最小的本征值是正的，但接近于零。第二种情况必须在理论上给出一个解决方案，但在数值上是困难的。如果只是直观地在我的矩阵的对角线上添加一个小常数来解决问题。但是这种方式并不严谨，因为它稍微修改了解决方案。如果您必须计算一个非常高精度的解决方案，请尝试对改进的 Cholesky 分解进行一些研究。

其它你可能感兴趣的问题
上一篇如何报告比例的不对称置信区间？ 下一篇评估相关性的重要性