从协方差矩阵的估算领会MATLAB矩阵编程思维

摘要：本文以使用一个随机向量的样本数据来估算(我们算的都是估算值)两个随机向量的协方差矩阵为例，分别从基于“逐点”、向量和矩阵三个层次编程实现对协方差矩阵的估算，展示MATLAB矩阵编程的思维和优点。

欢迎指正错误，欢迎交流。转载引用请注明出处(http://www.cnblogs.com/xikeguanyu/articles/5516782.html)。

1. 前言

MATLAB，这个从大学开始接触的软件，随着使用次数的增加，越来越感慨它的便利性，尤其在算法仿真验证上，使用它来编程进行验证，真的很方便。有点打广告的感觉？但真不是打广告的^_^。虽然说在控制硬件方面如实时采集信号等工程应用中使用MATLAB进行编程的比较少，但是其在数据处理和分析、建模上广为应用，再者它很方便有效且容易入门，因此学习用来做仿真、验证是很有必要的。

算法的验证，往往会涉及到数据，数据即证据，算法的可行性需要证据来验证！算法往往是一些公式、步骤，可能是连续的符号，而计算机处理数据必须是离散时间的，因此编程来验证时，自然涉及到对于数据的组织即数据结构的问题。对于很多问题，我们或许已经感受到，使用矩阵或向量来描述，分析处理会更加简洁方便。

MATLAB是以矩阵数值计算为特长的软件，在矩阵运算上真的很有效。有见到网友说，以矩阵思维编程，在MATLAB上可以达到高的运行速度，然后推荐使用矩阵编程而不是“for 循环”。

有见到某程序语言书籍上说，“程序 = 数据结构 + 算法”，学习中一些算法的验证，使用MATLAB进行编程来验证，作者再次体会到了这一说法。

本文通过举例说明MATLAB的矩阵编程思维及编程中的数据结构设计，分别从基于“逐点”、向量和矩阵三个层次编程实现对协方差矩阵的估算，向大家展示MATLAB矩阵编程的思维和优点。希望对大家有所收获和启示，若能让你喜欢上MATLAB，甚至激发其他兴趣，那将是对本文最大的鼓励和本文的一种重要价值体现。欢迎指正错误，欢迎交流。转载引用请注明出处。

1. 问题

给定M维随机向量(或M元随机变量)的观测值X，X为M*N的矩阵，即有N个观测值(每个观测值是一个M维的列向量，再具体，该列向量中的每个元素分别对应M维随机向量中的随机变量)。使用MATLAB估算M维随机向量(或M元随机变量)的协方差矩阵Cov。

2. 原理与方法

2.1 协方差的概念和意义

随机变量X1、X2间的协方差记为Cov(X1, X2)定义如下：

Cov(X, Y) = E[(X- E(X)(Y – E(Y)))]                  (1)

其中E(*)表示数学期望。

协方差来自于概率论，反映两个随机变量之间的相关性(线性相关)。(目前个人理解为正负)注意，从式(1)中，可以认为，只要X与Y同时比各自的均值大，即X增大时Y也增大，即相对于各自的均值的变化方向相同，则X、Y的协方差为正值，反之相反则取值为负。两个随机变量的协方差的取值大小与各自的观测值有关，因此其协方差的绝对值的大小并不能反映出其相关程度的高低。考虑到这一点，对协方差做某种归一化处理(也就是“相对”)，便可以由归一化的大小反映出相关程度，也就有了后面的相关系数。

(相关系数是以比值形式出现的，即“相对”的考查方式，类似与相对/比较的应用很多，比如生活中的比较，商品价格比高低；两辆同方向行驶的车的速度，只有相对速度大了才会更快地拉开距离；误差分析中的相对误差，PCM系统的非均匀量化就是采用了相对误差从而保护了小信号......)

2.2 协方差的计算

公式(1)只是数学公式，使用计算机来估算协方差，需要将数学公式翻译为程序语言，或者至少是“离散时间的数学形式”。

设随机变量X、Y的观测值分别为N个，即X={xi|i=1,2,…,N}，Y={yi|i=1,2,…,N}。

则对于X、Y间的协方差Cov估算公式为：

其中mx和my分别表示随机变量X、Y的N个观测值的平均值。即

2.3 多维随机变量(随机向量)间的协方差——协方差矩阵的计算

设X是M维随机列向量，即X = [X1,X2,…,XM]T，则随机向量X的协方差矩阵Cov的估算公式如式(4)所示。

其中，，m_i表示第i个随机变量的N个观测值的平均值。计算完后将得到一个M*M的对称矩阵。

2.4 使用MATLAB估算协方差矩阵

编写程序实现对算法的验证，需要将数学公式转化为时间离散化(至于由于计算机的存储器有限，需要幅值离散化，很多情况下尤其在MATLAB的64位double类型数据下基本可以不用考虑)的公式。式(2)、(3)和(4)就是已经离散时间化的公式，然后就可以编程了。以下分别从基于“逐点”、向量和矩阵三个层次编程实现对协方差矩阵的估算，向大家展示MATLAB矩阵编程的思维和优点。

计算M(M≥2)维随机向量的协方差公式，其观测值存储于一个M*N的矩阵X中，即共有N个(随机向量的)观测值(有M个随机变量的观测值，每个随机变量有N个观测值)。计算结果存储在一个M*M的方阵Cov中。

2.4.1基于逐点计算

(1). 计算公式

从式(4)中我们知道，协方差矩阵式一个M*M的对称方阵，因此基于逐点计算，需要两个“for”循环。然后随机变量X_i的观测值使用矩阵索引表示为X(i, k),k=1,2,…,N，即有N个观测值。因此计算公式如式(5)所示。

其中mi或mj如式(6)所示。

用MATLAB语言表示为。

(2). 计算程序

MATLAB计算程序如下：

——————————————————————————————

 1 Cov = zeros(M, M);   % 存储矩阵初始化
 2 for i = 1:M
 3     mi = mean(X(i, :));  % 计算mxi
 4     for j = 1:M
 5         mj = mean(X(j, :));  % 计算mxj
 6         sumt = 0;
 7         for k = 1:N           
 8             xi = X(i, k);  % 取到xi点      
 9             xj = X(j, k);  % 取到xj点          
10             sumt = sumt + (xi - mi)*(xj - mj);  % 完成累加计算  
11             % sumt = sumt + (X(i, k) - mi)*(X(i, k) - mj);  % 建议以上3行直接换为此句       
12         end
13         Cov(i, j) = sumt/(N-1);  % sumt/(N-1)完成一个值的计算 
14     end
15 end

——————————————————————————————

简要分析：从上面程序看，最大的特点就是使用了3个“for”循环！但循环中间部分即循环体却只有一句，效率不高。

2.4.2基于向量计算

(1). 计算公式

基于向量来计算，这里我很有必要深入叙述一下。再次提醒，协方差反映的是两个随机变量之间的相关性，即对于随机向量而言，反映的是向量各维的数据间的相关性。因此，我们为了计算，取的向量不是随机向量的观测值，而是该随机向量的某一维的所有观测值。即取X(i, :),i=1,2,…,M，是一个N维的行向量。相应的，求均值也是对某一维的N个观测值求平均，而不是按列即对向量的观测值(M维列向量)求平均。因此计算公式如式(7)所示。

细心的同学也可以发现从式(5)可以直接推到式(7)。

其中mi或mj同式(6)。<a, b>表示向量a,b之间的内积。注意X(i, :)是N维行向量，X(i, :) – mi等价于X(i, :) – mi*ones(1, N)，即X(i, :)的每个元素减去mi，这在MATLAB中是允许的。

(2). 计算程序

MATLAB计算程序如下：

——————————————————————————————

1 Cov = zeros(M, M);   % 存储矩阵初始化
2 for i = 1:M
3     mi = mean(X(i, :));  % 计算mxi
4     for j = 1:M
5         mj = mean(X(j, :));  % 计算mxj
6         Cov(i, j) = (X(i, :) - mi) * (X(i, :) - mj)\' ;  % 完成一个值的计算
7     end
8 end

——————————————————————————————

简要分析：从上面程序看，仅使用了2个“for”循环！循环中间部分即循环体虽一句，但执行相当于一个“for”循环，效率大幅提高。程序长度/行数也减小，变得简洁。

2.4.3基于矩阵计算

(1). 计算公式

由于协方差矩阵是对称矩阵(随机变量X和Y间的协方差Cov(X, Y)与顺序无关，即Cov(X, Y) = Cov(Y, X))，考查随机向量各维间的关系，可以列一个M*M的表格来表示要考虑到的协方差Cov(Xi, Xj)，如表1所示。

表1 M维随机向量的协方差矩阵(X_i是随机变量)

于是结合矩阵运算的规则，得到基于矩阵计算随机向量X的协方差矩阵的公式式(8)所示：

其中有

repmat和mean均是MATLAB计算函数。repmat(A, m, n)表示把A重复为m*n的矩阵，每个元素都是A。

(2). 计算程序

根据式(9)便可方便地得到MATLAB计算程序如下：

——————————————————————————————

1 Xm = repmat(mean(X, 2), 1, N);
2 % Cov = (X - repmat(mean(X, 2), 1, N)) * (X - repmat(mean(X, 2), 1, N))\'/(N-1); % 略费时间
3 Cov = (X - Xm) * (X - Xm)\'/(N-1);  % 得到协方差矩阵

——————————————————————————————

简要分析：从上面程序看，程序不再有“for”循环！ 核心语句只有2行！程序非常简洁。

2.4.4时间消耗分析

为比较三种算法的执行效率，这里使用了时间测量。测量时间使用tic/toc组合。测试条件为50个4维列向量，即X为4*50的矩阵，将三种算法放在一个程序文件中运行。还与MATLAB自带协方差计算函数cov进行了比较，cov在这里使用需要将X进行转置(详见help cov)。

时间消耗如表2所示。注意，表中逐点、向量栏时间包括Cov矩阵的初始化用时。

表2 自写三种计算方式和系统自带函数计算协方差矩阵时间消耗对比(单位：s[秒])

由表可知，基于“逐点”、向量和矩阵三个层次编程实现对协方差矩阵的估算，其中矩阵方法用时最小，约为其他两种方法的一半，系统自带的cov函数用时最短(应该是使用更底层语言编写、并且已经编译，执行效率更高)，约是矩阵法的一半，因此矩阵法耗时最少。另一方面由“逐点”、向量到矩阵法，程序越来越简洁。

3. 结论

由以上分析，可看出由“逐点”、向量到矩阵法，程序越来越简洁。执行效率矩阵法最高，耗时约是“逐点法”的一半。

对于多维(多类、多组)同类型数据处理，建议使用矩阵编程，总结了至少有3个优点：

(1)省时；

(2)程序少(敲代码少)，于是容易修改、纠错…

(3)便于数据组织。尤其对于多维、多组同类型数据处理，要对数据的不同维分别处理，这样用矩阵存储、组织临时数据将非常方便。如彩色图像处理中对多个彩色图像文件的各个颜色分量处理，这时结合“for”循环很方便。但是以矩阵方式组织数据，要求数据结构(维数/尺寸)相同。

4. 参考资料

1、http://blog.csdn.net/wzwind/article/details/50624398

2、http://www.zhihu.com/question/20852004