最小二乘法及其应用

1、第三章第三章 最小二乘法及其应用最小二乘法及其应用 最小二乘法是求解最优化问题的一种有效而方便的方法。信号处理中有许多问题可归结为最优化问题，因此最小二乘法是信号处理的重要工具之一。 希尔伯特空间中线性逼近问题的求解方法称为最小二乘法。通常它有三种不同的表现形式：投影法、求导法和配方法。下面来分别说明。3-1 最小二乘法的三种形式最小二乘法的三种形式 设X为希尔伯特空间， 为X中一组归一化正交元素，x为X中的某一元素。在子空间 中求一元素m。使得12 ,e e 12 ,Mspan e e0|min |mMxmxm（3-1-1）由于M中元素可表为 的线性组合，问题转化成为求 ，使得12,e e

2、12,a a 1|=min kkkxa e（3-1-2） 第二章中的投影定理指出了最优系数 应满足12,a a 1 ,1,2,kkmkxa eem（3-1-3） 由此即得 。也就是说，当且仅当 取为x关于归一化正交系 的傅立叶系数时，式(3-1-2)成立。kmk=1( ,)=(a, )=amkmx eee ka12 ,e e k( ,)ckkax e（3-1-4） 这种求解方法称为投影法，它是最小二乘法的第一种表现形式。第二种方法是求导法，仍以上面的问题为例来说明。记泛函2121(,) | kkkf a axa e为了能用求导法求此泛函的极小值，将它表为12112211(,)(, ) |2 k

3、kmmkkkkkkkfaaxa exa exa ca其中 。于是最优的 应满足即( ,)kkcx e12,a a =0 m=1,2,mfa20, m=1,2mmmmcaac或 ， 下面再用第三种方法即配方法来求解：2212112222111122211(,)|2 |2 |() min, 1,2 kkkkkkkkkkkkkkkkkkkkkf aaxa caxcca caxcacack（3-1-5） 以上三种方法都称为最小二乘法。在实际应用中，他们各有各的优势和缺陷，我们并不能通过简单的比较来说明他们谁优谁劣，因为衡量一种方法好坏的标准是多方面的。因此，在不同的场合根据不同的需要和可能，灵活选择和

4、使用合适的方法，是掌握最小二乘法的关键。 利用令导数等于零来求函数的极值是一种方便的方法。但是对于多元函数，有时由于变元太多而使表达式相当繁复，为此，本节介绍用向量-矩阵的形式来简化求导过程。 下面举例个例子来具体说明。例3-2-1 求矛盾方程组Ax=b的最小二乘 解（可参阅第二章的相关例题）3-2 向量向量-矩阵求导及配方法矩阵求导及配方法 解：解：求Ax=b的最小二乘解就是求 的极小点。由于2( ) |g xAxb( ) () ()2TTTTTg xAx bAx bx A Axb Ax b b下面先给出两个需要用到的向量求导公式：()()TTxxb xx b （3-2-1）()()TTxx

5、 AxAAx()2, ATxx AxAx其中 为对称阵（3-2-2）当A不时对称阵时，式（3-2-10）应该为（3-2-3）利用式（3-2-1）和（3-2-2）可以立即得到220TTxgA AxA b（3-2-4）TTA AxA b 这就是书中例2-4-1中所得到的法方程若使用配方法，则有：11( )2 () () () () minTTTTTTTTTTTTTTTTg xx A Axb Ax b bA AxA bA AA AxA bb b b A A AA bA AxA b可以看出，1min()TTTTgb bb A A AA b 本例中介绍的两个向量求导公式中，提到了对于向量x求导的梯度算符

