Loading [MathJax]/jax/output/CommonHTML/fonts/TeX/AMS-Regular.js
2022-09-12 15:12阅读: 3917评论: 0推荐: 0

多元高斯分布的协方差矩阵

看《概率机器人》时总是遇到这个协方差矩阵,但对这家伙一直没有个直观认识,因此集中学一下,同时顺带捡起来一点点线代知识。

定义

设待讨论的随机变量XiXi构成随机向量X=[X1,X2,,Xn]TX=[X1,X2,,Xn]T,则协方差矩阵定义为一个n×nn×n矩阵ΣΣ,其第iijj列的元素为

Σij=Cov(Xi,Xj)Σij=Cov(Xi,Xj)

现实中的协方差矩阵是基于对随机变量的采样数据进行计算而得到的,假设每一次对nn个随机变量进行采样,进而构成随机向量,并且共进行了kk次采样,第i次的结果记为xi=[xi1,xi2,,xin]Txi=[xi1,xi2,,xin]T,则协方差矩阵计算如下

Σ=1kni=1(xiˉx)(xiˉx)TΣ=1kni=1(xi¯x)(xi¯x)T

从中还可以看出,如果对随机向量进行线性变换TT,则协方差矩阵会变为TΣTTTΣTT(注意是后边那个做的转置)。

性质

对称

由定义直接可以得到。

线代拾遗

  • 实对称矩阵的特征向量相互正交。
    假设实对称矩阵AA满足Au1=λ1u1Au1=λ1u1Au2=λ2u2Au2=λ2u2,则有

    uT2Au1=λ1uT2u1=(uT1Au2)T=(λ1uT2u1)T=λ1uT2u1uT2Au1=λ1uT2u1=(uT1Au2)T=(λ1uT2u1)T=λ1uT2u1

    因为λ1λ2λ1λ2,所以必有uT2u1=0uT2u1=0,即两特征向量正交。

  • 各列向量为单位向量且互相正交的矩阵(也称为正交矩阵)的逆就是其本身的转置。

    假设AA的各列正交且为单位向量,则计算一下即可发现ATA=IATA=I。后边解释协方差矩阵可以拆成一个矩阵和其转置相乘时会用到这一点。

半正定

根据定义,协方差矩阵可由下式计算

Σ=E[(xμ)(xμ)T]Σ=E[(xμ)(xμ)T]

故对于任意与xx具有相同个数元素的向量uu

uTΣu=E[(uTxμ)(xμ)Tu]=E{[uT(xμ)]2}0uTΣu=E[(uTxμ)(xμ)Tu]=E{[uT(xμ)]2}0

即可得多元高斯分布的协方差矩阵具有半正定的性质。

在浏览有关协方差矩阵的材料中,经常可以看到一些假定协方差矩阵正定(同时也会满秩)的操作,这其实不是一定的(尽管大多数情况下还是会正定),参考[1]。

现在展开说一说,对于nn维非零随机向量uu,做如下操作

uTΣu=1kni=1uT(xiˉx)(xiˉx)Tu=1kni=1[(xiˉx)Tu]T(xiˉx)Tu=1kni=1(zTiu)2uTΣu=1kni=1uT(xi¯x)(xi¯x)Tu=1kni=1[(xi¯x)Tu]T(xi¯x)Tu=1kni=1(zTiu)2

其中zi=xiˉxzi=xi¯x,故uTΣuuTΣu为0的充要条件为zTiu=0zTiu=0

现假设uTΣu=0uTΣu=0,且span(z1,z2,,zk)=Rn,则u可以表示为u=α1z1+α2z2++αkzk,继而可得到

uTu=α1zT1u+α2zT2u++αkzTku=0

这说明u=0,与前提条件矛盾,故span(z1,z2,,zk)=Rn时基于采样数据计算得到的协方差矩阵即会正定,也即rank(z1,z2,,zk)=n时即正定。当采样数据足够多时,这一条件是不难达成的,因此大多数实际情况中协方差矩阵应是正定的。

线代拾遗

正定矩阵的特征值均为正数。

假定正定矩阵A满足Au=λu,则uTAu=λuTu>0,因为uTu>0,故λ>0

可拆成一个矩阵和其转置相乘

假设采样数据足够多,则协方差矩阵正定,正定则也就满秩,因此可对其进行相似对角化:

Σ=UΛU1=UΛ12Λ12UT=(UΛ12)(UΛ12)T=TTT

有了这一性质,一般的多维高斯分布都可以看作是标准多为高斯分布的线性变换再加上一个均值向量,其中线性变换就是上边的T,均值向量就是μ

协方差矩阵特征向量、特征值的几何意义

