线性代数机算中的速度和精度问题.doc

肇节莁蚁膀肄虿蚀衿芀薅蚀羂肃薁虿膄莈蒇蚈袄膁莃蚇羆莆蚂蚆肈腿薈蚅膁莅蒄螅袀膇莀螄羃莃芆螃膅膆蚄螂袅蒁薀螁羇芄蒆螀聿蒀莂蝿膁节蚁蝿袁肅薇袈羃芁蒃袇肆肃荿袆螅艿莅袅羈肂蚄袄肀莇薀袃膂膀蒆袃袂莆莂袂羄膈蚀羁肇莄薆羀腿膇蒂罿衿莂蒈薆肁芅莄薅膃蒀蚃薄袃芃蕿薃羅葿蒅薂肇节莁蚁膀肄虿蚀衿芀薅蚀羂肃薁虿膄莈蒇蚈袄膁莃蚇羆莆蚂蚆肈腿薈蚅膁莅蒄螅袀膇莀螄羃莃芆螃膅膆蚄螂袅蒁薀螁羇芄蒆螀聿蒀莂蝿膁节蚁蝿袁肅薇袈羃芁蒃袇肆肃荿袆螅艿莅袅羈肂蚄袄肀莇薀袃膂膀蒆袃袂莆莂袂羄膈蚀羁肇莄薆羀腿膇蒂罿衿莂蒈薆肁芅莄薅膃蒀蚃薄袃芃蕿薃羅葿蒅薂肇节莁蚁膀肄虿蚀衿芀薅蚀羂肃薁虿膄莈蒇蚈袄膁莃蚇羆莆蚂蚆肈腿薈蚅膁莅蒄螅袀膇莀螄羃莃芆螃膅膆蚄螂袅蒁薀螁羇芄蒆螀聿蒀莂蝿膁节蚁蝿袁肅薇袈羃芁蒃袇肆肃荿袆螅艿莅袅羈肂蚄袄肀莇薀袃膂膀蒆袃袂莆莂袂羄膈蚀羁肇莄薆羀腿膇蒂罿衿莂蒈薆肁芅莄薅膃蒀蚃薄袃芃蕿薃羅葿蒅薂肇节莁蚁膀肄虿蚀衿芀薅蚀羂肃薁虿膄莈蒇蚈袄膁莃蚇羆莆蚂蚆肈腿薈蚅膁莅蒄 线性代数机算中的速度和精度问题在用笔算时，通常都用整数矩阵来演示和做题，几乎不涉及误差，也就没有精度问题。至于矩阵计算速度的低下，即使在二、三阶的低价矩阵中，就已暴露无遗，但那从来不是纯粹进行理论探讨的数学家所关心的内容，甚至很怕读者触及这个敏感问题，这是用笔算解矩阵方程的根本弱点。承认这个弱点必然导致计算机的引入和课程的改造。线性代数进入应用领域，特别是使用计算机后，速度和精度两个现实问题就不可避免地摆在我们面前。此外，还有一个计算稳定性（或可靠性）的问题。一 精度问题1双精度浮点数的表示方法及其精度要理解科学计算中的精度问题，必须首先建立一个概念：计算机对数的表述和运算都是有误差的。对于MATLAB而言，它用二进制双精度浮点格式来表达一个数。按照IEEE标准，表达一个数需要8个字节，也就是64个二进制位。双精度浮点数用下式表示 其中M是一个小于一大于1/2的二进制分数，称为尾数，占用52个二进制位表示；而指数E是一个带符号的二进制整数。占11个二进制位，总共可表示2048个整数，即可以表示从-1023到+1024的数集，它决定了数的动态范围。数的正负号反映在S上，它只占一位。总计1+52+11=64位。浮点数的量化步长可以代表它的相对精度。它是由M的位数决定的。52位二进制数的量化步长是2 52=2.220410 16。该数的动态范围则取决于指数部分。因为2 102310 307及2 +102410+308。所以MATLAB中数的动态范围为2.225110 -3081.797710+308。在MATLAB命令窗中，键入eps, realmin, realmax，系统就会给出上面的几个数。在做浮点运算时，除非出现大数相减，相对误差可以基本不变。由于多次运算造成误差的积累，MATLAB中的大部分运算结果的误差比eps 要大一些。因此正常情况下（即系统不给出警告时），计算结果至少可以有12位十进制有效数字的可信度。但是如果遇到两个很接近的大数相减，把有效数位中的前N位都消掉了，那么数字的精度就减少了N位。比如12345-12344=1，前两个数可能具有5位有效数位，但它们的差只有一位有效。2高斯消元法中的精度问题（主元的部分交换partial Pivoting）举一个简单的例子，设方程组如下，写成矩阵形式：用消元法化简增广矩阵C=A,B由此得知 x=10000*10000/10001=0.9999; y=10000/10001=0.9999，这是手算的精确解。假如我们采用的是一个很粗糙的计算机，只能保持三位有效数字，那么它遇到10001时，只会解读为10000。则上面的消元过程将成为：这个结果将解读为x=0, y=1，x的误差达到了100%，完全不能采用。如果计算机的精度是N位，其结果虽然不至于那么惨，但也将使有效数位减少4位。现在来看问题出在什么地方？有什么办法可以避免？从计算过程看，其实问题就出在回代过程中出现了相接近的两个大数相减。两个大数出现于把微小的主元0.0001化为1的消元过程中，可见消元过程中主元太小是造成误差的根源。如果在做消元法的时候，检查主元行中各个元素，把最大的主元通过列交换放到主元的位置，那样这个问题就可以解决，例如上题中先交换第一行中的两列，可得：由此得知 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=AB， 得到X= 0.9999；0.99993由矩阵奇异性造成的误差方程左端的系数矩阵A的行列式是否等于零决定了此方程有无解。所以许多人以为det(A)是否接近于零将决定解的误差大小。例如对这一方程进行消元计算的过程将如下：结果是x=-0.0001,y=1。且x与y的有效数位都减小到一位。造成这种误差的原因是方程的奇异性，也就是行列式接近于零。在本题中的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，可得到AX=b，这些黑体的字符都代表向量或矩阵，它们的模（或范数）应该满足许伐茨不等式，即，于是有：将此两不等式左右端分别相乘，得到：其中就是条件数，它表示了解的相对误差与引起它的系数的相对误差之比，条件数愈大，解的误差受系数误差的影响愈大。也表示了该矩阵的奇异程度比较大，计算的稳定性和可靠性较低。这里还有一个问题不明确，那就是向量和矩阵的范数如何定义的。原则上，如果AX=b，则矩阵和向量范数的定义必须满足许伐茨不等式，能够满足这个条件的范数的定义是很多的，这里只介绍MATLAB中所采用的。向量的范数通常就是它的长度，即向量中各个元素平方和的开方，在MATLAB中用norm(b,2)表示。这是一个最通用的定义，所以可用norm(b)默认，也还有其他的定义。比如norm(b,1)则为b中各元素的绝对值之和。norm(b,inf)则取b中各元素的最大绝对值，矩阵的范数存在着更多定义，在MATLAB中就有五、六种。这里只介绍用得最多的，也是默认的norm(A)= norm(A,2)的定义。对于具有实特征根的正定方阵，它指的是A的最大特征根的绝对值。为了避免矩阵的复杂运算，举一个简单的二阶对角矩阵为例，设:对角矩阵的范数就是其绝对值最大的特征值，所以有:而根据同样的道理，由此可知，条件数。可见条件数c也就是A的最大特征值与最小特征值之比。条件数就作为衡量矩阵的解对系数b误差敏感程度的定量标志。同样可以证明：，即矩阵的解X对矩阵系数A误差敏感程度也可以用条件数来标志。注：对于长方阵，且特征值可能有复数时，它指的是方阵(AT*A)的最大特征值绝对值的平方根，后一个定义概括了前者。这个问题的研究通常归结为对方阵(AT*A)特征根的研究，通常也称为奇异值分解，英文为SVD(Singular Value Decomposition)分解，MATLAB中的svd函数具有 以下的调用格式：u,s,v=svd(A)，它把矩阵A分解为三个矩阵的乘积u*S*v，其中u和v分别是m*m及n*n维的正交矩阵，而S则是一个对角矩阵，其元素是按绝对值降秩排列的ATA的特征值。norm(A,2) =max(svd(A) =(ATA)max则取出其中的最大值。例如A=1,1.0001;1,1;键入Nc=cond(A)，得到 Nc= 4.0002e+004键入u,s,v=svd(A)，得到这样算出的条件数成了无穷大，这是由于显示精度太低造成的，如果用format long，可得可求出：条件数Nc=s(1,1)/s(2,2) = 4.000200007503750e+004条件数中10的指数部分说明该方程的解的有效数要降低的位数，这里是四位，与前分析同。当MATLAB在解线性方程组时发现其矩阵的条件数达到1016以上时，它会发出警告提示，说明所求出的解可能不对，因为那时MATLAB的有效数位16已被用尽，请使用者警惕。在纯理论的线性代数中，行列式的特性是二进制式的，要就是零，要就是非零，非黑即白；。行列式在零附近的取值和变化对方程组的解会产生何种影响是非常重要的问题。通常用Hilbert矩阵来构成逐次接近于奇异的矩阵。例如键入A=hilb(5)，便有：A = 1 1/2 1/3 1/4 1/5 1/2 1/3 1/4 1/5 1/6 1/3 1/4 1/5 1/6 1/7 1/4 1/5 1/6 1/7 1/8 1/5 1/6 1/7 1/8 1/9 要检验矩阵逐次奇异化时条件数的变化情况，可键入：for n=1:15 n,A=hilb(n);det(A),cond(A),x=inv(A)*ones(n,1);endMATLAB给出：n = 1，det(A) = 1, cond(A) = 1,n = 2 det(A) = 0.0833 cond(A)= 19.2815n = 3 det(A) = 4.6296e-004 cond(A)= 524.0568n = 4 det(A) = 1.6534e-007 cond(A)= 1.5514e+004 n = 5 det(A) = 3.7493e-012 cond(A)= 4.7661e+005 n = 6 det(A) = 5.3673e-018 cond(A)= 1.4951e+007n = 7 det(A) = 4.8358e-025 cond(A)= 4.7537e+008n = 8 det(A) = 2.7371e-033 cond(A)= 1.5258e+010n = 9 det(A) = 9.7203e-043 cond(A)= 4.9315e+011n = 10 det(A) = 2.1644e-053 cond(A)= 1.6025e+013n = 11 det(A) =3.0273e-065 cond(A)= 5.2239e+014n = 12 det(A) = 2.8581e-078 cond(A)= 1.7945e+016Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.632766e-017.n = 13 det(A) = 4.4480e-092 cond(A)= 3.7544e+018Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.339949e-018.n = 14 det(A) = -3.9220e-107 cond(A)= 4.0746e+017Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.708191e-019.n = 15 det(A) = -2.1903e-120 cond(A)= 8.4880e+017Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.543404e-018.RCOND是条件数cond(A)的反演，称为逆条件数，大体上相当于条件数的倒数1/Nc，RCOND愈接近于1，矩阵计算愈稳定， RCOND愈小，解就愈不准，RCOND接近于eps时，解的精度就没有意义了。从上例可看出当条件数达到1016，或RCOND小于10-16时，MATLAB将发出警告。4关于接近零的程度的讨论：容差tol对线性代数计算结果的影响理想的数学零在实际工程计算过程中几乎是永远不会出现的，或者说，数学零的出现概率为零。而线性代数中的许多函数却是要靠零来定义的，比如行列式是否为零决定矩阵的奇异性，秩是由行列式不为零的最大的子矩阵的阶数决定的，向量的相关性、正交性也都与运算结果是否等于零有关，等等。为了能把工程问题与数学理论相衔接，在数学软件中设立了容差量（通常用tol表示tolerance）。人们可以根据工程问题的性质，给出tol的值，意思指所有绝对值小于tol的数值，都当做零来处理。现在拿A=hilb(8)作为例子，看看它的秩rank(A,tol)如何随容差不同而变化。tol10-1010-810-610-510-410-210-1100rank(A,tol)87654321行最简型函数的结果也与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)得到：U = 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 1.0000 0 0.0833 0.0889 0.0833 0.0762 0.0694 0.0635 0.0583 0.6667 0 0 0.0065 0.0110 0.0139 0.0156 0.0165 0.0170 0.4861 0 0 0 0.0011 0.0024 0.0034 0.0043 0.0049 0.2500 0 0 0 0 -0.0000 -0.0001 -0.0002 -0.0003 -0.0429 0 0 0 0 0 0 0 0.0000 0.0036 0 0 0 0 0 0 0 0 -0.0011 0 0 0 0 0 0 0 0 0.0001不管增广项，可以看出，如果把小于0.0001的数看做数学零，则U有5个非零行（只看主元是4个），也就意味着A的秩应为5(4)，依此类推。这个问题在介绍多个空间点是否位在一个平面上的例题中将有实际的应用。5用最小二乘法处理误差超定方程的解参阅1 8.2节或2 5.7节。二 速度问题运算时间通常用完成某一运算所需做的乘加法次数来判断，所以运算速度就与乘法次数成反比。消元法正向消元所需的乘法次数为，回代所需的次数少的多，为，在n很大时，通常就简单地按来估计。N=10时，为333次；n=25时为5208次。对现有的最新计算机（2008年测试的最高纪录是每秒1015次，测试所用的标准就是用LINPACK软件包解方程AX=b），这都可以在几十个微微秒中完成，不算什么。求行列式所需的时间与此相同。但如果不用行消元的方法，而按行列式的数序定义来求，同样是25阶的行列式，就要做3.7*1026次乘法，即使用上述最新的计算机，也要1011秒钟。一年是86400秒，就算105秒吧，那也要100万年。所以消元法是一个巨大的成就，虽然以高斯命名，其实是中国人在1000多年前的“九章算法”中最先提出的方法，可惜只留下了算法，没有留下作者的姓名。矩阵的构