张亚蕾,杨少静
1.仰恩大学数学系,福建 泉州 362014;2.河北科技学院公共课部,河北 保定 071000
现代社会当中许多科学和工程技术等实际问题都可以转换归纳成为非对称线性方程组的特征值和特征向量求解问题,也被称为非对称性方程组求解问题,例如医学上的数据分析问题、图像以及信号等问题的处理、化学反应问题等,最后都归结为解决矩阵的特征问题,通过求解大规模的非对称线性方程组
Rm×n(Cm×n)αi=λiαi
(1)
来解决实际问题.其中,Rn(Cn) 表示n维实(复)向量空间,Rm×n(Cm×n)表示m×n实(复)矩阵的全体,(λi,αi)i=1,2,...,n为Rm×n(Cm×n)的特征对( ‖αi‖=1).因此研究非对称线性方程组的特征求解问题的数值方法具有非常重要的作用,关于探讨算法设计的好坏,开发出更有效的算法也成为一个非常重要的问题.
为了求出更优,精度更高的特征对,需要对算法的精度和实用性进行反复地研究.一个好的算法不仅需要具备数值收敛、占用的存储空间小、运算量少,还需要易于计算机实现.传统方式解决此类问题应用Lanczos 算法,当原始矩阵规模庞大时,利用这种方法就存在着计算复杂,存储量过大,实际问题中我们产生的矩阵都是上万阶,甚至于几十万阶,因此它的速度和内存问题就会格外突出.Lanczos 算法在求解时采用的是正交投影方法[1],虽然只需要存储4个向量,但是其特征值的精度不够高.除此之外,随着算法生成的Krylov子空间维度的增大,向量会逐渐失去正交性,直接导致算法的收敛速度和收敛性降低.为了解决这一系列的问题,对原始的Lanczos 算法进行改进和精化,使得求解出的特征对精度更高,速度更快.改变算法的解题方式,从正交投影的方式转变成为双正交的方式,由于矩阵空间维度扩大会影响向量的双正交性,所以借助压缩技术,在得到高精度特征值的同时维持其双正交性,且能减少硬件的内存消耗[2].
首先对非线性方程组进行初始化处理,将方程组转换成等价的三角矩阵.非线性方程组的矩阵形式为
(2)
其中,大写字母V,W,Y表示矩阵,小写字母v,w表示向量.在三角化中,首先设定u为启动向量,对于i=1,2,...,n,计算式(3):
Wk+1(β1e1)=bAWk=Wk+1Yk+1(∶,k+1)
(3)
其中,β1是非负的且满足Wk+1∈Rk+1,k+1,Wk+1为非对称三角矩阵,Yk+1(∶k+1)表示的是将Wk+1的最后一列去除掉.通常情况下,要是想找到矩阵A的奇异三元组,则必须对Wk进行奇异值分解,接着用Wk的奇异值近似A的奇异值,将Wk的奇异值向量与Lanczos向量相结合并对A的奇异值向量进行逼近,则可进行三角化处理,设A为非对称线性方程的等价矩阵表现形式,且A∈Rm×n(Cm×n),b表示正向量且b∈Rm×n,k为整数.执行循环算法:
y=norm(b);v=b/y;g=A*v;x=v*g;g=g-x*v;T(1,1)=x;
fori=2:k;Y(i,i-1)=Y(i-1,i)=y;end
(4)
由多项式的连续性可以得出,位移点对于不需要奇异值的近似程度更好.用每一次计算出的那些不需要的Ritz值作为位移,即为准确位移.取那些不需要的近似奇异值(Ritz值)作为位移,由此得到重新启动方法.
(5)
(6)
通过算法计算出近似值与真实值进行对比,对比两种算法计算出的近似值与真实值之间的误差,从而得出精化度对比结果.误差估计以子空间迭代法为标准进行对比,计算是在IBM-PC机上进行的.两种算法都是采用的多个初始向量进行迭代计算,精化Lanczos算法对所有初始向量采用了平行迭代格式,且借助了逐个加入初始向量的迭代格式.从而在取相同截断值的情况下,精化Lanczos算法的初始向量实际参加迭代的次数更少.迭代次数越多则收敛越快,计算结果精度也就越高.针对误差值的对比分析进行数据统计并绘制对比曲线如图1所示.
图1 精细度对比曲线Fig.1 The curve of precision contrastion
结合图中曲线分析,受收敛作用的影响,两种算法的精度随着求解数量的增加而呈上升趋势.数量较少时精化Lanczos算法的精度较低,比Lanczos算法的精度低102,待数量增多后精化Lanczos算法的误差急剧上升,最后维持在102以上.
表1 Lanczos算法与精化Lanczos算法的收敛效果Tab.1 The convergence effect of Lanczos algorithm and refined Lanczos algorithm
从表1中可知,两种算法在运行过程中得到的收敛结果相似,但由于在算法运行当中的收敛效果主要体现在算法的运行速度上,由于精化算法自身具备收敛检验和重启操作,所以精化Lanczos算法计算相对更加简单快速.
时间消耗是检测一个算法质量和效率的重要参考依据,在保证准确率的情况下,算法计算的运行时间越短,该算法的质量也就越高,性能也就越强.由于精化Lanczos算法有收敛检验的功能,所以可以自身控制收敛速度,进而影响了计算的时间.针对Lanczos算法和精化Lanczos算法都可以顺利地进行大规模非对称线性方程组的求解,并得到一定误差范围内的准确结果,可以确保此两种算法为可运行的算法,即可进行计算时间的探究与分析,可以统计出两种算法计算方程组消耗的时间并整理成对比图,如图2所示.
从图2中结果来看,对于大规模非对称线性方程组的求解算例,两种Lanczos算法计算得到的参数估值均值出入不大,而从观测数据上可以看出,精化Lanczos算法的计算时间在每次计算时均小于Lanczos算法30秒的时间.从图中可以看出随着数据量的增加,利用Lanczos算法计算时所涉及系数矩阵的协因数阵的维度扩大,分解出的向量构造也逐渐复杂,而精化Lanczos算法提取了系数矩阵中重复出现的随机元素,进行收敛操作和检验,此时两种Lanczos算法的计算时间差会迅速增大,精化的Lanczos算法大大节省了时间消耗.因此,实验同样证明了对大量数据进行求解处理时,采用精化Lanczos算法可以明显地提高计算效率,减少计算时间.
图2 两种算法时间消耗对比图Fig.2 The comparison of time consumption between two algorithms
在精化Lanczos算法当中引用了隐式压缩重启技术[10],其操作分为两个步骤,一个是隐式压缩,该技术的特点是选择若干个位移,对Lanczos过程产生的低阶上矩阵进行分解,使得计算量和存储量大大减少了,从而加快了Lanczos算法的收敛速度.隐式压缩重启技术已广泛应用于特征值问题,例如隐式重启Lanczos算法、Lanczos压缩重启算法、循环收缩Lanczos算法、隐式重启块Lanczos算法、隐式重启全局Lanczos算法[10]等.该技术在全局的应用中首先令k为要求特征值的个数,i=k+p,考虑i=k+p步全局Lanczos算法过程,有
(7)
其中,e为单位向量,Is为最小奇异值的奇异向量.x为位移,对Wk+p做带位移的QR分解,有
Wk+p-xI=QR
(8)
其中,Q是正交矩阵,R是三角矩阵.经过p步带位移隐式重启之后,可以得出结果:
(9)
计算结果分为计算结果的准确性和残量值[5].设定CPU表示计算时间单位为秒,Iter表示重启次数,m为投影了空间的维数,允许误差界tol=10-5,全局算法的残量范数取res,初始矩阵或向量均为随机选取.调用Matlab R2010a中的eig命令求得矩阵A的4个最大特征值为
ϑ1=7.992 413 314 948 18
ϑ2=7.981 047 676 817 97
ϑ3=7.981 047 676 817 96
ϑ4=7.699 682 038 687 75
用精化Lanczos算法计算初始矩阵的4个最大特征值,使用两种解题算法进行计算得出计算结果如表2所示.
表2 两种算法计算结果Tab.2 The results of two methods
从计算结果中可以看出 Lanczos算法得出的结果为小数点后两位,而精化 Lanczos算法特征值为小数点后五位,与设定的最大特征值对比两种算法的计算结果均在准确范围之内,但精化的算法在结果的体现上更加准确,即精化 Lanczos算法的准确性更高.至于残量值根据残量定理:当精化向量收敛到原矩阵的特征值当中,则精化向量的残量收敛值为零.利用Ipsen引理可以得出精化向量和特征向量的残量值.
(10)
图3 残量对比结果图Fig.3 The figure of residual comparison
相对非对称线性方程组求解算法及应用而言,非对称线性方程组求解初始化策略的研究相对较少,通过提出精化Lanczos算法并与Lanczos算法的各项性能进行对比,可以充分地验证出精化算法的优势,从理论上说明精化Lanczos方法比Lanczos方法更优越,可以在实际问题的解决中得到广泛应用.