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

线性代数机算的速度和精度 机算中精度和速度的重要性 在用笔算时 通常都用整数矩阵来演示和做题 几乎不涉及误差 也就没有精度问题 至于矩阵计算速度的低下 即使在二 三阶的低价矩阵中 就已暴露无遗 但那从来不是纯粹进行理论探讨的数学家所关心的内容 甚至很怕读者触及这个敏感问题 这是用笔算解矩阵方程的根本弱点 承认这个弱点必然导致计算机的引入和课程的改造 线性代数进入应用领域 必定要使用计算机 速度和精度两个现实问题就不可避免地摆在我们面前 在这里 我们将注重精度问题 并只着重于MATLAB命令中包含的 或会对计算精度作出提示的问题做一介绍 1 双精度浮点数的精度 按照IEEE标准 表达一个数需要8个字节 也就是64个二进制位 双精度浮点数 用下式表示其中M是一个小于一大于1 2的二进制分数 称为尾数 占用52个二进制位表示 而指数E是一个带符号的二进制整数 占11个二进制位 总共可表示2048个整数 即可以表示从 1023到 1024的数集 它决定了数的动态范围 数的正负号反映在S上 它只占一位 总计1 52 11 64位 MATLAB中数的精度 浮点数的量化步长可以代表它的相对精度 它是由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 这是手算的精确解 计算机舍入误差造成消元误差 假如我们采用的是一个很粗糙的计算机 只能保持三位有效数字 那么它遇到10001时 只会四舍五入解读为10000 则上面的消元过程将成为 这个结果将解读为x 0 y 1 x的误差达到了100 完全不能采用 如果计算机的精度是N位 其结果虽然不至于完全无用 但也将使有效数位减少4位 主元太小是造成误差的根源 从计算过程看 其实问题就出在回代过程中出现了相接近的两个大数相减 两个大数出现于把微小的主元0 0001化为1的消元过程中 主元太小是造成误差的根源 如果在做消元法的时候 检查主元行中各个元素 把最大的主元通过列交换放到主元的位置 那样这个问题就可以解决 例如上题中先交换第一行中的两列 可得 在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 endend 3 由矩阵奇异性造成的误差 方程左端的系数矩阵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 可得到A X b 这些黑体的字符都代表向量或矩阵 它们的模 或范数 应该满足许伐茨不等式 即 于是有 将此两不等式左右端分别相乘 得到 条件数的概念 其中就是条件数 它表示了解的相对误差与引起它的系数的相对误差之比 条件数愈大 解的误差受系数误差的影响愈大 也表示了该矩阵的奇异程度比较大 计算的稳定性和可靠性较低 正定方阵的范数定义为其最大特征值的绝对值 故 A 是A的最大特征值 A 1 是A的最小特征值的倒数 所以c等于最大与最小特征值之比 二维矩阵条件数的概念也可以从几何图形上去理解 例如 在eigshow 8 3 2 1 图形中 一个数字例 举一个简单的二阶对角矩阵为例 设 对角矩阵的范数就是其绝对值最大的特征值 所以有 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 004 用Hilbert矩阵检验条件数 条件数中10的指数部分说明该方程的解的有效数要降低的位数 这里是四位 与前分析同 当MATLAB在解线性方程组时发现其矩阵的条件数达到1016以上时 它会发出警告提示 说明所求出的解可能不对 因为那时MATLAB的有效数位16已被用尽 请使用者警惕 通常用Hilbert矩阵来构成奇异程度逐次增加的矩阵A hilb 5 要检验矩阵逐次奇异化时条件数的变化情况 可键入 forn 1 15n A hilb n det A cond A x inv A ones n 1 end 用Hilbert矩阵检验的结果 程序运行后的屏幕显示 n 15det A 2 1903e 120cond A 8 4880e 017Warning 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函数的影响 A 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 0160 0170000 0010 0020 0030 0040 0050000 0 0 0001 0 0002 0 000300000000 00000000000000000000不管增广项 可以看出 如果把小于0 0001的数看做数学零 则U有5个非零行 只看主元是4个 也就意味着A的秩应为5 4 依此类推 这个问题在介绍多个空间点是否位在一个平面上的例题中将有实际的应用 5 用最小二乘法处理误差 超定方程的解 参阅 1 8 2节或 2 5 7节 工程问题都可以允许方程有误差 把一组解代入方程后 每个方程都有误差 要找误差在一定意义下的总和为最小的解 在这样的思路引导下 就产生了超定方程求解的方法 对于m n 即方程数多于变量数的情况 严格地满足Ax b的解x是不存在的 要找到近似的 就要引入误差e Ax b e 误差线性方程组的建立 引入误差向量e e Ax b写出其完全的矩阵形式如下问题是 找到解x 使e的长度或范数为最小 从向量空间的视点分析 研究例6 1的超定方程组 d 改写成简写为选择不同的x1和x2将得到不同的合成向量A x x1v1 x2v2 q q必定处于v1和v2张成的平面之内 而方程中的b则一般不会在这个平面内 所以精确解通常不可能存在 本例的向量空间图 这时最近似的解就应该是该平面上与b点最近的点所对应的坐标A xhat 它应该是b点向v1和v2张成的平面的投影 所以A xhat和b的连线应该和v1和v2张成的平面垂直 也就是说必须分别与v1和v2正交 如图8 6所示 最小二乘解的公式推导 A xhat和b的连线向量应该是这两个向量之差 即 它与v1和v2正交的要求可以分别表示为 和综合在一起可以写成 最后得到公式 最小二乘解的数字例 例8 求题6 1 d 方程组的最小二乘解 解 MATLAB程序ag808如下 A 1 1 1 1 1 2 b 1 3 3 xhat inv A A A be A xhat b norm e 运行此程序 得到 MATLAB中超定方程的解 在MATLAB中 把运算 ATA 1AT单独编成一个子程序 称为pinv函数 求最小二乘解的公式可以写成x pinv A b 与 适定方程 的解x inv A b非常相似 只是pinv函数并不要求A是方阵 最小二乘解也可用 运算符表示 这就把 欠定方程 适定方程 和 超定方程 用统一的运算格式 x A bMATLAB会自动根据系数矩阵A的行数m和列数n 来判断采用哪个方法和程序 不过对于欠定方程 这个式子只给出了一个特解 没给通解 二 机算速度问题 运算时间通常用完成某一运算所需做的乘加法次数来判断 消元法正向消元所需的乘法次数为 n3 n 3 回代所需的次数少 为n2 2 在n很大时 近似按n3 3估计 N 10时 为333次 n 25时为5208次 对现有的最新计算机 2008年测试的最高纪录是每秒1015次 测试所用的标准就是用LINPACK软件包解方程AX b 这都可以在几十个微微秒中完成 求行列式所需的时间与此相同 但如果按行列式的数序定义来求 25阶的行列式 就要做3 7 1026次乘法 用最新的计算机 也要100万年 所以消元法是一个巨大的成就 虽然以高斯命名 其实是中国人在1000多年前的 九章算法 中最先提出的方法 矩阵形式对计算速度影响 矩阵的构造形式对计算速度影响也很大 最显然的是对角矩阵 它完全可以按标量计算 算一个行列式只需要做n 1次乘法 解一