大数据计算方法 课件 chap3-方程求解与矩阵分解

1大数据计算方法Chap3：线性方程组求解与矩阵分解基于LU分解的方程求解条件数与舍入误差分析对称正定矩阵的Cholesky分解QR分解与正交三角化线性最小二乘问题解法Outline2基于LU分解的方程求解方程求解的高斯消去法部分主元LU分解的推导排列阵、上/下三角阵方程部分主元LU分解的程序有关Matlab命令3

求解线性方程组

4

高斯消去法5用高斯消去法求解单个线性方程组有时需要求解仅右端向量不同的多个方程组线性方程组求解器6常微分方程组的例子：用向后欧拉法(backwardEuler)求解线性电路瞬态仿真问题，假设时间步长为常数

t举例7,1

1F写出A,B固定不变反复用高斯消去法求解不同右端项的方程组效率低能否节省计算量、复用中间结果?LU分解8x1b1x2b2使用LU分解!bxyy只做1次反复做实际中，还要

考虑选主元部分主元(列主元)求解线性方程组的一个例子部分(列)主元的高斯消去过程

9

(第1个方程乘0.3,-0.5,加到下面的方程上)”乘子”:当前列其他系数除以主元，即-0.3,0.5互为相反数

部分(列选)主元的高斯消去过程(第1个方程乘0.3,-0.5,下加)消元时,乘子为-0.0410最后一个方程

L=

部分(列选)主元LU分解A矩阵A与L,U的关系：UPA=LU11M2P2M1A=U由乘子填充;选主元行交换称为部分主元LU分解行交换阵P2消去矩阵M1,M2M2”乘子”:第一次为-0.3,0.5；第二次为-0.04

排列阵P:由单位阵经一系列的行交换而得.每行(列)上有且仅有一个1，其余为0例:向量p表示单位阵各行的新顺序，简洁地表示排列阵I(p,:)就是排列阵P如何计算PA？求解系数为排列阵的方程>>B=A(p,:)排列阵与上(下)三角阵PA

~

对A作行重排AP~对A作列重排>>p=[4132];>>

I=eye(4);P=I(p,:)(P为正交阵)121:n的一个重排%时间短,省内存I(:,p)中p代表列的顺序因此就是PTpermutationmatrixx(p(1))=b(1)x(p(2))=b(2)...>>x(p)=b;

系数矩阵为上(下)三角阵的方程矩阵的LU分解得:上三角阵U，单位下三角阵L解Ux=b的两种”回代”(backsubs)排列阵与上(下)三角阵x=zeros(n,1);forj=n:-1:1x(j)=b(j)/U(j,j);

i=1:j-1;b(i)=b(i)-x(j)*U(i,j)endx=zeros(n,1);fork=n:-1:1

j=k+1:n;x(k)=(b(k)-U(k,j)*x(j))/U(k,k)end从b中依次减去U某列的倍数用U的行与解出的部分x作内积13b;x(这样,b不改变)L=LU分解的算法14算法3.1不选主元的LU分解

乘子原地存储：lupp程序“原地存储”方式最后两行算出L和Ufunction[L,U,p]=lupp(A)[n,n]=size(A);p=(1:n)';fork=1:n-1[r,m]=max(abs(A(k:n,k)));m=m+k-1;

if(A(m,k)~=0)

if(m~=k)A([km],:)=A([mk],:);p([km])=p([mk]);

end

i=k+1:n;A(i,k)=A(i,k)/A(k,k);j=k+1:n;A(i,j)=A(i,j)-A(i,k)*A(k,j);

endend最大值,位置索引交换两行算乘子更新未消去部分(kijversion)15L=tril(A,-1)+eye(n);U=triu(A);部分主元LU分解程序尽量用向量/矩阵运算保证算法不中断防止小主元导致的误差传播例子:在一个只有3位十进制精度的机器上求解如果不选主元,直接用高斯消去法为什么必须选主元?

16对上例,如果采用列主元高斯消去法为什么必须选主元?17算法的稳定性不同!非常接近准确解!其他选主元策略全主元：从未消去部分选最大矩阵元素，然后交换行、列行主元：从当前行选最大元素，然后交换列带阈值选主元：选某个大于最大元素倍(如=0.2)的元素Matlab中相关命令lu:

部分主元LU分解基于LU分解求解方程Ax=b:\(backwardslash):

X=A\B求解线性方程组AX=B基于部分主元LU分解的算法:mybslash.m是”\”的简化版(习题3.6)部分主元LU分解程序18稀疏阵稠密阵Lapack(BLAS)UMFPACK/CHOLMODPA=LU

LUx=PbSolveLy=PbSolve

Ux=yO(n3)O(n2)O(n2)A(准确度与运行时间的tradeoff)条件数与舍入误差分析复习矩阵范数敏感性分析与矩阵条件数考虑计算误差的分析算法稳定性19范数与矩阵条件数(复习)

(曼哈顿范数)(“最大”范数)(欧氏距离)

>>x=(1:4)/5>>norm1=norm(x,1)>>norm2=norm(x)>>norminf=norm(x,inf)20范数与矩阵条件数(复习)

01

21范数与矩阵条件数(复习)

