奇异值分解SVD与总体最小二乘法TLS

---

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维对角阵，且Σ矩阵的主对角线上的元素是非负的，并按此顺序排列：σ₁₁≥σ₂₂≥…≥σ_hh≥0，式中h=min(m,n)。

矩阵U和V称为酉阵，若它们分别满足U^-1=U*和V^-1=V*。

一个实的酉阵U又叫正交阵，即有U^-1=U^T。对角元素σ_kk称为矩阵A的奇异值，U和V分别叫做矩阵A的左奇异阵和右奇异阵，而U=[u₁,u₂,…,u_m]和V=[v₁,v₂,…,v_n]的列向量u_i和v_j分别称作矩阵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∈R^m×(n+1)使得B+D是非满秩的(如果满秩，则只有平凡解z=0)。

SVD可用于实现上述目的。计算矩阵B的SVD分解B=UΣV*，且奇异值仍按递减次序排列：σ₁≥σ₂≥…≥σ_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会给出多个解。不妨令：

σ₁≥σ₂≥…≥σ_k>σ_k+1≈…≈σ_n+1

且v_i是子空间S=span(v_k+1,…,V_n+1)中的任一列向量，那么，TLS解x将由：

x=y_i/a_i

给出，其中v_i=[a_i/y_i]，其中y_i是向量，a_i是标量，i=k+1,…,n+1，且a_i≠0。换句话说x^共有n+1-k个解，它可以由v_k+1,…,v_n+1中的任何一个列向量求出。在这种情况下，某种意义上可能的唯一解有两种：最小范数解和最佳最小二乘近似解。

2.1 最小范数解

计算TLS问题的最小范数解的算法由以下四个步骤组成：

观察步骤(3)可以发现，y可以通过使Q的第一列为V₁第一行的复数共轭来获得(还有其它各种办法，但这是最简单的一种)。由于确定解x_TLS只需要使用：

所以只要计算Q的第一列即可。于是可将x_TLS表示成(可以自己推一推，确实是这样)：

a=0的情况对应于

因此，在a=0或a值太小的情况下，应该减小原定的指数p，也就是扩大矩阵V₁的列数，直到得到一个非零或不接近零的行向量

指数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。其中v_kⁱ定义为：

上图中的v(i,j)是矩阵V的第i行第j列上的因素。u_kⁱ也是这样定义的。根据最小二乘的原理，方程组4-1-25的最小二乘解等价于使测度函数：

最小化。定义(p+1)×(p+1)维矩阵(记为4-1-28)：

则测度函数可以写为：

f(a)的最小化就是求解∂f(a)/∂a=0，于是可给出下列线性方程组(记为4-1-30)：

其中e₁=[1,0,…,0]，归一化常数α(阿尔法)的选择应使a的第一个分量等于1。由定义式(4-1-28)及(4-1-25)可知，显然有：

上式不含任何左奇异向量u_j。对于方程组(4-1-30)，令S^-(p)是S^p的逆矩阵，则解a仅取决于逆矩阵S^-(p)的第一列。最终解由：

给出。下面归纳求最佳最小二乘近似解的算法的步骤：

其中(3)的提到的(4-1-25b)应该为：

上图中的v(i,j)是矩阵V的第i行第j列上的因素。

---

3 参考文献

现代信号处理-张贤达

---

END