Eigen 求协方差矩阵

方差和协方差

方差：

方差是在概率论和统计方差衡量随机变量或一组数据时离散程度的度量。概率论中方差用来度量随机变量和其数学期望（即均值）之间的偏离程度。统计中的方差（样本方差）是每个样本值与全体样本值的平均数之差的平方值的平均数。在许多实际问题中，研究方差即偏离程度有着重要意义。-----百度词条

方差是衡量源数据和期望值相差的度量值。

公式：

s2=∑i=1n(Xi−X‾)2n−1s^2=\frac{\sum_{i=1}^{n}(X_i-\overline{X})^2}{n-1}s2=n−1∑i=1n(Xi−X)2

其中：

X‾=∑i=1nXin\overline{X}=\frac{\sum_{i=1}^nX_i}{n}X=n∑i=1nXi

协方差

可以看到方差是表示数据整体和期望也就是均值的偏差程度，方差越大表示数据越参差不齐，或者说波动越大，在控制中可以表述为噪声或干扰越大，其可信度越低。

方差只能表示数组自身的信息也就是一个随机变量的信息，如果需要知道两个随机变量之间的关系需要用协方差

协方差表示的是两个变量的总体的误差，这与只表示一个变量误差的方差不同。如果两个变量的变化趋势一致，也就是说如果其中一个大于自身的期望值，另外一个也大于自身的期望值，那么两个变量之间的协方差就是正值。如果两个变量的变化趋势相反，即其中一个大于自身的期望值，另外一个却小于自身的期望值，那么两个变量之间的协方差就是负值。—百度词条

通俗讲协方差表示两个随机变量的关联度和关联方向，即变化趋势是相同还是相反或者不相关

公式：

Cov(X,Y)=∑i=1n(Xi−X‾)(Yi−Y‾)n−1Cov(X,Y)=\frac{\sum_{i=1}^{n}(X_i-\overline{X})(Y_i-\overline{Y})}{n-1}Cov(X,Y)=n−1∑i=1n(Xi−X)(Yi−Y)

协方差矩阵

但是协方差只能用来表述一维数据，或者说包含变量包含一个维度的数据，实际现实中往往有多维变量需要处理，每维度都需要计算，这时就需要计算多个协方差，为了便于计算和管理，用协方差组成的矩阵来表示，即协方差矩阵。

样本组成矩阵，其中的数据按行排列与按列排列求出的协方差矩阵是不同的，这里默认数据是按行排列（我各人喜欢按列排列即用列向量表示一个样本）。每一列是一个样本或观测结果，那么每一行就是一个随机变量。

Xm×n=[a11a12⋯a1na21a22⋯a2n⋮⋮⋱⋮am1am2⋯amn]=[c1c2⋯cm]TX_{m\times n}=\begin{bmatrix}a_{11}\hspace{0.5cm}a_{12}\hspace{0.5cm}\cdots\hspace{0.5cm}a_{1n}\\a_{21}\hspace{0.5cm}a_{22}\hspace{0.5cm}\cdots\hspace{0.5cm}a_{2n}\\\vdots\hspace{0.5cm}\vdots\hspace{0.5cm}\ddots\hspace{0.5cm}\vdots\\a_{m1}\hspace{0.5cm}a_{m2}\hspace{0.5cm}\cdots\hspace{0.5cm}a_{mn} \end{bmatrix}=[c_1\hspace{0.5cm}c_2\hspace{0.5cm}\cdots\hspace{0.5cm}c_m]^{T}Xm×n=⎣⎢⎢⎢⎡a11a12⋯a1na21a22⋯a2n⋮⋮⋱⋮am1am2⋯amn⎦⎥⎥⎥⎤=[c1c2⋯cm]T

有n个数据，每个数据有m个随机变量（维度），那么协方差矩阵为：

∑=1n−1(cov(c1,c1)cov(c1,c2)⋯cov(c1,cm)cov(c2,c1)cov(c2,c2)⋯cov(c2,cm)⋮⋮⋱⋮cov(cm,c1)cov(cm,c2)⋯cov(cm,cm))\sum=\frac{1}{n-1}\left(\begin{matrix}cov(c_1,c_1) &cov(c_1,c_2)&\cdots&cov(c_1,c_m)\\cov(c_2,c_1)&cov(c_2,c_2)&\cdots&cov(c_2,c_m)\\\vdots&\vdots&\ddots&\vdots\\cov(c_m,c_1)&cov(c_m,c_2)&\cdots&cov(c_m,c_m)\end{matrix} \right)∑=n−11⎝⎜⎜⎜⎛cov(c1,c1)cov(c2,c1)⋮cov(cm,c1)cov(c1,c2)cov(c2,c2)⋮cov(cm,c2)⋯⋯⋱⋯cov(c1,cm)cov(c2,cm)⋮cov(cm,cm)⎠⎟⎟⎟⎞