意义1:投影到一维后具有极值标准差

先说结论:随机向量在新协方差矩阵特征向量方向上的投影(这是一个新的一维随机变量),相比于在其他方向上的投影,取到了方差的极值,并且该极值就是该特征向量对应的特征值。

注意:如果把各个方向上一倍标准差上取点并连线,得到的并不是一个椭圆(因为《概率机器人》习题中的描述我一度误以为这应该是个椭圆)。

然后从理论上证明一下:设XN(μ,Σ),则采样数据在单位ω方向上的投影的方差为:

σ2=1nki[ωT(xiμ)]2=1nki[ωT(xiμ)][ωT(xiμ)]T=1nkiωT(xiμ)2ω=ωT[1kni=1(xiˉx)(xiˉx)T]ω=ωTΣω

求在ωTω=1的条件上令σ2取极值的ω,可以通过拉格朗日乘数法构造一个目标函数:

L(ω,λ)=ωTΣω+λ(1ωTω)

对其求偏导并令偏导为0,可以得到取极值的条件:

L(ω,λ)ω=2ωTΣ+2λωT=0Σω=λω

即有ωΣ的特征向量。将该结论带回上式,可以得到此时还满足σ2=λ

最后在二维下看一个例子,尝试修改不同的transform变量,会发现表示特征向量的箭头总是指向数据最分散或最集中的方向,且长的向量总是指向数据分散的方向,短的向量总是指向数据集中的方向。

import numpy as np
import matplotlib.pyplot as plt
# 生成白噪声
x_white = np.random.normal(0, 1.0, 1000)
y_white = np.random.normal(0, 1.0, 1000)
plt.scatter(x_white, y_white, 2)
plt.xlim((-5, 5))
plt.ylim((-5, 5))
plt.grid()
plt.gca().set_aspect('equal')
plt.show()


png

def linear_transform(T, x_, y_):
x = T[0, 0] * x_ + T[0, 1] * y_
y = T[1, 0] * x_ + T[1, 1] * y_
return x, y
# 对白噪声进行线性变换得到一般的二维高斯分布
transform = np.array([[6, 3], [7, 1]])
x, y = linear_transform(transform, x_white, y_white)
plt.scatter(x, y, 2)
# 画新协方差矩阵transform * transform^T的特征向量
eig, eig_vecs = np.linalg.eig(np.matmul(transform, transform.T))
plt.arrow(0, 0, eig_vecs[0, 0] * np.sqrt(eig[0]), eig_vecs[1, 0] * np.sqrt(eig[0]), width=0.1, head_width=0.5, color='red')
plt.arrow(0, 0, eig_vecs[0, 1] * np.sqrt(eig[1]), eig_vecs[1, 1] * np.sqrt(eig[1]), width=0.1, head_width=0.5, color='purple')
plt.grid()
plt.gca().set_aspect('equal')
plt.show()

png

意义2:误差椭圆长短轴的方向与大小

多元高斯分布的概率密度函数的等高线才是误差椭圆(error ellipse / uncertainty ellipse)。

先以二维为例,高维应该可以推而广之。观察多为高斯分布的定义式,等高线应该满足:

(xμ)TΣ1(xμ)=Const

为了方便起见,不如就取Const=1,且设z=xμ。另外,设协方差矩阵可以相似对角化到Σ=UΛU1,其中U是正交矩阵,则协方差矩阵的逆矩阵为Σ1=UΛ1U1=UΛ1UT。于是上式化简为:

zTUΛ1UTz=1zT[u1u2][λ11λ12][uT1uT2]z=1(uT1xλ1)2+(uT2xλ2)2=1

因为u1u2是相互垂直的单位向量(正交矩阵的性质),所以uT1xuT2x就是向量x在这两个特征向量方向上的投影,也即[uT1uT2]z是以这两个特征向量为新的基对x做了一次线性变换。在这一组新的基下,椭圆方程即是标准的x2/a2+y2/b2=1的形式,因此也可以看出长短轴就是λ1λ2,长短轴的方向就是特征向量的方向。

参考

神奇又好玩的协方差矩阵

如何直观地理解「协方差矩阵」?

A geometric interpretation of the covariance matrix

PCA(主成分分析)原理推导

Drawing Ellipse from eigenvalue-eigenvector

本文作者:Harold_Lu

本文链接:https://www.cnblogs.com/harold-lu/p/16686272.html

版权声明:本作品采用知识共享署名-非商业性使用-禁止演绎 2.5 中国大陆许可协议进行许可。

posted @   Harold_Lu  阅读(3917)  评论(0编辑  收藏  举报
点击右上角即可分享
微信分享提示
评论
收藏
关注
推荐
深色
回顶
收起