线性代数机算的速度和精度

线性代数机算的速度和精度机算中精度和速度的重要性•在用笔算时，通常都用整数矩阵来演示和做题，几乎不涉及误差，也就没有精度问题。至于矩阵计算速度的低下，即使在二、三阶的低价矩阵中，就已暴露无遗，但那从来不是纯粹进行理论探讨的数学家所关心的内容，甚至很怕读者触及这个敏感问题，这是用笔算解矩阵方程的根本弱点。承认这个弱点必然导致计算机的引入和课程的改造。•线性代数进入应用领域，必定要使用计算机，速度和精度两个现实问题就不可避免地摆在我们面前。在这里，我们将注重精度问题，并只着重于MATLAB命令中包含的、或会对计算精度作出提示的问题做一介绍。1．双精度浮点数的精度•按照IEEE标准，表达一个数需要8个字节，也就是64个二进制位。双精度浮点数η用下式表示••其中M是一个小于一大于1/2的二进制分数，称为尾数，占用52个二进制位表示；••而指数E是一个带符号的二进制整数。占11个二进制位，总共可表示2048个整数，即可以表示从-1023到+1024的数集，它决定了数的动态范围。数的正负号反映在S上，它只占一位。总计1+52+11=64位。ESM2)1(12/1MMATLAB中数的精度•浮点数的量化步长可以代表它的相对精度。它是由M的位数决定的。52位二进制数的量化步长是2-52=2.2204×10–16。该数的动态范围取决于指数部分。因为2–1023≈10–307及21024≈10308。所以MATLAB中数的动态范围为2.2251×10-308～1.7977×10+308。在MATLAB命令窗中，键入eps,realmin,realmax，系统就会给出上面的几个数。MATLAB中运算的精度•在做浮点加法、乘法、除法运算时，相对误差可以基本不变。由于多次运算造成误差的积累，MATLAB中的大部分运算结果的误差比eps要大一些。因此正常情况下（即系统不给出警告时），计算结果至少可以有12位十进制有效数字的可信度。但是减法有时会有问题，如果遇到两个很接近的大数相减，把有效数位中的前N位都消掉了，那么数字的精度就减少了N位。比如12345-12344=1，这两个数具有5位有效数位，但它们的差只有一位有效。2．高斯消元法中的精度问题——主元交换法（partialPivoting）•设方程组如下，写成矩阵形式：•用消元法化简增广矩阵C=[A,B]•由此得知x=10000/10001=0.9999;y=10000/10001=0.9999，这是手算的精确解。0.000110.00011x1*0-110xyxyyAX=B0.0001111100001000011000010000-110-110010001100001100001000010(1-1/10001)0110000/100010110000/10001C=[A,B]计算机舍入误差造成消元误差•假如我们采用的是一个很粗糙的计算机，只能保持三位有效数字，那么它遇到10001时，只会四舍五入解读为10000。则上面的消元过程将成为：•这个结果将解读为x=0,y=1，x的误差达到了100%，完全不能采用。如果计算机的精度是N位，其结果虽然不至于完全无用，但也将使有效数位减少4位。0.0001111100001000011000010000-110-1100100001000011000010000100011011C=[A,B]主元太小是造成误差的根源•从计算过程看，其实问题就出在回代过程中出现了相接近的两个大数相减。两个大数出现于把微小的主元0.0001化为1的消元过程中，主元太小是造成误差的根源。如果在做消元法的时候，检查主元行中各个元素，把最大的主元通过列交换放到主元的位置，那样这个问题就可以解决，例如上题中先交换第一行中的两列，可得：0.00011110.0001110.00011-1101-100-1.0001-110.000111011/10001011/1.0001011/1.0001C=[A,B]在rref函数中的换主元措施•由此得知x=10000/10001=0.9999;y=0.9999，这和精确解完全一致，说明只要进行主元最大的列交换。计算精度是可以保证的。•在rref命令中，语句[p,k]=max(abs(A(i:m,j)));就是为了在第i行中，找到绝对值最大的元素A(i,k)，以后再进行列交换。•在MATLAB命令窗中键入：•A=[0.0001,1;-1,1],B=[1;0],X=A\B，•得到•X=[0.9999；0.9999]键入typerref得出的主要程序代码•while(i=m)&(j=n)•%Findvalueandindexoflargestelementintheremainderofcolumnj.•[p,k]=max(abs(A(i:m,j)));k=k+i-1;•if(p=tol)•%Thecolumnisnegligible,zeroitout.•A(i:m,j)=zeros(m-i+1,1);•j=j+1;•else•%Remembercolumnindex•jb=[jbj];•%Swapi-thandk-throws.•A([ik],j:n)=A([ki],j:n);•%Dividethepivotrowbythepivotelement.•A(i,j:n)=A(i,j:n)/A(i,j);•%Subtractmultiplesofthepivotrowfromalltheotherrows.•fork=[1:i-1i+1:m]•A(k,j:n)=A(k,j:n)-A(k,j)*A(i,j:n);•end•i=i+1;•j=j+1;•end•end3．由矩阵奇异性造成的误差方程左端的系数矩阵A的行列式是否等于零决定了此方程有无解。所以通常的概念是行列式是否接近于零将决定解的误差大小。例如•对这一方程进行消元计算的过程将如下：•结果x=-0.0001,y=1。且x与y的有效数位都减小到一位。因为在解题过程中出现过大数相减的问题。1.0001111.0001x1*0.9999110.9999xyxyyAX=B11.0001111.0001110-0.0001110.99990-0.0001-0.0001011C=[A,B]要用误差的相对值来评价•造成这种误差的原因是方程的奇异性，即系数行列式接近于零。在本题中的det(A)=-0.0001，引起解的精度下降4位有效数。但这种用绝对误差为标准的方法不科学，若把A的所有元素都加1000倍，方程性质没变，行列式却大了1000倍，显得不那么奇异了，所以行列式绝对值的大小不能严格地说明奇异性。•需要仔细地分析矩阵方程AX=b中，系数矩阵b或A的相对误差Δb/b或ΔA/A会引起解X多大相对变化ΔX/X，从而判定A的奇异程度。一般地分析这个问题相当困难，但可以由简单的二维和对角矩阵实例中，找到一些基本规律。相对误差的分析•先考虑系数b的误差Δb，由AX=b及A(X+ΔX)=b+Δb，可得到AΔX=Δb，这些黑体的字符都代表向量或矩阵，它们的模（或范数）应该满足许伐茨不等式，即，于是有：•将此两不等式左右端分别相乘，得到：,-1ΔXAΔbbA*X**c-1-1ΔXbA*XAΔbΔXΔbΔbA*AXbbx*yx*y条件数的概念•其中•就是条件数，它表示了解的相对误差与引起它的系数的相对误差之比，条件数愈大，解的误差受系数误差的影响愈大。也表示了该矩阵的奇异程度比较大，计算的稳定性和可靠性较低。•正定方阵的范数定义为其最大特征值的绝对值，故||A||是A的最大特征值，||A-1||是A的最小特征值的倒数，所以c等于最大与最小特征值之比。•二维矩阵条件数的概念也可以从几何图形上去理解：例如，在eigshow([8,3;-2,1])图形中：c-1A*A一个数字例•举一个简单的二阶对角矩阵为例，设:•对角矩阵的范数就是其绝对值最大的特征值，所以有||A||=6，而根据同样的道理||A-1||=1/3，由此可知，条件数c=||A||*||A-1||=6/3=2。可见条件数c也就是A的最大特征值与最小特征值之比。它是衡量矩阵的解对系数b误差敏感程度的定量标志。•例如A=[1,1.0001;1,1];•键入[u,s,v]=svd(A)，•键入Nc=cond(A)，得到Nc=4.0002e+00461/61/3-100A=,A=,030则用Hilbert矩阵检验条件数•条件数中10的指数部分说明该方程的解的有效数要降低的位数，这里是四位，与前分析同。当MATLAB在解线性方程组时发现其矩阵的条件数达到1016以上时，它会发出警告提示，说明所求出的解可能不对，因为那时MATLAB的有效数位16已被用尽，请使用者警惕。•通常用Hilbert矩阵来构成奇异程度逐次增加的矩阵•A=hilb(5)：•要检验矩阵逐次奇异化时条件数的变化情况，可键入：•forn=1:15•n,A=hilb(n);det(A),cond(A),•x=inv(A)*ones(n,1);•end用Hilbert矩阵检验的结果•程序运行后的屏幕显示•….•n=15det(A)=-2.1903e-120cond(A)=8.4880e+017•Warning:Matrixisclosetosingularorbadlyscaled.•Resultsmaybeinaccurate.RCOND=1.543404e-018.•RCOND是条件数cond(A)的反演，称为逆条件数，大体上相当于条件数的倒数1/Nc，RCOND愈接近于1，矩阵计算愈稳定，RCOND愈小，解就愈不准，RCOND接近于eps时，解的精度就没有意义了。从上例可看出当条件数达到10^16，或RCOND小于10^-16时，MATLAB将发出警告。4．关于接近零的程度的讨论：容差tol对线性代数计算结果的影响•理想的数学零在实际工程计算过程中几乎是永远不会出现的，或者说，数学零的出现概率为零。而线性代数中的许多函数却是要靠零来定义的，比如行列式是否为零决定矩阵的奇异性，秩是由行列式不为零的最大的子矩阵的阶数决定的，向量的相关性、正交性也都与运算结果是否等于零有关，…等等。为了能把工程问题与数学理论相衔接，在数学软件中设立了容差量（通常用tol表示tolerance）。人们可以根据工程问题的性质，给出tol的值，意思指所有绝对值小于tol的数值，都当做零来处理。tol对rank函数的影响tol10^-1010^-810^-610^-510^-410^-210^-110^0rank(A,tol)87654321A=hilb(8)的秩rank(A,tol)随容差不同而变化。行最简型函数的结果也与tol有关，为了说明这点，用A=hilb(8)作为矩阵的系数构成方程A*X=ones(8,1),故增广矩阵为：C=[hilb(8),ones(8,1)]，下面求C的行最简型。为了看出行简化中元素绝对值的大小，我们将不执行把主元化为1的行倍乘操作。这样的命令在MATLAB中没有，只存在于教材《工程线性代数（MATLAB版）》的程序集dsk06中。其中ref1是前向消元的结果，ref2则是回代后的结果。键入：U=ref1([hilb(8),ones(8,1)])tol对rref函数的影响U=1.00.5000.3330.2500.2000.1670.1430.12500.0830.0890.0830.0760.0690.0640.058000.0060.0110.0140.0160.