SVD(奇异值分解)算法及其评估

1SVD（奇异值分解）算法及其评估本文第一部分对SVD进行了简单的介绍，给出了定义和奇异值分解定理；第二部分简要地列举了SVD的应用；第三部分则构造和分析了各种求解SVD的算法，特别对传统QR迭代算法和零位移QR迭代算法进行了详细完整的分析；第四部分给出了复矩阵时的处理办法；第五部分是对各种算法的一个简要的总结。一、SVD简介定义1.1设mnAR×∈，TAA的特征值的非负平方根称作A的奇异值；A的奇异值的全体记作()Aσ[1]。当A为复矩阵mnC×时，只需将TAA改为HAA，定义1.1仍然成立。定理1.1（奇异值分解定理）设mnAR×∈，则必存在正交矩阵[]1,...,mmmUuuR×=∈和[]1,...,nnnVvvR×=∈使得000rTrnrrUAVmr−Σ=−，…(1.1)其中11(,...,),...0rrrdiagσσσσΣ=≥≥[2]。当A为复矩阵mnC×时，只需将定理中,UV改为酉矩阵，其它不变，定理1.2仍然成立[1]。称分解式(1.1)为矩阵A的奇异值分解，通常简称为SVD。iσ是A的奇异值，向量iu和iv分别是第i个左奇异向量和第i个右奇异向量。从A的奇异值分解，我们可以得到A的一些非常有用的信息，下面的推论就列举其中几条最基本的结论[1]：推论1.2设mnrAC×∈，则(1)A的非零奇异值的个数就等于()rrankA=；(2)1,...,rnvv+是()A的一组标准正交基；(3)1,...,ruv是()A的一组标准正交基；(4)1rHiiiiAuvσ==∑称为A的满秩奇异值分解；其中()A，()A分别指得是A的零空间和值域。为了方便，我们采用如下表示奇异值的记号：()AiAσ＝的第i大奇异值；max()AAσ＝的最大奇异值；min()AAσ＝的最小奇异值。现在来考察矩阵奇异值的几何意义[1,2]，不妨设mn=，此时有{}2:,,1nnnEyCyAxxCx=∈=∈=2是一个超椭球面，它的n个半轴长正好是A的n个奇异值12...0nσσσ≥≥≥≥，这些轴所在的直线正好是A的左奇异向量所在的直线，它们分别是对应的右奇异向量所在直线的象。一般地我们假设mn≥，（对于mn的情况，我们可以先对A转置，然后进行奇异值分解，最后对所求得的SVD分解式进行转置就可以得到原式SVD分解式），此时我们对(1.1)进行化简将U表示为：12(,),UUU=…(1.2)则可以得到更加细腻的SVD分解式[2,3]：1TAUV=Σ…(1.3)其中1U具有n列m维正交向量，V和(1.1)式中的定义相同；12(,,...,)ndiagσσσΣ=，并且12...0nσσσ≥≥≥≥为矩阵A的奇异值。二、SVD应用在现代科学计算中SVD具有广泛的应用，在已经比较成熟的软件包LINPACK中列举的应用有以下几点[3]:(1)确定矩阵的秩(rank)假设矩阵A的秩为r，那么A的奇异值满足如下式子11......0rrnσσσσ+≥≥===；反之，如果0rσ≠，且1...0rnσσ+===，那么矩阵A的秩为r，这样奇异值分解就可以被用来确定矩阵的秩了。事实的计算中我们几乎不可能计算得到奇异值正好等于0，所以我们还要确定什么时候计算得到的奇异值足够接近于0，以致可以忽略而近似为0。关于这个问题，不同的算法有不同的判断标准，将在给出各种算法的时候详细说明。(2)确定投影算子假设矩阵A的秩为r,那么我们可以将(1.3)式中的1U划分为以下的形式()(1)(1)112,UUU=其中(1)1U为mr×的矩阵，并且(1)1U的列向量构成了矩阵A的列空间的正交向量基。容易得到(1)(1)11TAPUU=为投影到矩阵A的列空间上的正交投影算子；而(1)(1)2222(,)(,)TAPUUUU⊥=则是到矩阵A列空间的正交补空间上的投影。同理如果将V划分为以下的形式()12,VVV=其中1V为nr×的矩阵，则1V的列向量构成矩阵A的行空间的正交向量基。且11TARVV=为投影到矩阵A的行空间上的正交投影算子；而22TARVV⊥=则是到矩阵A行空间的正交补空间上的投影。(3)最小二乘法问题(LS问题)3促进人们研究SVD并且应用SVD比较早的应该是最小二乘法问题了，在前一份关于QR分解的报告已经对LS问题进行了一些介绍，这里继续讨论SVD在LS问题中的应用。LS问题即相当于，设mnAR×∈(),mmnbR∈，求nxR∈使得22min{:}.nAxbAvbvR−=−∈……(2.1)假设已知矩阵A有式子(1.1)得到的SVD分解式为TUVΣ，U和V分别为,mn阶正交方阵，而∑为和A具有相同维数的对角矩阵，那么我们可以得到[4]：()()......(2.2)()TTTAxbUVxbUVxUUbUyc−=Σ−=Σ−=Σ−其中TyVx=，TcUb=。因为U是一正交矩阵，所以222()AxbUycyc−=Σ−=Σ−，从而把原最小二乘法问题化为求使2ycΣ−最小的y这一最小二乘法问题，因为Σ为对角矩阵，所以使得新的这一最小二乘法问题简单的多，接着将对此仔细分析[4]。假设矩阵A的秩为r，则有：11000rryyyσσΣ=，11112rrrrrmycycyccccσσ++−−Σ−=−−−可知/,(1,2,...,)iiiycirσ==使得ycΣ−达到它的最小长度1/221miirc=+∑，并且可见当rm=时，上面的这一长度为0，也就是当矩阵A的列张成m空间时最小二乘法问题可以无误差地求解。而当rn时，1,...,knyy+可以任意取，而不影响ycΣ−的长度。我们将对∑转置并且对非零的对角元素求逆所得到的矩阵定义为+∑，那么yc+=Σ的前r个元素将等于/,(1,2,...,)iicirσ=，并且其余的元素为0。并且由TyVx=，TcUb=，容易得到：......(2.3)TxVUb+=∑由此得到的是LS问题的最小范数解。而文献[3]中还给出了一般通解的形式如下：2......(2.4)TxVUbVw+=∑+其中2V如前定义，而w是任意的nr−维向量。4(4)广义逆问题(pseudo-inverse)记TAVU++=∑，从(2.3)式我们可以看出，最小二乘法的解为xAb+=，和一般的线性方程组Axb=的解为1xAb−=相类似，所以我们当我们已知矩阵A的奇异值分解TAUV=Σ后可以定义A的广义逆为TAVU++=∑。(5)条件数如果已知矩阵A的秩为r，那么在式子(2.3)中，解x随着矩阵A的扰动而改变的剧烈程度有多大呢？这可以用矩阵的条件数来衡量，条件数的定义如下：1()/......(2.5)rrAκσσ=以上只是针对SVD的应用，而简单地介绍了LS问题，广义逆；而在文献[5,6]中则对这两个问题有详细的说明，在以后的报告中也将进行更进一步的分析，并且综合其它方法来全面分析和处理这些问题。除了这些传统的应用以外，在图像压缩和大型数据库的数据恢复中，SVD也具有广泛的应用[7]。三、各种SVD算法及其特点首先来对SVD算法的发展来做简单的回顾[11,12]：关于SVD算法的研究最早可以追溯到1873年Beltrami所做的工作，这中间在理论方面进行了大量的工作，这个历史过程可以参考Stewart的文献[8]。但是直到1965年Golub和Kahan才在SVD的数值计算领域取得突破性进展[9]，并且于1969给出了比较稳定的算法[10]（以下简称传统QR迭代算法），这也是后来在LINPACK中所采用的方法[3]。它的中心思想是用正交变换将原矩阵化为双对角线矩阵，然后再对双对角线矩阵迭代进行QR分解。六十年代一份没有出版的技术报告中，Kahan证明了双对角线矩阵的奇异值可以精确地计算，具有和原矩阵元素的相对的精确度；进一步，1990年Demmel和Kahan给出了一种零位移的QR算法(zero-shiftQRalgorithm)，这种算法计算双对角矩阵的奇异值具有很高的相对精度[13]，并且由此得到的奇异向量也具有很高的精度[14]。Fernando和Parlett在1994年将qd算法应用到奇异值的计算上，从而得到了一种全新的比zero-shiftQRalgorithm更加精确和快速的计算奇异值的算法[15,16]。而Demmel和Veselic在文献[17]中则说明了用Jacobi方法与其它方法相比计算所得到的奇异值和奇异向量具有更高的精度，可惜Jacobi算法比DK算法速度要慢的多；在文献[18]中对Jacobi方法进行了改进使得其在速度上几乎和DK算法相当。和Strassen算法类似，SVD问题也有快速的分而制之算法，在文献[19,20]对其有详细的介绍，由此得到的算法在总计算量上比现有的LAPACK软件包中所采用的方法要少的多。在文献[21,22]中给出的bisection算法也可以对双对角线矩阵计算得到具有相对精度的全部奇异值。以下就开始对各种算法原理进行详细说明，并分析它们的计算量，计算的精确度，以及所占得内存。53.1：传统QR迭代算法[1,2,3]设()mnARmn×∈≥，可知奇异值分解可从实对称矩阵TCAA=的Schur分解导出[1]，因此我们自然想到先利用对称QR方法来实现C的Schur分解，然后借助C的Schur分解来实现A的奇异值分解，然而这样做有两个缺点：一是计算TCAA=要很大的计算量；二是计算TCAA=会引入较大的误差。因此Golub和Kahan在1965年提出了另一种十分稳定的方法，其基本思想就是隐含地应用对称QR算法于TAA上，而并不需要将TCAA=计算出来。方法第一步是：将A二对角化，即求正交矩阵1U和1V，使得110nTmnBUAV−=……(3.1.1)其中12230000000000000000nnBδγδγγδ=……(3.1.2)分解式(3.1.1)可以用Householder变换来实现，将A分块为[]1111nAvA−=先计算m阶Householder变换1P使得111111(,)mPveReRδδ=∈∈并且形成：111111TuPAmA=−再计算n－1阶Householder变换1H使得1112121(,)nHueReRγγ−=∈∈并形成：[]112212nAHvA−=然后对2,3,...,2kn=−依次进行：(a)计算m-k+1阶Householder变换kP使得111(,)mkkkkkPveReRδδ−+=∈∈并且形成：11TkkkkuPAmA=−(b)计算n－k阶Householder变换kH使得61111(,)nkkkkkHueReRγγ−++=∈∈并形成：[]1111kkkknkAHvA++−−=进行到2kn=−之后，再计算2mn−+阶Householder变换矩阵1nP−使得2111111(,)mnnnnnPveReRδδ−+−−−−=∈∈并形成：1111nnnnPAvmnγ−−=−+然后计算1mn−+阶Householder变换矩阵nP使得111(,)mnnnnnPveReRδδ−+=∈∈。现令()1,,2,...kkkPdiagIPkn−==(),,1,2,...2kkkHdiagIHkn==−112...nUPPP=，1122...nVHHH−=12230000000000000000nnBδγδγγδ=则有：110nTmnBUAV−=即实现了分解(3.1.1)。将A二对角化以后，下一步任务就是对三对角矩阵TTBB=进行带Wilkinson位移的对称QR迭代，这一步也可以不通过明确地将T计算出来而进行。进行QR迭代的第一步是取矩阵TTBB=的右下角22×主子阵：22111221nnnnnnnnδγδγδγδγ−−−−++靠近22nnδγ+最近的特征值作为位移µ，这一步不需将TTBB=计算出来。第