已知协方差矩阵，基于Gram-Schmidt构造求随机向量新息表示的MATLAB实现

Got it, I totally get why you're stuck—most tutorials fixate on raw data matrices, but we absolutely don't need actual samples to do this. The Gram-Schmidt process relies on inner products (which is exactly what covariance entries are, since your random vector is zero-mean), so we can build the lower triangular matrix ( B ) directly from your 3×3 covariance matrix ( R ).

Let's break this down step by step for your 3-variable case, then generalize it (and give you MATLAB code to boot):

Core Idea #

We want ( \mathbf{y} = B\boldsymbol{\varepsilon} ), where:

• ( B ) is lower triangular
• ( \boldsymbol{\varepsilon} ) is a zero-mean vector of uncorrelated components (we'll make their variances 1 for simplicity, so ( E[\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}^T] = I ))

Each entry of ( B ) comes from projecting the corresponding ( y_i ) onto the space spanned by the previous ( y_1, ..., y_{i-1} ), then taking the residual norm—all using covariance values instead of raw data.

Step-by-Step Calculation for 3 Variables #

Let ( R[i,j] = E[y_i y_j] ):

1. First component:
   • ( y_1 = b_{11}\varepsilon_1 ). To make ( E[\varepsilon_1^2] = 1 ), we set ( b_{11} = \sqrt{R[1,1]} ).
2. Second component:
   • We need ( y_2 = b_{21}\varepsilon_1 + b_{22}\varepsilon_2 ), with ( E[\varepsilon_2\varepsilon_1] = 0 ).
   • ( b_{21} = \frac{R[2,1]}{b_{11}} ) (this is the projection of ( y_2 ) onto ( \varepsilon_1 ))
   • ( b_{22} = \sqrt{R[2,2] - b_{21}^2} ) (the norm of the residual after projection)
3. Third component:
   • ( y_3 = b_{31}\varepsilon_1 + b_{32}\varepsilon_2 + b_{33}\varepsilon_3 ), with ( E[\varepsilon_3\varepsilon_1] = E[\varepsilon_3\varepsilon_2] = 0 ).
   • ( b_{31} = \frac{R[3,1]}{b_{11}} )
   • ( b_{32} = \frac{R[3,2] - b_{21}b_{31}}{b_{22}} ) (subtract the projection onto ( \varepsilon_1 ) first, then project onto ( \varepsilon_2 ))
   • ( b_{33} = \sqrt{R[3,3] - b_{31}^2 - b_{32}^2} )

MATLAB Implementation #

Here's a function that implements this logic (works for any size covariance matrix, not just 3×3):

function B = gram_schmidt_from_cov(R)
    n = size(R, 1);
    B = zeros(n, n);
    
    % Handle first row/column
    B(1,1) = sqrt(R(1,1));
    
    % Iterate over each subsequent variable
    for i = 2:n
        % Compute off-diagonal entries (projections)
        for j = 1:i-1
            sum_proj = 0;
            % Subtract projections onto all previous orthogonal directions
            for k = 1:j-1
                sum_proj = sum_proj + B(i,k)*B(j,k);
            end
            B(i,j) = (R(i,j) - sum_proj)/B(j,j);
        end
        % Compute diagonal entry (residual norm)
        sum_sq = sum(B(i,1:i-1).^2);
        B(i,i) = sqrt(R(i,i) - sum_sq);
    end
end

Quick Verification #

Test it with a sample covariance matrix to make sure it works:

% Example 3x3 covariance matrix
R = [4, 2, 1;
     2, 5, 3;
     1, 3, 6];

B = gram_schmidt_from_cov(R);
% Check if B*B' equals R (should match up to numerical precision)
disp('Covariance matrix from B*B''');
disp(B*B');
disp('Original covariance matrix R');
disp(R);

Bonus Note #

If you just need the lower triangular matrix ( B ) and don't need to explicitly walk through Gram-Schmidt steps, MATLAB's built-in Cholesky decomposition does exactly this (and is more numerically stable):

B = chol(R, 'lower');

This gives the same result as our manual Gram-Schmidt implementation, since Cholesky decomposition is essentially a numerically optimized version of the Gram-Schmidt process for positive semi-definite matrices.

内容的提问来源于stack exchange，提问作者Resting Platypus