group 1 _4_

实验报告

——第四次上机实验 group 1

实验内容：

1、 对比试验借助法方程组，QR 因子分解和SVD 方法解决最小二乘问题。

2、 要求分析合理，结合算例，注重文档质量。

实验原理：

自从高斯和勒让德在1800年左右发明最小二乘数据拟合以来，它的应用遍及数学科学的各个分支。最小二乘问题即如何求解一个矩阵的超定方程组Ax=b;最小二乘的思想是通过最小化剩余向量的2-范数来求解方程组的。

最小二乘问题的方法通常有三种：法方程组法、QR 因子分解法、SVD 法。 ① 法方程组法.

它是解最小二乘问题的经典方法。剩余向量r= b=Ax,且r⊥range(A),则A*r=0

即 A*Ax=A*b，若A 满秩，A*A 是正定的Hermite 矩阵，可以用楚列斯基因子分解法进行简化，得到A*A=R*R,其中R 是上三角矩阵，则方程组化为R*R=A*b。接下来利用下算法来解方程组。

算法11.1

1、 形成矩阵A*A 和向量A*b。

2、 计算楚列斯基因子分解A*A=R*R。

3、 对w 解下三角方程组R*w=A*b。

4、 对x 解上三角方程组Rx=w。 总工作量：≈3

23

1n mn +

flop ②QR 因子分解

QR 因子分解法是自20世纪60年代开始流行的解最小二乘问题的“现代经典”方法，是基于约化QR 因子分解的。有格拉姆?施密特正交化或更更常用的豪斯霍尔德三角形化可以构造

因子分解A=R Q ??。正交投影算子P 可写成P=*??Q Q ，这样与有y=Pb=*??Q Q b。因为y∈range(A),所以方程组Ax=b,有一个准确解，组合QR 因子分解得到，b Q x R

*?=?，为一上三角方程组，若A 满秩，则它是非奇异的，容易有向后带入求解。 算法11.2

1、 计算约化QR 因子分解A=R Q ??。

2、 计算向量b Q

*?。 3、 对x 解上三角方程组b Q x R

*?=?。 总工作量：≈3

2

3

2?2n mn flop ③SVD 法

P 可以表示为P=?U U ??，则y=Pb=b U U

???,可以得到b U U x V U ?=?Σ?????,及b U x V ?=?Σ?? 算法11.3

1、 计算约化SVD ?Σ=V U A ??。

2、 计算向量b U

??。 3、 对w 解对角方程组b U w ?=Σ

??。 4、 令x=Vw。

总工作量：≈3

2

11+2n mn flop

算法的比较

从速度方面考虑时，法方程组法是最好的。但四由于舍入误差的出现，解法方程组不总是稳定的。从而QR 方法是经常使用的。然而当A 接近秩亏损时，QR 方法不稳定，在这种情况下，一般使用SVD 法。

实验1

在解决最小二乘问题的方法（法方程组、豪斯霍尔德三角化、格拉姆-斯密特正交化法、SVD 法）中，验证法方程组法一般是不稳定的。 程序见stability.m

实验思路：

构造100×15的范德蒙德矩阵A，其列是幂1，τ，……，14

τ

在t 的点上的样本，右边的b