6、 ,我们还可以引入对矩阵 求导的梯度算符 ：xijAaA1112112nAnnnnaaaaaa （3-2-5）需要说明的是，算符 只有作用在关于 的标量函数上才有意义。例如对于二次型由于 ，故Aija1112,1(,)nTnnijiji jf aaaa x xx Ax（3-2-6）ijijfx xa,1,.,()TTAiji jnx Axx xxx（3-2-7） 在课本中，给出了一些常用的向量-矩阵求导公式，在实际应用中可供大家查阅。 设有如图3-3-1所示的系统T。当输入n个数据 时，输出为y，且有下列线性关系：3-3 应用举例应用举例3-3-1 系统辨识系统辨识12,nx xx1 122n

7、nya xa xa x12,na aa（3-3-1）其中 为未知，需要通过对输入输出的观测值来确定这组参数。 现设进行了m次观测，观测值为 和 12( ),( ),( )nx kx kx k1x2xnxy图3-3-1 多输入单输出系统( ),1,2,y k km则问题成为求 使之满足12,na aa1 122( )( )( )( ) 1,2,nna x ka xka xky kkm（3-3-2）若记 及则方程（3-3-2）成为1 (1), (2),., ( ) ,.,TTnyyyy maaa11(1)(1)( )( )nnxxXx nx n Xay（3-3-3）当方程（3-3-3）无解时，问题

8、就转化为求矛盾方程组的最小二乘解。可以得到1()TTaX XX y（3-3-4）进而考察多输入多输出的情形。关系式为其中1122 1,2, iiininTTya xa xa xipyx A或（3-3-6）（3-3-5）1112,., ,.,TTTpnyy yyxx xx1111pnnpaaAaa 现设输入和输出的第k次观测值分别是 ( )( ),1,2,TTxkk km和y则系统的辨识问题就是求A使之满足(1)(1)()() Y = X ATTTTyxAymxm或 简 记 为（3-3-8）（3-3-7）其中Y为 矩阵，X为 矩阵。当上述方程无解时，问题就转化成为求A使下列非负定矩阵达到极小：m

9、pm n-minTJ AYXAYXA（ ）（） （）（3-3-9）问题（3-3-9）可以用配方法来求解：1111 - - minTTTTTTTTTTTTTTTTTTJ AY YA X YY XAA X XAY YY X X XX YAX XX YX X AX XX YAX XX Y（ ）（）（）（） （）（）（3-3-10）其中假定 可逆。TX X 这个问题不能用求导法来求解，因为目标函数J（A）不是标量而是矩阵。要用求导法来求解该问题，需要引入矩阵范数的概念。矩阵 的范数定义为m pYR211|pmTijijYtr Y Yy （）（3-3-11）事实上，可以把式（3-3-11）理解成向量 的

10、范数。1112121,Tppmmpyyyyyy， 这样，我们可以把多输入多输出线性系统的辨识问题叙述为求 矩阵A，使得np21( ) | () ()minTJAYXAtr YXAYXA（3-3-12）式（3-3-12）的形式与（3-3-9）类似，但应注意在此处 是标量函数。她可以完全类似于式（3-3-10）那样来配方而求解，也可体用求导法来求解。由于1( )JA1( )() 2 ()()TTTTJ Atr Y Ytr Y XAtr A X XA（3-3-13）利用课本中表3-2-2中的公式5和7，得到11( )220 ()TTATTJAX YX XAAX XX Y 即（3-3-14） 数据压缩

11、是指在传输或存储信号时对信号数据量进行压缩。实际中的信号往往都是维数很高的随机数据向量。各种数据间的相关性也很大，简单的随意压缩会导致数据严重失真。按照最优化原则设计的数据压缩技术可以解决通讯和数据传输系统的信道容量不足和计算机存储容量不足的问题，因而是一种从容量方面提高系统使用效率的重要技术。3-3-2 数据压缩数据压缩下面向大家介绍一种有效的数据压缩方法。其思想是对信号作正交变换，根据失真最小的原则在变幻域进行压缩。其框图如下所示。1TT压缩xyyx 设 为n维随机变量，n的值很大。经过正交T变换后，得到变幻域的n维向量1 ,Tnxxx1,TnyyyyTx（3-3-15）其中变幻矩阵T的列

