0 前言
1、奇异值分解,Singular Value Decomposition,SVD
2、总体最小二乘法,Total Least Squares,TLS
3、为叙述方便,使用了与原书——《现代信号处理-张贤达》中相同的公式序号,读起来可能有点不方便,请见谅。
1 奇异值分解SVD
奇异值分解(SVD)主要用于求解线性方程组。
【定理4.1】:令A是一个m×n维复数矩阵,则分别存在一个m×m维和一个n×n维酉矩阵U和V,使得:
A=UΣV*
上式中的上标*表示矩阵的共轭转置,Σ是一个m×n维对角阵,且Σ矩阵的主对角线上的元素是非负的,并按此顺序排列:σ11≥σ22≥…≥σhh≥0,式中h=min(m,n)。
矩阵U和V称为酉阵,若它们分别满足U-1=U*和V-1=V*。
一个实的酉阵U又叫正交阵,即有U-1=UT。对角元素σkk称为矩阵A的奇异值,U和V分别叫做矩阵A的左奇异阵和右奇异阵,而U=[u1,u2,…,um]和V=[v1,v2,…,vn]的列向量ui和vj分别称作矩阵A的左奇异向量和右奇异向量。
非零的奇异值对应于非负Hermitian阵AA*和A*A的特征值的正平方根;而且U(或V)的列对应于非负Hermitian阵AA*(或A*A)的适当排列的正交特征向量。
有两种方法用于确定m×n维复数矩阵A的有效秩,分别为归一化比值法与归一化奇异值法。
1.1 归一化比值法
考虑这样一个问题:求秩为k的m×n维矩阵,使得它在Frobenious范数意义上可最佳逼近A,这儿假定k≤rank(A)。
m×n维矩阵差A-B的Frobenious范数定义为:
现在的问题转化为:【寻找一个m×n维且秩为k的矩阵B,从而使上述范数为最小】。这一逼近问题的解可以用下面的定理来描述。
【定理4.2】:在Frobenious范数意义下能最佳逼近m×n维矩阵A的唯一m×n维,且秩k≤rank(A)的矩阵由
A(k)=UΣkV*
给定,其中U和V由本节即【1 奇异值分解】最开始的定理中给出,而Σk是通过在Σ内令除k个最大的奇异值以外的所有其它奇异值都等于零后得到的对角阵。这一最佳逼近的质量由:
描述。定理4.2表明,A(k)逼近A的程度取决于(h-k)个最小的奇异值的平方和。当k接近h时,这一求和将变得越来越小,并在k=h时最终趋于零。考虑下面的归一化比值:
在实际应用中,通常选择一个接近于1的数作为门限值(例如0.995),并把v(k)大于该门限值的最小k值定为矩阵A的有效秩。
1.2 归一化奇异值法
除去使用归一化比值v(k)外,也可以利用归一化的奇异值来确定矩阵A的有效秩。把
叫做归一化的奇异值。与使用归一化比值v(k)的情况相反,利用
确定有效秩时,选择一个接近零的数作为门限值(例如0.05),并把
小于此门限的最大k值取为矩阵A的有效秩。
2 总体最小二乘法TLS
在TLS问题中,考虑矩阵方程(e、E统一被称为“扰动”):
(A+E)x=b+e
的求解。上面的式子可写为下式(记为4-1-14):
于是上图中的齐次方程的总体最小二乘解可以表示为:求解向量z使得:
【||D||F=min】与【(b+e)∈range(A+E)】
其中:
且:
TLS方法使来自A和b的噪声扰动影响最小。TLS问题可以归结为:求一具有最小范数的扰动矩阵D∈Rm×(n+1)使得B+D是非满秩的(如果满秩,则只有平凡解z=0)。
SVD可用于实现上述目的。计算矩阵B的SVD分解B=UΣV*,且奇异值仍按递减次序排列:σ1≥σ2≥…≥σmin(m,n+1)≥0。其中增广矩阵B(前面已列出,下面只是重复列出而已)为
有如下情况:
一、当m<n+1时,方程(4-1-14)是欠定的,存在无穷多个解,而TLS方法可给出最小范数解;
二、当m>n+1时,方程(4-1-14)是超定的。这里仅考虑超定方程组的总体最小二乘问题的求解。在超定方程组的TLS求解中,又有两种可能的情况:
1、第一种情况是σn明显地比σn+1大。此时,存在唯一的TLS解,由:
给出,式中v(i,n+1)是右奇异阵V的第n+1列向量中的第i个元素。
2、第二种情况是B的后面若干个奇异值是重复的(实际应用中常表现为多个奇异值很接近)。此时,TLS会给出多个解。不妨令:
σ1≥σ2≥…≥σk>σk+1≈…≈σn+1
且vi是子空间S=span(vk+1,…,Vn+1)中的任一列向量,那么,TLS解x将由:
x=yi/ai
给出,其中vi=[ai/yi],其中yi是向量,ai是标量,i=k+1,…,n+1,且ai≠0。换句话说x^共有n+1-k个解,它可以由vk+1,…,vn+1中的任何一个列向量求出。在这种情况下,某种意义上可能的唯一解有两种:最小范数解和最佳最小二乘近似解。
2.1 最小范数解
计算TLS问题的最小范数解的算法由以下四个步骤组成:
观察步骤(3)可以发现,y可以通过使Q的第一列为V1第一行的复数共轭来获得(还有其它各种办法,但这是最简单的一种)。由于确定解xTLS只需要使用:
所以只要计算Q的第一列即可。于是可将xTLS表示成(可以自己推一推,确实是这样):
a=0的情况对应于
因此,在a=0或a值太小的情况下,应该减小原定的指数p,也就是扩大矩阵V1的列数,直到得到一个非零或不接近零的行向量
指数p的确定实际上就是前面【1.1 归一化比值法】、【1.2 归一化奇异值法】两节中所述的有效秩确定,因此步骤(2)中可改用上述两节的归一化比值或归一化奇异值去确定p。
2.2 最佳最小二乘近似解
首先说明下,下述矩阵B是【2 总体最小二乘法TLS】一节最开始的式4-1-14的增广矩阵,有:
令(p+1)×1维向量a为:
于是【2 总体最小二乘法TLS】一节的原方程组(4-1-14)的求解就变成了下列n+1-p个方程组的求解:
不难证明:
记上式为4-1-25。其中vki定义为:
上图中的v(i,j)是矩阵V的第i行第j列上的因素。uki也是这样定义的。根据最小二乘的原理,方程组4-1-25的最小二乘解等价于使测度函数:
最小化。定义(p+1)×(p+1)维矩阵(记为4-1-28):
则测度函数可以写为:
f(a)的最小化就是求解∂f(a)/∂a=0,于是可给出下列线性方程组(记为4-1-30):
其中e1=[1,0,…,0],归一化常数α(阿尔法)的选择应使a的第一个分量等于1。由定义式(4-1-28)及(4-1-25)可知,显然有:
上式不含任何左奇异向量uj。对于方程组(4-1-30),令S-(p)是Sp的逆矩阵,则解a仅取决于逆矩阵S-(p)的第一列。最终解由:
给出。下面归纳求最佳最小二乘近似解的算法的步骤:
其中(3)的提到的(4-1-25b)应该为:
上图中的v(i,j)是矩阵V的第i行第j列上的因素。
3 参考文献
END