从协方差矩阵的估算领会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)所示。
其中，，mi表示第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”循环。然后随机变量Xi的观测值使用矩阵索引表示为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维随机向量的协方差矩阵(Xi是随机变量)
于是结合矩阵运算的规则，得到基于矩阵计算随机向量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函数 |
|
1 |
0.001187 |
0.001277 |
0.000582 |
0.000279 |
|
2 |
0.001197 |
0.001305 |
0.000618 |
0.000868 |
|
3 |
0.001414 |
0.001264 |
0.000574 |
0.000282 |
|
平均值 |
0.001266 |
0.001282 |
0.0005913 |
0.0004763 |
由表可知，基于“逐点”、向量和矩阵三个层次编程实现对协方差矩阵的估算，其中矩阵方法用时最小，约为其他两种方法的一半，系统自带的cov函数用时最短(应该是使用更底层语言编写、并且已经编译，执行效率更高)，约是矩阵法的一半，因此矩阵法耗时最少。另一方面由“逐点”、向量到矩阵法，程序越来越简洁。
3. 结论
由以上分析，可看出由“逐点”、向量到矩阵法，程序越来越简洁。执行效率矩阵法最高，耗时约是“逐点法”的一半。
对于多维(多类、多组)同类型数据处理，建议使用矩阵编程，总结了至少有3个优点：
(1)省时；
(2)程序少(敲代码少)，于是容易修改、纠错…
(3)便于数据组织。尤其对于多维、多组同类型数据处理，要对数据的不同维分别处理，这样用矩阵存储、组织临时数据将非常方便。如彩色图像处理中对多个彩色图像文件的各个颜色分量处理，这时结合“for”循环很方便。但是以矩阵方式组织数据，要求数据结构(维数/尺寸)相同。
4. 参考资料