是函数))exp(sin(τ4在这些点上的样本。并且对b 进行了正规化，由于之前得到的x（15）=2.006787914503308e +003，正规化之后得到的x(15)=1,使实验容易进行。

并用下程序来计算A 的稳定性： X =A \b ; y=A*x ;

K appa=c on d (A)

Th eta=asin(norm(b-y)/norm(b)) E ta=norm(A)*norm(x)/norm(y) 结果如下：

K appa =2.271777594999816e +010 Th eta =3.746111141881315e-006 E ta =2.103559861988276e +005

得到y 和x 对b 和A 扰动的条件数： y

x

b 1.00000000000702 1.079968122641626e +005 A

2.271777595015756e +010

3.190866253243956e +010

运行程序得到如下结果： h o u se h ol d er 三角化

x(15)=1.00000025012134 增广矩阵法

x(15)=1.00000025010651 gs 正交化

x(15)= 0.97625233534956 增广矩阵法

x(15)=1.00000014033264 法方程组

x(15)=-0.72375775186580 SVD 法

x(15)=1.00000025012316

从上述结果看精度大约有710的误差，这是由于机机ε被A 的条件数放大所致，并非算法的不稳定所致。豪斯霍尔德法及其增广矩阵法得到的结果大致相同，且两者算法稳定。格拉姆-斯密特正交化法的结果误差较大，但其增广矩阵法能得到较好的解，其算法稳定。而法方程组法得到的解极不稳定，由定理19.3可以知道满秩最小二乘法的法方程组法在什么条件下是稳定的。SVD 法算法是稳定的。

实验2

用法方程组、QR 法、SVD 法解最小二乘问题Ax=b（A∈n m C ×,b∈m C ）

⑴当A是满秩矩阵且m>n,对于最小二乘问题Ax=b求解（见c ompare）

A = 1 1 1 1

1 2 3 4

1 3 6 10

1 4 10 20

1 5 15 35

b = [ 1 -2 3 0 -1];

法方程组法

x =7.20000000000049 -15.00000000000097 11.42857142857212 -3.00000000000017

QR因子分解

x =7.19999999999996 -14.99999999999996 11.42857142857140 -2.99999999999999 SVD方法

x =7.19999999999977 -14.99999999999957 11.42857142857110 -2.99999999999992

c on d(A)=5.630271651605274e+002

则三种算法得到的解均稳定,且得到的解大致一样，三种方法都是用满秩的最小二乘问题的

求解，且法方程组法最快。

exp(sin(τ4的样本点值，求Ax=b的最小二乘问题⑵构造100×15的范德蒙德矩阵A,b是))

的解。(见compare2.m)

经计算得

c on d(A)=2.271777594999816e+010

ran k(A)=15

法方程组法

x =1.0e+004 *

0.00009998584565 0.00040403820475 0.00060778935205 0.00365997278261 -0.039304465401670.19995747736018-0.72094730187206 1.51472500770520 -1.586951617325070.22023517182503 1.43554520939307-1.78520689013300 1.01612622879876 -0.29377513818611 0.03487145186049

QR因子分解

x =1.0e+005 *

0.00001000012035 0.00003993699213 0.00008441307716 -0.00012209295732 0.00145351255131 -0.01613043466124 0.08654970525716 -0.32446021338388 0.83987446597834 -1.45561207128515 1.67627615374629-1.26681558769009 0.60410162054904 -0.16483602611027 0.01959130937631

SVD方法

x=1.0e+005 *

0.000010000126280.00003993475785 0.00008454335416 -0.00012516266370 0.00149202153511 -0.016423944776480.08801381252476 -0.32945470085255 0.85181468145895 -1.47580780609404 1.70032727562214 -1.28653905488980 0.61470258459311 -0.16819737268221 0.02006787955022

可以看出当A的条件数很大时，法方程组法不稳定。这时应使用QR分解法或SVD法。

⑶当A是秩亏损时，法方程组法和QR法均不稳定。（见程序c ompare2.m）

A = [ 35 1 6 26 19 24;

3 32 7 21 23 25;

31 9 2 22 27 20;

8 28 33 17 10 15;

30 5 34 12 14 16;

4 36 29 13 18 11]

b=[1 -2 3 0 -1 0]'

ran k（A）=5

其中，

λ(A)=1.0e+002 *

1.11000000000000 0.27000000000000 -0.27000000000000 0.09797958971133 -0.09797958971133

ran k(A)<6 则A是秩亏损矩阵。

法方程组法

x =

1.65554978054980 1.59999422499424 -0.82552263802265 -1.27754677754679

-1.51365788865791 0.37019230769231

QR因子分解

x =1.0e+013 *

-7.39811027083440 -7.39811027083440 3.69905513541720 7.39811027083443

7.39811027083442 -3.69905513541725

SVD方法

X=1.0e+014 *

-1.02375488540363 -1.02375488540364 0.51187744270182 1.02375488540364

1.02375488540364 -0.51187744270182

可以看出，前两种方法不稳定，此时应用SVD解法来求解问题。