例如slam中的定位使用卡尔曼滤波（现在基本不用了，目前流行使用非线性优化和图优化算法，但是思想是值得学习的）机器人状态由位置和速度来表示XT=[p,v]T=[x,y,z,vx,vy,vz]TX^T=[p,v]^{T}=[x,y,z,v_x,v_y,v_z]^{T}XT=[p,v]T=[x,y,z,vx,vy,vz]T ，如果他们各自都有三个维度。卡尔曼滤波中的一重要的参数就是Pk−P_k^-Pk− ，表示位置和速度的先验协方差。每个变量有三维度，合成的向量XTX^TXT有6个维度，所以协方差矩阵是一个6×66\times66×6的矩阵

可以表示为：Pk−=[∑pp∑pv∑pv∑vv]P_k^-=\begin{bmatrix}\sum{pp} & \sum{pv}\\\sum{pv} & \sum{vv}\end{bmatrix}Pk−=[∑pp∑pv∑pv∑vv]

协方差矩阵的计算

如果按照定义把每个协方差都计算出来然后组成矩阵太过麻烦了，好在利用矩阵计算可以大大简化。

把样本组成矩阵(假设有n个数据，每个数据有m个维度,XiX_iXi表示随机变量)

可以每行一个样本也可以每列一个样本
- 按列 Xm×n=[a11a12⋯a1na21a22⋯a2n⋮⋮⋱⋮am1am2⋯amn]m×n=[X1X2⋯Xm]TX_{m\times n}=\begin{bmatrix} {a_{11}}\hspace{0.5cm}a_{12}\hspace{0.5cm}\cdots\hspace{0.5cm}a_{1n}\\a_{21}\hspace{0.5cm}a_{22}\hspace{0.5cm}\cdots\hspace{0.5cm}a_{2n}\\\vdots\hspace{0.5cm}\vdots\hspace{0.5cm}\ddots\hspace{0.5cm}\vdots\\a_{m1}\hspace{0.5cm}a_{m2}\hspace{0.5cm}\cdots\hspace{0.5cm}a_{mn}\end{bmatrix}_{m\times n}=[X_1\space X_2\space\cdots\space X_m]^TXm×n=⎣⎢⎢⎢⎡a11a12⋯a1na21a22⋯a2n⋮⋮⋱⋮am1am2⋯amn⎦⎥⎥⎥⎤m×n=[X1 X2 ⋯ Xm]T
- 按行 Xn×m=[a11a12⋯a1ma21a22⋯a2m⋮⋮⋱⋮an1an2⋯anm]n×m=[X1X2⋯Xm]X_{n\times m}=\begin{bmatrix}a_{11}\hspace{0.5cm}a_{12}\hspace{0.5cm}\cdots\hspace{0.5cm}a_{1m}\\a_{21}\hspace{0.5cm}a_{22}\hspace{0.5cm}\cdots\hspace{0.5cm}a_{2m}\\\vdots\hspace{0.5cm}\vdots\hspace{0.5cm}\ddots\hspace{0.5cm}\vdots\\a_{n1}\hspace{0.5cm}a_{n2}\hspace{0.5cm}\cdots\hspace{0.5cm}a_{nm}\end{bmatrix}_{n\times m}=[X_1\space X_2\space\cdots\space X_m]Xn×m=⎣⎢⎢⎢⎡a11a12⋯a1ma21a22⋯a2m⋮⋮⋱⋮an1an2⋯anm⎦⎥⎥⎥⎤n×m=[X1 X2 ⋯ Xm]
计算每个维度的均值
- 按列 Mean=[M1M2⋯Mm]TMean=[M_1\hspace{0.5cm}M_2\hspace{0.5cm}\cdots\hspace{0.5cm}M_m]^TMean=[M1M2⋯Mm]T
- 按行 Mean=[M1M2⋯Mm]Mean=[M_1\hspace{0.5cm}M_2\hspace{0.5cm}\cdots\hspace{0.5cm}M_m]Mean=[M1M2⋯Mm]
每个样本减区各个维度的均值，即求(Xi−X‾)(X_i-\overline X )(Xi−X)
- 按列temp=[a11−M1a12−M1⋯a1n−M1a21−M2a22−M2⋯a2n−M2⋮⋮⋱⋮am1−Mmam2−Mm⋯amn−Mm]m×ntemp=\begin{bmatrix}a_{11}-M_1\hspace{0.5cm}a_{12}-M_1\hspace{0.5cm}\cdots\hspace{0.5cm}a_{1n}-M_1\\a_{21}-M_2\hspace{0.5cm}a_{22}-M_2\hspace{0.5cm}\cdots\hspace{0.5cm}a_{2n}-M_2\\\vdots\hspace{0.5cm}\vdots\hspace{0.5cm}\ddots\hspace{0.5cm}\vdots\\a_{m1}-M_m\hspace{0.5cm}a_{m2}-M_m\hspace{0.5cm}\cdots\hspace{0.5cm}a_{mn}-M_m\end{bmatrix}_{m\times n}temp=⎣⎢⎢⎢⎡a11−M1a12−M1⋯a1n−M1a21−M2a22−M2⋯a2n−M2⋮⋮⋱⋮am1−Mmam2−Mm⋯amn−Mm⎦⎥⎥⎥⎤m×n
- 按行temp=[a11−M1a12−M2⋯a1m−Mma21−M1a22−M2⋯a2m−Mm⋮⋮⋱⋮an1−M1an2−M2⋯anm−Mm]n×mtemp=\begin{bmatrix}a_{11}-M_1\hspace{0.5cm}a_{12}-M_2\hspace{0.5cm}\cdots\hspace{0.5cm}a_{1m}-M_m\\a_{21}-M_1\hspace{0.5cm}a_{22}-M_2\hspace{0.5cm}\cdots\hspace{0.5cm}a_{2m}-M_m\\\vdots\hspace{0.5cm}\vdots\hspace{0.5cm}\ddots\hspace{0.5cm}\vdots\\a_{n1}-M_1\hspace{0.5cm}a_{n2}-M_2\hspace{0.5cm}\cdots\hspace{0.5cm}a_{nm}-M_m\end{bmatrix}_{n \times m}temp=⎣⎢⎢⎢⎡a11−M1a12−M2⋯a1m−Mma21−M1a22−M2⋯a2m−Mm⋮⋮⋱⋮an1−M1an2−M2⋯anm−Mm⎦⎥⎥⎥⎤n×m
求1n−1∑i,jm(Xi−X‾i)(Xj−X‾j)\frac{1}{n-1}\sum_{i,j}^{m}(X_i-\overline X_i)(X_j-\overline X_j)n−11∑i,jm(Xi−Xi)(Xj−Xj)两两随机变量求协方差,组成矩阵即协方差矩阵
- 按列∑=1n−1temp∗tempT\sum=\frac{1}{n-1}temp*temp^T∑=n−11temp∗tempT
- 按行∑=1n−1tempT∗temp\sum=\frac{1}{n-1}temp^T*temp∑=n−11tempT∗temp
注意n不能为1

