一种改进的调和Arnoldi算法及其应用

肖小花

(西安理工大学 理学院， 陕西 西安 710054)

0 引言

20世纪50年代以来，求解大规模矩阵特征值问题主要采用的方法为投影类方法.目前，投影类方法中出现了很多比较有代表性的算法，其中Arnoldi-type算法[1]被公认为是投影类方法中最成功的方法之一.传统Arnoldi方法[2]虽然对求解大规模矩阵的端部特征对有着很好的效果，但是对矩阵的内部特征值问题常常失效.后来paige等人提出的调和Arnoldi方法已经成为求解大规模矩阵特征对，尤其是内部部分特征对的最重要的方法之一.

本文根据投影子空间中所含的特征信息越丰富，该子空间上Ritz对的近似程度越好的指导思想，结合调和Ritz值收敛而调和Ritz向量却难以收敛的情况，充分利用Arnoldi过程产生的第m+1个基向量提供的有用信息，在m+1维的Krylov子空间中来求解近似向量，使Ritz对的残量范数达到极小.并对这种方法进行了理论分析，给出了数值实验.理论分析显示这种方法的可行性；数值实验显示这种方法的有效性.最后将本文的算法应用于K-L变换的变换矩阵求解中，K-L变换的核心过程是计算特征值和特征向量.由于待处理矩阵维数高的情况，一般的方法很难求出其特征值及特征向量，甚至无法求出.而K-L变换[5]的一些优化处理过程复杂，很难满足实时性的要求，使得用K-L变换进行图像压缩难以广泛应用.而本文的算法正好就是一种较快求解大规模矩阵特征值及特征向量的方法，实验表明将其用于K-L变换来进行图像压缩是可行的.

文中Km(A,v1)=span{q0,Aq0,…,Am-1q0}表示m维的Krylov子空间，上标*表示矩阵或向量的共轭转置，I表示N阶单位矩阵，em表示m阶单位矩阵IM的第M列，σmin(G)表示G的最小奇异值.文中的范数‖*‖如无特殊说明均为2范数.

1 调和Arnoldi方法

大规模矩阵特征值问题：

Aφi=λiφi(1)

的数值求解，其中A为n阶实或复矩阵，(λi,φi)i=1,2,…,n为A的特征对(‖φi‖=1).

给定m维Krylov子空间Km(A,v1)由Arnoldi过程产生它的一组标准正交基，这一过程的矩阵表达式为：

(2)

(3)

(4)

2 改进的调和ARNOLDI算法

考虑大规模矩阵特征值问题：

Aφi=λiφi(5)

(6)

则有矩阵表达式：

(7)

给(7)式两边同时左乘(A-τI)得到：

(8)

(0代表零行向量)

所以有：

从上面的讨论，可得以下结论:

根据以上过程，下面给出整个算法如下:

(1)给定子空间维数m，需要计算特征向量的个数k以及要求达到的精度tol，选择一个单位初始向量q0以及位移点τ；

(5)重新启动：构造新的初始向量q0，转向第(2)步.

注:算法第(5)步采用显式重新启动策略，即重新启动后的初始向量q0取作：

其中α为规范化因子.

3 数值算例与实验结果

3.1 数值算例

图1

图2

数值算例1．以图1中的图像矩阵作为待计算特征值和特征向量的矩阵，矩阵的规模为1000×1 000.按本文方法计算得到矩阵按实部最大的3个特征值依次近似为：λ1≈12.3285，λ2≈9.5228,λ3≈8.2896.表1给出了计算结果，其中m表示子空间的维数，iter表示重新启动的次，time表示运算时间(秒).mv表示矩阵与向量乘积的个数.tol=10-7,位移τ=1，初始向量q0是按均匀分布生成的随机向量.

表1 数字图像转换生成的1 200×1 200矩阵

数值算例2．此例中是以图2中的图像矩阵作为待计算特征值和特征向量的矩阵,矩阵的规模为2 400×2 400.按本文方法算得矩阵按实部最大的三个特征值依次近似为：λ1≈13.358 1，λ2≈10.807 9，λ3≈7.394 2.表2给出了计算的结果，其中m、iter、 time、mv的表示同数值算例1.位移τ=1，精度tol=10-6，初始向量q0是按均匀分布生成的随机向量.

表2 对2 400阶方阵用不同方法计算的结果比较

由以上两例的计算结果可以得出：当子空间维数越小时，本文算法达到收敛迭代的步数越少，优越性比较明显.随着子空间维数的增加，子空间中含有的特征信息比较丰富，则两种方法的收敛速度都比较快且比较接近.

3.2 实验

在用K-L变换对图像进行压缩的实验中，需考虑计算机处理器的性能，如果计算机处理器的性能比较好，则用本文的算法就可以直接对整幅图像进行K-L变换来压缩图像；如果计算机处理器的性能不是太好，则可以把图像矩阵分成大小相同的几个矩阵，分别对每个矩阵进行K-L变换，再把得到的每个压缩图像矩阵相应的合成为压缩后的整幅图像矩阵.对于将图像分成若干小块而对每个小块分别进行K-L变换的方法[9]，一般选择大小为8×8的块.