条件数矩阵

条件数22F-范数:

(不是诱导范数)01

矩阵条件数23线性变换对

单元圆的扭曲

(排列阵)

(在p范数下)

算法稳定性与解的误差24

J.H.Wilkinson矩阵计算的舍入误差分析(绝对值)

(增长因子)

算法稳定性与解的误差

交换行,消元(采用截断舍入)解出:若不知道准确解,计算残差

25

注意：残差或相对残差不总反映结果的准确度某种相对残差相对残差小是算法稳定的必要条件

算法稳定性与解的误差26cond~106

矩阵条件数

算法稳定性再次说明:算法稳定性与问题敏感性共同影响解的误差

Cholesky分解对称正定矩阵的分解分解算法与应用27Cholesky分解定理1.如果实对称阵A的前n-1个顺序主子式

0

(LDLT分解)：其中L为单位下三角阵，D为对角阵2.若矩阵A同时正定，则存在非奇异下三角阵L，若限定L

的对角元>0,则此分解唯一In2x2case,forexample,whichimplies

28Cholesky分解算法可根据LU分解算法的kij版本进行修改，考虑对称性；原地存储:仅使用A的下三角部分，L的结果将其覆盖1.Fork=1ton1’.2.Fori=k+1ton3.4.Forj=k+1toi5.6.End7.End8.End注意：L不是单位

下三角=AFork=1tonFori=k+1tonEndEnd第2种算法29(不用选主元)每步计算都求得L中的某个元素!对称正定矩阵的Cholesky分解算法计算量分析(稠密阵)：仅需要n3/6次乘除法和差不多的加减法，是LU分解的一半计算量对于对称正定矩阵，可证明：所有n个平方根运算的操作数都为正数，算法可行Cholesky算法是稳定的(对舍入误差不敏感)，不需要选主元Matlab命令为chol,缺省得到上三角矩阵(R'*R=A),

可用于检测矩阵正定性解Ax=b:>>R=chol(A)>>L=chol(A,‘lower’)30(仅使用A的上三角部分元素)A=RTR

RTRx=bSolve

RTy=bSolveRx=y复杂度:O(n3)O(n2)O(n2)应用实例

000

31QR分解与正交三角化QR分解Gram-Schmidt正交化过程Householder变换Givens变换32“六大矩阵分解”A=QRA为m

n矩阵(先假设mn)Q为正交阵，R为上(右)三角阵正交阵的特点？元素值;正交性;转置;乘积;保2、F范数精简形式:Q为正交阵的前n列(列正交阵)若m<n定理:

对实矩阵𝑨,一定存在QR分解;若𝑨列满秩,且要求𝑹的对角元都>0,精简QR分解唯一矩阵的QR分解33AAAQQA=QRRR>>[Q,R]=qr(A);>>[Q,R]=qr(A,0);(用Cholesky分解唯一性证明)基于Gram-Schmidt正交化的算法O

计算过程保证了{qi}中任意两元素正交34fork

=1tonend基于Gram-Schmidt正交化的算法ClassicalG-S

procedurefork

=1tonforj

=1tok-1endendModifiedG-S

procedure计算量:

~mn2用剩下的向量做投影，而不是用原始的,提高准确度得到精简QR分解若A不列满秩，算法中断正交性丧失问题fork

=1tonforj

=1tok-1endend相邻的qk与qk+1之间正交性更好重正交化：35

基于Householder变换的算法

S

36

基于Householder变换的算法

可能同符号数相减!

修改后无抵消37

(2m次乘法,4mflops)

验证：基于Householder变换的算法

38完整QR分解算法39

(稳定的算法！)flops

完整QR分解算法40

基于Givens变换的算法

将二阶Givens阵

嵌入n阶单位阵:

41

基于Givens变换的算法例:

用Givens旋转变换进行消元每个Givens旋转阵用参数c,s刻画每次旋转变换仅影响向量两个元素则:

解:先针对第1,3分量构造二阶旋转矩阵,

接着对第1,4分量旋转

,仅影响矩阵的两行

42基于Givens变换的算法做若干次Givens变换实现对稠密矩阵用Givens旋转来正交三角化,比基于Householder变换的方法计算量约多50%与Householder变换需存u向量不同，Givens旋转需存储c,s参数值，存储用量约是前者的两倍适合处理稀疏矩阵A(填入现象，通过列交换减少填入)A实现第1列消元实现第2列消元

(不影响第1列)...得到R实现AE=QR求解稀疏线性最小二乘回归43线性最小二乘问题解法概述法方程方法利用QR分解解线性最小二乘列不满秩的情况44

求解线性最小二乘Amn

奇异45Ax

b

被拟合函数值基函数采样点

求解线性最小二乘正交变换保2范数设

n行

取等号条件:46n行

Ax

b

不满秩情况

=Ax

b47QR分解近似0

矩阵R1近奇异怎么求某一个解？解不唯一!(Householder变换)(这种满秩可能由数值误差造成,将其看成秩为1更合理)Rankrevealing问题做正交变换时适当地交换列在未消去的子阵中找2-范数最大的列,交换到当前列若第k步消去后,第k+1列的2-范数为0,则后续列均为0,且

k=rank(A