使用Eigen库计算

#include <iostream>
#include<Dense>
using namespace Eigen;
using namespace std;
int main()
{// //  （0.43，1.1，0,0.32,0.21,0)// （0.4，1.28，0,0.34,0.17,0)// （0.39，1.2，0,0.31,0.25,0)// （0.45，1.4，0,0.22,0.13,0)Eigen::MatrixXf mat_1(6, 4);//按列MatrixXf mat_2(4, 6);//按行MatrixXf cov_1 = MatrixXf::Zero(6,6);MatrixXf cov_2 = MatrixXf::Zero(6,6);mat_2<< 0.43, 1.1, 0, 0.32, 0.21, 0,0.4, 1.28, 0, 0.34, 0.17, 0,0.39, 1.2, 0, 0.31, 0.25, 0,0.45, 1.4, 0, 0.22, 0.13, 0;mat_1 = mat_2.transpose();//1 求各个维度均值MatrixXf mean_1 = mat_1.rowwise().mean();MatrixXf mean_2 = mat_2.colwise().mean();//2 每个样本减去均值for (size_t i = 0; i < mat_1.cols(); i++){cov_1.col(i) = mat_1.col(i) - mean_1;}for (size_t i = 0; i < mat_2.rows(); i++){cov_2.row(i) = mat_2.row(i) - mean_2;}//3 计算cov_1 = cov_1 * cov_1.transpose()/(cov_1.cols()-1);cov_2 = cov_2.transpose() * cov_2 / (cov_2.rows() - 1);cout << cov_1 << endl<<endl << cov_2;}

运行结果

可以看到是按行计算和按列计算一样的