实验1对大小为200×200×8 bit的灰度图像shanshui.jpg进行压缩.比较了以下两种方法：第一种是本文方法即将本文的算法直接用于图像矩阵进行K-L变换；第二种是将图像分成若干小块，而对每个小块分别进行变换[10]的方法，每个小块选为8×8.在实验中，本文的算法选取Krylov子空间 的维数m=30 ,V表示算法的执行速度(以分为单位).对特征值的保留个数分别取为1、2、4、8、16，将特征向量量化为8位二进制数，用两种方法分别进行测试的情况如表3所示.

表3 200×200×8 bit图像的压缩比和峰值信嗓比及算法执行速度

用本文算法对shanshui.jpg进行压缩的原图和保留特征值个数分别为1，2，4，8，16的压缩后的重建图像如图3所示.

将图像分块的方法对shanshui.jpg进行压缩的原图和保留特征值个数分别为1，2，4，8，16的压缩后的重建图像如图4所示.

图3 shanshui.jpg原图及保留特征值个数为1、2、4、8、16的重建图像

图4 shanshui.jpg原图及保留特征值个数为1、2、4、8、16的重建图像

图5 fengye.jpg原图及保留特征值个数为1、2、4、8、16的重建图像

图6 fengye.jpg原图及保留特征值个数为1、2、4、8、16的重建图像

实验2对大小为1 000×1 000×8 bit的灰度图像fengye.jpg进行压缩,比较了以下两种方法：第一种是本文方法即将本文的算法直接用于图像矩阵进行K-L变换,将图像矩阵分成25个200×200的方阵；第二种是将图像分成若干小块而对每个小块分别进行K-L变换,每个数据块选为8×8.在实验中，本文的算法选取Krylov子空间的维数m=30 ,V表示算法的执行速度(以分为单位).对特征值保留个数分别取为1，2，4，8，16时，将特征向量量化为8位二进制数，用两种方法分别进行测试的情况如表4所示.

表3 200×200×8 bit图像的压缩比和峰值信嗓及算法执行速度

用本文算法对fengye.jpg进行压缩的原图和保留特征值个数分别为1，2，4，8，16的压缩后的重建图像如图5所示.

将图像分块的方法对fengye.jpg进行压缩的原图和保留特征值个数分别为1，2，4，8，16的压缩后的重建图像如图6所示.

从以上两个实验的结果可以看出：从压缩比、峰值信嗓比以及算法的执行时间上，用本文的方法来使用K-L变换进行图像压缩要优于将图像数据矩阵分成若干小块而在每个小块上使用K-L变换的方法.

4 结束语

本文对调和Arnoldi方法进行了改进.针对往往调和Ritz值收敛而相应的调和Ritz向量近似程度很差的情况，保持调和Ritz值不变，充分利用基向量Vm+1提供的信息求解调和Ritz向量即改进的调和Ritz向量.理论分析和数值实验结果均表明了该方法的可行性及有效性，并将其用于K-L变换进行图像压缩.将这种方法引入K-L变换来进行图像压缩，解决了K-L变换中变换矩阵过大而求解困难的问题.对大规模矩阵特征值问题数值求解方法的讨论和研究，将有助于K-L变换在图像压缩中的广泛应用.

[1] Jia Z．Arnoldi type algorithms for large unsymmetric multipleeig envalue problems[J]．J Comput Math，1999，17：257-274．

[2] Saad Y．Variations on Arnoldi′s method for computing eigenelements of large unsymmetric matrices[J]．Linear Algebra Appl，1980，34：269-295．

[3] Sorensen D C．Implicit application of polynomial filters in a k-step Arnoldi method[J]．Siamj Matrix Anal Appl,1992,13：357-385．

[4] Jia Z．The convergence of harmonic Ritz values，harmonic Ritz vectors and refined harmonic Ritz vectors[J]．Math Comput,2005，74：1 441-1 456．

[5] 赵荣椿，赵忠明，崔苏生．数字图像处理导论[M]．西安：西北工业大学出版社，1999.

[6] Morgan R B．Zeng M．Harmonic projection methods for large non-symmetric eigenvalue problems[J]．Numer Linear Algebra Appl, 1998，5：33-55．

[7] Jia Z．The refined harmonic Arnoldi method and an implicitly restarted refined algorithm for computing interior eigenpairs of large matrices[J]．Appl Numer Math，2002，42：489-512．

[8] Morgan R B. Implicitiy restarted GMRES and Arnoldi methods for nonsymmetrics-tems of equations[J].Siamj Matrix Anal Appl,2000,21:1 112-1 135.

[9] 王文峰.K-L变换的研究及其在图像压缩编码中的应用[D].沈阳：沈阳理工大学，2008.