稳定性是数值分析的一个基本问题。
--L N. Trefethen
一个问题定义为由数据的向量空间 X 到解空间 Y 的一个函数 f:X->Y。相应地,一个算法可以看成是两个相同空间之间的另外一个映射 f{bar}:X->Y。注意,前者大部分情况下是一个连续系统,而由于计算机浮点数表示的原因后者是离散系统(即里面表示的数字是可数的,而针对浮点数而言,它不仅可数,而且是有限个数的)。离散系统要表达出连续系统必然要进行舍入。因而,f^{bar}的结果势必要受到舍入误差的影响。数值稳定性要解决的是一个算法,是否能够使用离散系统取得“正确答案”[1]。
显然,一个好的算法应该保证对于给定的 x,考虑计算的相对误差(||f(x)-f{bar}(x)||)/||f(x)||,自然地,我们期望相对误差很小,由于计算机浮点数精度的限制,它有个界限,不妨记作e_{mach},如果对每个x,有(||f(x)-f{bar}(x)||)/||f(x)|| = O(e_{mach}),我们就可以说算法 f{bar} 对问题 f 是准确的。进一步地,由于 f{bar} 的定义域是离散的,如果对于每个x,(||f(x)-f{bar}(x{bar})||)/||f(x)|| = O(e_{mach}),对某些满足||x-x{bar}||/||x{bar}||的x成立,则说算法 f{bar}是稳定的。另外,如果f{bar}(x{bar})=f(x)对于满足||x-x{bar}||/||x{bar}||=O(eps_{mach})的x成立,则说算法是向后稳定的。值得注意的是,有些算法是稳定的但不是向后稳定的,如计算sin(x)或cos(x)。给出这么多铺垫,涉及到的符号和术语也很多,连我自己都觉得有些绕,如果只是了解大概内容,上面的内容可以忽略,而直观的解释参见下面叙述的数值例子。
假设一个算法向后稳定,且关于问题的条件数较小的话,那么可以得到准确的结果。这个结论由定理保证:假设一个向后稳定的算法在具有条件数k的问题f:X->Y,则相对误差满足:||f{bar}(x)-f(x)||/||f(x)|| = O(k*eps_{mach})。
实践已经证明,数值代数的大多数算法而言,向前误差分析比起向后误差分析更难于实施(注1)。除此之外,很多时候,向后误差分析能更合理地反映出算法的误差。给出了一个数值例子[1],关于一个随机矩阵Q和R(注2),然后计算它们的乘积,不妨记作A,使用Householder三角化方法[2]对它进行QR分解,得到数值解Q1和R1,接着分析数值解和理想解之间的误差,结果表明数值的结果仅有2位或者3位有效数字(考虑相对误差),而实际计算过程中采用double型浮点数(即有16位有效数字)。而数值解的两个矩阵的乘积Q1*R1和矩阵A之间的误差能够达到15位的有效数字。另外,如果对Q和R进行微小扰动得到Q2和R2,两者的乘积Q2*R2和A之间的相对误差却仅仅能够达到了3位有效数字。从结果上来看,Q2(或R2)比Q1(或R1)更接近于Q(或R),但是乘积Q2*R2却比乘积Q1*R1更接近于A(Q1和R1中的误差被Wilkinson称为“魔鬼相关”的)。这个例子表明,使用向后误差分析比向前误差分析更合理。
关于向后稳定性的两点结论:
一,Householder三角形化以及相应的使用它来进行QR因子分解求解线性方程组Ax=b向后误差稳定的,前者用数学语言表达出来的形式为:令矩阵A的QR因子分解是由Householder三角形化进行数值计算得到因子Q1和R1,则有Q1*R1=A+dA, ||dA||/||A||=O(eps_{mach})对某个dA成立。
二,关于矩阵特征值有如下的Weyl定理:设 A 和 E 是n*n阶对称阵,设 alpha_1>=alpha_2>=...>=alpha_n 是 A 的特征值,而 alpha_{bar}_1>=alpha_{bar}_2>=...alpha_{bar}_n 是 A_{bar}=A+E 的特征值,则 |alpha_i - alpha_{bar}_i|<=||E||_2。虽然该定理是向前误差分析的结果,但是常常利用该定理可以确定QR迭代法(向后误差稳定)算法计算得到的数值特征值的误差界[3]。
注1:一般而言,大的向前误差可能是一个病态问题或者是一个不稳定算法的结果。
注2:作为一个法则,随机三角形矩阵的列空间序列作为矩阵元素的函数是极端病态的。
参考文献:
[1] 数值线性代数 Chap14-17,L N. Trefethen,David Bau, lll 著,陆金甫,关治译,人民邮电出版社,2006年
[2] 矩阵计算(第三版),Gene H.Golub,
[3] 应用数值线性代数 Chap5,J W. Demmel 著,王国荣译,人民邮电出版社,2007年