12、向量为满足12,n 12,TnT （3-3-16）Tkmkm （3-3-17）现在对数据进行压缩，即保留y的m个分量其余的n-m个分量用预先选定的常数 代，替，得到一个新的向量，再由 通过逆变换 得到x的估计kb11,Tmmnyyybby1T1xTy我们的问题是：如何选取变幻矩阵T和常数 ，y的那些分量被压缩掉，才能使 最接近x，即均方误差最小：kbx221|() minniiixxExx（3-3-18）下面来注意解决这些问题。由式（3-3-15）和式（3-3-16）可知1121,nniiixTyyy 111mniiiiii mxTyyb由此式可知，要使 达到最小，应有2211|()|() n

13、niiiiii mi mE ybEyb（3-3-19）2 ()0kkkE ybb （3-3-20）即,1,kkbE ykmn 这就使常数 应选取的数值。由此式得kb TTkkkbExE x（3-3-21）其中 是随机向量x的协方差阵。由式（3-3-22）可考虑如何选择T的行向量使 达到极小。注意到 是归一化正交的。为求 在条件 下的极值，令代入式（3-3-19），得到11()() nnTTiiiiixii mi mEybybC（3-3-22）xCii1Tii 1(1)nTiii m 这表明 是协方差阵 的特征值，而 是相应的特征向量，很显然,T取成x的卡享南-洛厄维变换是最合理的，此时的最小均

