posts - 217, comments - 61, trackbacks - 0, articles - 0
   :: 首页 :: 新随笔 :: 联系 :: 聚合  :: 管理

协方差矩阵

Posted on 2011-07-05 18:17 魔のkyo 阅读(2787) 评论(1)  编辑 收藏 引用 所属分类: Math
首先一个维随机变量的均值记为E[X]
两个随机变量的协方差定义为cov(X,Y) = E[(X-E[X])(Y-E[Y])]
两个随机向量的协方差矩阵由随机向量分量间的协方差定义而成。
协方差矩阵元素Aij = cov(Xi,Yj)
协方差矩阵的阶=随机向量的维数,而与随机向量的个数无关。
协方差矩阵是对称矩阵。
其主对角线元素为点的各个分量的方差,可以反映点集在各个轴向的离散程度。
其他元素反映了不同分量的相关性,当两个分量距各自均值发生同向偏离时,协方差(即相应的协方差矩阵元素)趋向正值,当两个分量距各自均值发生异向偏离时,协方差趋向负值。
通过计算协方差矩阵的特征值和特征向量,可以反映出点集离散的主要方向。若特征向量为最大特征向量,则特征值即为沿着该轴的顶点数据具有的最大方差值;若特征向量为最小特征向量,则特征值即为沿着该轴的顶点数据具有的最小方差值。
特征值和特征向量的计算可以使用Jacobi方法,可以参考《实时碰撞检测算法技术》Page63~65,《Real-Time Collision Dectection》Page95~97

实用意义:
估计一个模型(一个点集)的最小包围球。
估计一个模型(一个点集)的最小OBB。
估计一个2D图形(二维点集)的最小面积包围矩形。
求2D点集的拟合直线。
求3D点集的拟合平面。

//3维向量协方差矩阵
void CovarianceMatrix(__out Matrix3f &covarianceMatrix, __in const Vector3f pt[], __in int numPts)
{
    
float invNum = 1.0f / (float)numPts;

    Vector3f sum(
0.0f);
    
for(int i=0;i<numPts;i++)
    {
        sum 
+= pt[i];
    }
    Vector3f u 
= sum * invNum;
    
for(int i=0;i<3;i++)
    {
        
for(int j=i;j<3;j++)
        {
            covarianceMatrix[i][j] 
= 0.0f;
            
for(int k=0;k<numPts;k++)
            {
                covarianceMatrix[i][j] 
+= (pt[k][i] - u[i]) * (pt[k][j] - u[j]);
            }
            covarianceMatrix[i][j] 
*= invNum;
        }
    }
    
for(int i=0;i<3;i++)
    {
        
for(int j=0;j<i;j++)
        {
            covarianceMatrix[i][j] 
= covarianceMatrix[j][i];            
        }
    }
}

// 2-by-2 Symmetric Schur decomposition. Given an n-by-n symmetric matrix
// and indices p, q such that 1 <= p < q <= n, computes a sine-cosine pair
// (s, c) that will serve to form a Jacobi rotation matrix.
//
// See Golub, Van Loan, Matrix Computations, 3rd ed, p428
void SymSchur2(__in const Matrix3f &a, __in int p, __in int q, __out float* c, __out float* s)
{
    
if (abs(a[p][q]) > 0.0001f)
    {
        
float r = (a[q][q] - a[p][p]) / (2.0f * a[p][q]);
        
float t;
        
if (r >= 0.0f)
            t 
= 1.0f / (r + sqrt(1.0f + r*r));
        
else
            t 
= -1.0f / (-+ sqrt(1.0f + r*r));
        
*= 1.0f / sqrt(1.0f + t*t);
        
*= t * *c;
    }
    
else
    {
        
*= 1.0f;
        
*= 0.0f;
    }
}

// Computes the eigenvectors and eigenvalues of the symmetric matrix A using
// the classic Jacobi method of iteratively updating A as A = J∧T * A * J,
// where J = J(p, q, theta) is the Jacobi rotation matrix.
//
// On exit, v will contain the eigenvectors, and the diagonal elements
// of a are the corresponding eigenvalues.
//
// See Golub, Van Loan, Matrix Computations, 3rd ed, p428
void Jacobi(__inout Matrix3f &a, __out Matrix3f &v)
{
    v.SetIdentity();

    
float prevoff;
    
// Repeat for some maximum number of iterations
    const int MAX_ITERATIONS = 50;
    
for (int k = 0; k < MAX_ITERATIONS; k++)
    {
        
// Find largest off-diagonal absolute element a[p][q]
        int p = 0, q = 1;
        
for (int i = 0; i < 3; i++)
        {
            
for (int j = 0; j < 3; j++)
            {
                
if (i != j && abs(a[i][j]) > abs(a[p][q]))
                {
                    p 
= i;
                    q 
= j;
                }
            }
        }
        
// Compute the Jacobi rotation matrix J(p, q, theta)
        
// (This code can be optimized for the three different cases of rotation)
        float c,s;
        SymSchur2(a, p, q, 
&c, &s);
        Matrix3f J;
        J.SetIdentity();
        J[p][p] 
= c; J[p][q] = s;
        J[q][p] 
= -s; J[q][q] = c;
        
// Cumulate rotations into what will contain the eigenvectors
        v = v * J;
        
// Make ’a’ more diagonal, until just eigenvalues remain on diagonal
        a = (J.Transpose() * a) * J;
        
// Compute "norm" of off-diagonal elements
        float off = 0.0f;
        
for (int i = 0; i < 3; i++)
        {
            
for (int j = 0; j < 3; j++)
            {
                
if (i != j)
                {
                    off 
+= a[i][j] * a[i][j];
                }
            }
        }
        
/* off = sqrt(off); not needed for norm comparison */
        
// Stop when norm no longer decreasing
        if (k > 2 && off >= prevoff)
            
return;
        prevoff 
= off;
    }
}
需要注意的是,
Jacobi(__inout Matrix3f &a, __out Matrix3f &v);
参数a要求是对称矩阵。
函数返回时a是一个与传入矩阵相似的对角矩阵,v中的特征向量是按列存储的。
三个特征向量分别是
Vector3f(v[0][0],v[1][0],v[2][0])
Vector3f(v[0][1],v[1][1],v[2][1])
Vector3f(v[0][2],v[1][2],v[2][2])

Feedback

# re: 协方差矩阵  回复  更多评论   

2013-12-18 18:19 by ee
很好的
只有注册用户登录后才能发表评论。