14、方误差为则由式（3-3-22），有1(1)nTTiXiiiii mC 于是 1,220kxkkkkmnxkkkCC （3-3-23）kkCk根据此式，我们可先把非负定矩阵 的特征值按大小次序排列，然后根据实际问题对均方误差的要求选择m，使得 小于指定的误差（即选择满足此条件的最小的m），把n-m个较小的特征值 所对应的 换成 。min11nnTixiii mi mC（3-3-24）xCmin1nii miiy iibE y12n特别地，如果取Ex=0(通常的信号经过预处理后可满足此条件），则取 ，至此，前面提出的问题便全部解决。 0ib 在上述解法中，卡享南-洛厄维变换被选用并不是偶然的，因为

15、这种变换消除了原始信号x的诸分量间的相关性，从而使数据压缩能遵循均方误差最小的准则实施。上述数据压缩方法告诉我们应该压缩掉y中那些方差大的分量，这称为数据压缩的方差准则。 由于卡享南-洛厄维变换需要知道矩阵的特征值和特征向量，其计算量非常大，因此在实际应用中通常都使用固定程式的有限正交变换。尽管这些变换不是最佳的，但实践和理论表明它们也都能在较大程度上消除随机向量诸分量间的相关性，而且由于它们具有快速算法，因而是实用的。卡享南-洛厄维变换是理论上的最佳变换，它可作为理论研究的工具，也可用来衡量其他变换优劣的标准。 设有随机向量 ，它含有真实信号 及噪声 ，如下述模型所示：3-3-3 维纳滤波维

16、纳滤波1 ,Tnxxx1 ,Tnsss1 ,TnvvvxHsv（3-3-25）其中H为已知 矩阵，它也表示一种干扰。现在将x输入到一个如下图所示的系统n n1TTxyA1sTAy要使相应的输出 成为真实信号s的最佳估计，即均方差最小：求滤波器矩阵A，其中T是正交变换。s2|() ()minTssE ssss（3-3-26） 这个问题称为维纳(Wiener)滤波器设计问题。正交变换T的作用是把滤波问题转化到变换域处理。由系统框图知于是11sTAyTATx（3-3-27）利用表3-2-2中的公式2和3，得到11() () (2TTTTTTEsTAysTAyEs sy A Tsy A Ay2 2 T

17、TAE TsyE Ayy 因此最佳滤波器矩阵A应满足TTAE yyTE sy（3-3-28）当 可逆时TE yy1 TTATE syE yy（3-3-29）设信号与噪声的均值为零：Ex=0；Ev=0；真实信号与噪声不相关： ；并记真实信号s与噪声v的自相关矩阵为则0TTE svE vsTTE ssPE vvR () ()TTTTTTTTTTTTE syE sxTE s HsvTPH TE xxHPHRE yyT HPHR T于是式（3-3-29）便成为这就是维纳滤波器矩阵。1()TTTATPHHPHRT（3-3-30） 下面再用配方法来解决此问题。为此，用矩阵的迹来表示 是有效的。1111()

18、 () ()() ()() TTTTTTTTTE s ss strE s s s strE s T Ay s T Aytr E ssT AE yy AT E sy AT T AE ys又一次得到了式（3-3-29），而且得到了最小误差：1111111 ( ) min TTTTTTTTTTTTtr E sstr E syE yyE ystr TAE syE yyE yyTAE syE yyTAE syE yy1min ()TTTTTTtr E sstr E syE yyE ystr PP HPHR P（3-3-31）注意到 的表达式中不出现正交变换T，由此可知维纳滤波与正交变换的选取无关。于是由

19、式（3-3-30）联想到若取T使得矩阵A呈对角型，则将使矩阵乘向量的计算量由 降低到n。这种想法是有实际意义的。如在图象处理中，（3-3-25）是一个典型的失真模型。若选T为傅立叶变换，则A便成为对角型矩阵。这样的维纳滤波器称为对角线维纳滤波器。min2n 模式识别问题是设计一种分类器，使之能自动地将类别未知的对象进行归类。需要归类的对象称为模式。 本节来介绍模式分类器的一种设计方法。一个需要识别的对象可以用向量来表示。下面的分类器的设计是对压缩后的模式向量进行的。 3-3-4 模式识别模式识别分类器的设计是通过下列方式（称为训练或学习）来完成的。被识别的对象共有k类： 。首先选取一些类别已知

20、的模式向量（称为训练模式）1,kCC,1,2,1,2,ijiixCjN ik求属于同一类 的模式向量的平均模式：iC11,1,2,iNiijjixxikN那么，对于任一类别未知的模式向量x，判决准则是 称为判决函数。根据判决准则(3-3-33)实现分类的分类器称为最小距离分类器。由于实际问题中属同一类的模式往往分布的很分散和凌乱，直接采用 这种线性函数作为判决界容易造成误判。0min | |iiixxxxxC若，则判定（3-3-32）由于222| -2|Tiiixxxx xx若记2( )2|Tiiig xx xx则判决准则等价于若00max( )( )iiig xgxxC，则判定（3-3-33

21、）( )ig x( )ig x为了克服上述缺点，可采用如下所述的最小二乘映射技术。 取定k个向量 。希望存在一个线性变换A（矩阵）把属于同一类 的训练模式 ，都变为 ：1,kvviC1,iiiNxxiv,1,1,ijiiAxvjN ik于是问题转化为求矩阵A，它能够把 都变到 的附近，越近越好。记ijxiv2|ijijiAxv则矩阵A关于第i类训练模式的平均偏差为11iNiijjiN总偏差为2111111|1 (2)iiNkkiijiiijiNkTTTTTijijijiiiijiAxvNx A Axx A vv vN于是A的选取应使 ，故A应满足min0 A即（3-3-34）111(22)0i

22、NkTTijijiijijiAx xv xN因此1111111()()iiNNkkTTiijijijijijiiAv xx xNN（3-3-35）求出A以后，对于任一类别未知的模式x，可如下述进行判决归类：00min | |iiiiAxvAxvxC若，则判定（3-3-36）这种分类判决的实现称为最小二乘最小距离分类器。它的主要思想是利用最小二乘映射把属于同一类的训练模式尽可能地聚集在一起。从而克服了最小距离分类器的缺点。 在这里所有 都是可以任意选择的，可以通过选择特殊的 ，简化判决函数的获得过程。在课本中就给出了一种简化方法。其中取 为K维单位向量iviviv0,0,1,0,0Tiv 可以看

23、出，最小二乘最小距离分类器的判决函数仅依赖于最优矩阵A。一旦求出A，立即就可实现这种分类器。在实际应用中此方法的设计技巧是选取合适的判决向量 。iv 设X为希尔伯特空间， 是X中的一组线性无关元（不一定正交）。对某一 ，求 ，使得3-4 法方程法方程1,nyyxX1 ,nxMspan yy| min |m Mxxxm（3-4-1）换言之，求系数 ，使得1,naa211(,) |minnnkkkf aaxa y（3-4-2）这就是最优系数 应满足的方程，它是一个线性代数方程组。具体写出来是利用投影定理， 应满足ka1,1,kkmkxa yymn即1(,)( ,)1,nkmkmkyyax ymn（

24、3-4-3）1,naa1121111122222212( ,)(,)(,)( ,)( ,)(,)(,)( ,)( ,)(,)(,)( ,)nnnnnnnny yy yy yax yy yy yy yax yy yy yy yax y（3-4-4）矩阵Y就是元素组 的格拉姆矩阵。它是非负定的，且当 线性无关时Y可逆。方程（3-4-4）称为Wiener-Hopf方程。一般地，一个最优化解应满足的方程称为该最优化问题的法方程。因此Wiener-Hopf方程是最优化问题（3-4-2）的法方程 。或写成矩阵-向量形简洁形式：Ya（3-4-5）1 ,nyy1,nyy 现设 是最优系数，即满足（3-4-3）

25、，这时的最小误差为：1,naa2211111| (,) (,) ( ,)(,)nkkknnkkmmkmnkkknkkkxa yxa yxa yxa yxx xayx（3-4-6）把这个方程合并到方程组（3-4-4）中去，成为其中detG表示矩阵G的行列式。2111212112112221122( ,)(,)(,)0( ,)( ,)(,)(,)0( ,)( , )(, )(, )0( , )nnnnnnnnnny y ay y ay y ax yy y ay y ay y ax yy x ay x ay x ax x用克莱姆法则即可求得最小误差为211det(, )det(,)nnG yyxG

26、yy（3-4-7） 特别地，如果 为归一化正交系，则 为单位矩阵，式（3-4-4）成为 ，式（3-4-7）成为1 ,ny y1 ,nGyy( ,)kkax y2221|nkkxa 上述问题利用求导法也可以解。由式（3-4-2），有1121112(,) |2( ,)(,) |2nnkkmmkmnnnkkkmkmkkmTTfxa yxa yxx yayya axY 利用求导公式， 应满足220fY 此即式（3-4-5）若用配方法，则有21121min( ) |()() min0|TTTfxYYYYYfxY易知它与式（3-4-6）是等价的。下面，我们再来讨论一个随机序列的预测问题，来进一步说明法方程

27、的意义考虑二阶矩有限的希尔伯特空间 中的序列 。记子空间12 ,x x ,11,k Nk Nk NkMspan xxx（3-4-8）现在的问题是：用 中的元素,k NM()1NNkmk mmxa x来估计 ，使得均方误差最小。也就是求系数 使kx1,Naa2()2|() minNNkkkkxxExx这个问题称为随机序列的预测问题。（3-4-9）根据投影定理， 应是 在子空间中的投影，即 满足()Nkxkx（3-4-10）,k NM1,Naa1,1,Nkmk mk lmxa xxlN根据空间 中的正交性定义，上式即为1,1,Nk mk lmkk lmE xxaE x xlN这就是最佳预测的法方程

28、。其中 是该平稳序列的自相关。它满足 。（3-4-11）又称为关于平稳序列 预测问题的Yule-Walker方程。其分量形式为如果随机序列 是平稳的，则式（3-4-10）成为12 ,x x 1,1,Nm lmlmrarlN（3-4-11）rrm rmrE xx12 ,x x 0011111022212NNNNNNrrrrarrrrarrrra（3-4-12）相应的预测误差为()2()()()1(, ) | (,) (,) NkkNNkkkkNkkkNkkkk mmmN kxxxxxxxxxE x xE x xa01 Nmmmrr a把这个方程合并到方程（3-4-12）中去，可以写成001101

29、1121(, )00NNNNNrrrN krrrarrra 记其系数矩阵为 。利用克莱姆法则可得NR1det(, )detNNRN kR（3-4-13）（3-4-14） 对于给定的序列x和y，求一个离散时不变系统h，使得当输入为x时相应的输出恰为y。这个问题的数学描述是：给定x和y，求h使得3-4 最小二乘滤波最小二乘滤波x hy（3-5-1）一般来说，这个问题不一定是有解的，因此将问题改成：已知x，y，求h使得| minx hy （3-5-2）我们假定x，y，h都是因果的，能量有限的于是x，y，h ，式（3-5-2）中的范数是 意义下的范数。这个问题称为最小二乘滤波问题。我们先用求导法来解。

30、由2l2l00nn kkn kkkkx hxhxh知目标函数为212200( ,) | ()nn kknkf h hyx hyxh 于是最优解 应满足12 ,hh h002()0nn kkn mnkmfyxhxh即000n kn mknn mknnxxhy x若记00,1,2,1,2,xnnnxynnnxxyx它们分别称为x的自相关序列和x与y的互相关序列。则式（3-5-3）成为（3-5-3）这就是最小二乘滤波的法方程。 下面再利用投影法来求解最小二乘滤波问题。在空间 中讨论。对于引入移位算子T：00,1,2,xxym kkmkhm（3-5-4）2l012,xx xl012010,0,0, T

31、xx xT xx x 则000110021 120200120, kkkx hh x h xh x h xh xh xh xhTxh T xh T x这样，就把卷积表示成 中一组元素 的线性组合，这是应用投影法的关键。记子空间 。问题便转化为求 ，使 ，也就是求y在M中的投影。根据投影定理， 应满足 2 ,x Tx T x 2l2 ,Mspan x Tx T xkh0| minkkkyh T x0,0,1,2,kmkkyh T xT xm即00(,)0,0,1,2,(,)( ,),0,1,2,kmkkkmmkkyh T x T xmh T x T xy T xm（3-5-5）kh因此式（3-5

32、-5）就是法方程（3-5-4）。接下来我们介绍两个实际应用的例子。其中0120101 1220( ,)(,0,0,)mmmmmxym nnmny T xyy yx xy xyxyxyx 个同理0(,)kmxn kn mm knT x T xxx 仍然考虑因果序列。设 是原始声音信号序列。由于回声的影响，在n时刻接收到的信号 除了 以外还叠加有 时刻之前的回音：例例 3-5-1回声干扰与最小二乘逆滤波 nssnuns0nnnn uusrs（3-5-6）其中r是衰减因子， 。现在要求设计一个系统 能够消除这种回声干扰，即| | 1r nhhu hs（3-5-7）我们把式（3-5-6）写成卷积形式：

33、其中 为usx （3-5-8）nxx01,0,0,nnxrnn其他将式（3-5-8）代入式（3-5-7）得这等价于sx hs x h（3-5-9）其中 是单位脉冲序列。于是问题转化成给定x，要求设计一个系统h，使之满足式（3-5-9）。这是一个特殊的最小二乘滤波问题，称为最小二乘逆滤波。利用法方程（3-5-4），取 ，即得y0,0,1,2xxm kkmkrhrm（3-5-10）而根据序列x的因果性，有00,00,0 xmn mnmnxmrxxm如果把方程（3-5-10）中的m截断到N为止，则成为000111011000 xxxNxxxNxxxNNNhxrrrhrrrhrrr（3-5-11）这就是最小二乘逆滤波的“N截断”法方程，它适合于实际计算