矩阵程序大作业—说明文档.doc

LU分解、QR分解Householder reduction、Givens reduction1. 输入输出参数的含义Factorization_and_Reduction.m中输入输出参数的含义：输入参数：A：待分解的矩阵 mode：工作模式的选择 mode=0：函数进行LU分解 mode=1：函数进行QR分解 mode=2：函数进行Household reduction mode=3：函数进行Givens reduction输出参数：flag_right：分解或约减成功标志，为0时表示失败，为1时表示成功 B1,B2,B3，在不同模式下的含义如下 mode=0：B1=L：下三角阵； B2=U：上三角阵； B3=P：permutation matrix； 满足PA=LU。 mode=1：B1=Q：R(A)的单位正交基为列组成的矩阵； B2=R：上三角阵； B3=0：没有实际意义，默认值为0； 满足A=QR。 mode=2：B1=P：由反射得到的单位正交阵 (unitary matrix or orthogonal matrix)； B2=T：上梯形阵； B3=P_T：P的转置； 满足PA=T和A=P_T*T。 mode=3：B1=P：由旋转得到的单位正交阵 (unitary matrix or orthogonal matrix)； B2=T：上梯形阵； B3=P_T：P的转置； 满足PA=T和A=P_T*T。2. 实验步骤和原理1）LU分解（1）输入待分解的矩阵：（2）判断输入矩阵A是否为方阵，如果不是，不进行LU分解。（3）根据矩阵A的大小，初始化矩阵L，U，P和标签label（用于记录行交换后行的位置，将其初始化为）。（4）j从1到n循环执行以下步骤（j表示第j个主元）：在第j列以下位置中寻找第一个非零元素作为第j个主元，即从中寻找第一个非零元素作为主元。如果找到的非零主元不在第j行，而是在第i行，分别交换A和L中第i行和第j行，交换label中第i个和第j个元素。将A的第j行赋给U的第j行，并将L中位置元素赋1。将A中主元以下的元素变为0，并将消元时主元行需要乘的系数赋给L的相应位置元素，如：将主元以下的变为0时，L中位置的元素。（5） 对矩阵P进行赋值，即将P的第i行的第label(i)个元素赋为1，其他元素为0。（6） 判断矩阵A是否为奇异阵：查看矩阵U的主元位置是否存在0元素，如果存在则U为奇异阵且A也为奇异阵，此时A不能进行LU分解。（7） 对分解结果进行验证，计算是否成立，若成立，则分解正确。2） QR分解（1）输入待分解的矩阵：（2）k从1到n循环执行以下步骤：如果k=1，（表示矩阵A的第一列）； 否则，（表示矩阵A的第k列，表示矩阵Q的第i列，）。矩阵Q的第k列：给矩阵R的第k行赋值：（3） 循环结束即可得到Q矩阵和R矩阵。下面验证以上QR分解的结果是否正确：验证Q矩阵的各列是否相互正交且各列的模为1，即验证是否成立。验证是否成立。若以上两个条件均满足，则QR分解正确。3） Householder reduction（1）输入待处理的矩阵：（2）初始化矩阵P为单位矩阵I。（3）k从1到n循环执行以下步骤：令矩阵Ak是P*A的第k行k列下方右方构成的矩阵，如k=2，Ak矩阵如下图红框所示：Ak令构造m*m的矩阵更新P：用P*A的第k+1行和k+1列下方右方的元素构成的矩阵更新Ak（4） 以上循环结束以后，即可得到矩阵P，矩阵T=P*A。（5） 下面验证以上Householder reduction的结果是否正确： 验证P是否是单位正交矩阵(unitary matrix or orthogonal matrix)，即验证是否成立。验证是否成立。 若以上两个条件均满足，则Householder reduction分解正确。4）Givens reduction（1）输入待处理的矩阵：（2）初始化矩阵P为单位矩阵I。（3）k从1到n循环执行以下步骤：利用旋转矩阵，使得原矩阵k列中k行元素以下的元素为0 则j从k+1到m循环执行以下步骤： 令矩阵Ak等于P*A ， 令，然后用c去替代中第k行k列的元素，用s去替代中第k行j列的元素，用-s去替代中第j行k列的元素，用c去替代中第j行j列的元素。 更新P： 更新Ak：（4）以上循环结束以后，即可得到矩阵P，矩阵T=P*A。（5）下面验证以上Givens reduction的结果是否正确：验证P是否是单位正交阵(unitary matrix or orthogonal matrix)，即验证是否成立。验证是否成立。 若以上两个条件均满足，则Givens reduction分解正确。3. 实验结果截图1) 当输入矩阵A为方阵例：（1） LU分解显示结果中变量含义：L：下三角阵；U：上三角阵；P：permutation matrix。（2） QR分解显示结果中变量含义：Q：R(A)的单位正交基为列组成的矩阵；R：上三角阵。（3） Householder reduction显示结果中变量含义：P_H：由反射得到的单位正交阵(unitary matrix or orthogonal matrix)；T_H：上梯形阵；P_T_H：P_H的转置。（后缀_H用于与Givens的结果区分）（4） Givens reduction显示结果中变量含义：P_G：由旋转得到的单位正交阵(unitary matrix or orthogonal matrix)；T_G：上梯形阵；P_T_G：P_G的转置。（后缀_G用于与Householder的结果区分）2) 当输入矩阵A不是方阵例：（1）LU分解（2）QR分解显示结果中变量含义：Q：R(A)的单位正交基为列组成的矩阵；R：上三角阵。（3）Householder reduction显示结果中变量含义：P_H：由反射得到的单位正交阵(unitary matrix or orthogonal matrix)；T_H：上梯形阵；P_T_H：P_H的转置。（后缀_H用于与Givens的结果区分）（4）Givens reduction显示结果中变量含义：P_G：由旋转得到的单位正交阵(unitary matrix or orthogonal matrix)；T_G：上梯形阵；P_T_G：P_G的转置。（后缀_G用于与Householder的结果区分）4. 程序代码1） RUN_FR.mclc;clear all;%赵晓梅 201428014628066 2014.11.19%提供几个待分解的矩阵% A=4,-3,4;2,-14,-3;-2,14,0;1,-7,15;%待分解矩阵1非方阵A=0,-20,-14;3,27,-4;4,11,-2;%待分解矩阵2方阵fprintf(待处理的矩阵为n);A%在Command Window中显示待处理矩阵fprintf(*n); L,U,P,flag_right_0 = Factorization_and_Reduction( A,0 );%LU分解if flag_right_0=1%如果分解正确，显示结果 fprintf(显示LU分解结果：n); L %下三角阵； U %上三角阵； P %permutation matrix；endfprintf(*n); Q,R,B3,flag_right_1 = Factorization_and_Reduction( A,1 );%QR分解if flag_right_1=1%如果分解正确，显示结果 fprintf(显示QR分解结果：n); Q %R(A)的单位正交基为列组成的矩阵； R %上三角阵；endfprintf(*n); P_H,T_H,P_T_H,flag_right_2 = Factorization_and_Reduction( A,2 );%Householder reductionif flag_right_2=1%如果Householder reduction正确，显示结果 fprintf(显示Householder reduction结果（后缀_H用于与Givens区分）：n); P_H %由反射得到的单位正交阵； T_H %上梯形阵； P_T_H %P_H的转置；endfprintf(*n); P_G,T_G,P_T_G,flag_right_3 = Factorization_and_Reduction( A,3 );%Givens reductionif flag_right_3=1%如果Givens reduction正确，显示结果 fprintf(显示Givens reduction结果（后缀_G用于与Householder区分）：n); P_G %由反射得到的单位正交阵； T_G %上梯形阵； P_T_G %P_G的转置；end2） Factorization_and_Reduction.m%-%函数名称：Factorization_and_Reduction%函数功能：对矩阵进行LU分解，QR分解，Household reduction，Givens reduction%输入参数：A：待分解的矩阵% mode：工作模式的选择% mode=0：函数进行LU分解% mode=1：函数进行QR分解% mode=2：函数进行Householder reduction% mode=3：函数进行Givens reduction%输出参数：flag_right：分解或约减成功标志，为0时表示失败，为1时表示成功% B1,B2,B3，在不同模式下的含义如下% mode=0：B1=L：下三角阵；% B2=U：上三角阵；% B3=P：permutation matrix；% 满足PA=LU。% mode=1：B1=Q：R(A)的单位正交基为列组成的矩阵；% B2=R：上三角阵；% B3=0：没有实际意义，默认值为0；% 满足A=QR。% mode=2：B1=P：由反射得到的单位正交阵(unitary matrix or orthogonal matrix)；% B2=T：上梯形阵；% B3=P_T：P的转置；% 满足PA=T和A=P_T*T。% mode=3：B1=P：由旋转得到的单位正交阵(unitary matrix or orthogonal matrix)；% B2=T：上梯形阵；% B3=P_T：P的转置；% 满足PA=T和A=P_T*T。%基本信息：赵晓梅 201428014628066 2014.11.19%-function B1,B2,B3,flag_right = Factorization_and_Reduction( A,mode )B1=0;B2=0;B3=0;flag_right=0;%当不满足分解或约减条件时，为保证有输出，将参数初始化为0m,n=size(A);%获取矩阵A的行数和列数%-if mode=0%对矩阵进行LU分解 fprintf(LU分解n); A0=A;%分解过程中会对A中的元素进行改变，为方便验证P*A=L*U，用A0保存A的原值 if m=n fprintf(输入矩阵不是方阵，不能进行LU分解n); return; end L=zeros(n,n);%L,U,P初始化 U=zeros(n,n); P=zeros(n,n); label=(1:1:n);%标签，标记进行的行变换 t=zeros(1,n);%用于行交换时的中间变量 for j=1:n%j表示第j个主元 for i=j:n if A(i,j)=0%在第j列，寻找不为0的主元 break; end end if i=j%如果不为0的主元不在位置（j,j）,则对主元所在的第i行和第j行进行交换 t=A(j,:); A(j,:)=A(i,:); A(i,:)=t; l=label(j);%label相应位置进行交换 label(j)=label(i); label(i)=l; t=L(j,:);%L的相应位置交换 L(j,:)=L(i,:); L(i,:)=t; end U(j,:)=A(j,:);%交换以后的A的第j行赋给U的第j行 L(j,j)=1;%L的（j,j）位置赋1 for i=(j+1):n%将主元以下变成0 L(i,j)=A(i,j)/A(j,j); A(i,:)=A(i,:)-A(j,:)*L(i,j); end end for i=1:n P(i,label(i)=1; end flag_right=1; for i=1:n if U(i,i)=0;%若U的主元位置存在0，则U为奇异阵，A也为奇异阵 flag_right=0; end end if flag_right=0; fprintf(A是奇异矩阵不能进行LU分解n); else %验证LU分解是否正确 fprintf(验证结果是否正确n); if P*A0=L*U fprintf(P*A=L*U，LU分解正确进行n); else fprintf(P*A不等于L*U，LU分解错误n); end end B1=L; B2=U; B3=P;end%-if mode=1%QR分解 fprintf(QR分解n);% m,n=size(A); Q=zeros(m,n);%初始化Q、R的大小 R=zeros(n,n); for k=1:n u=A(:,k); if k1 for i=1:k-1%矩阵A中的第k列分别减去它在Q的第1k-1列上的投影得到u u=u-Q(:,i)*A(:,k)*Q(:,i); end end v=norm(u); Q(:,k)=u/v;%u单位化后赋给Q的第k列 for j=k:n%矩阵R第k行的元素分别为A的第j列构成的向量在Q的第k列构成的向量上的投影的大小 R(k,j)=Q(:,k)*A(:,j); end end %验证分解结果是否正确 fprintf(验证结果是否正确n); I=zeros(n,n); for i=1:n I(i,i)=1; end eps=10(-10); if norm(Q*Q-I)eps%验证Q的各列是否正交且各列的模为1，即Q的转置乘以Q是否是单位阵 flag_right_1=1; fprintf(1. Q的转置乘以Q是单位阵,Q各列相互正交且各列的模为1n); else flag_right_1=0; fprintf(1. Q的转置乘以Q不是单位阵,Q不满足各列正交且各列的模为1n); end if norm(A-Q*R)eps flag_right_2=1; fprintf(2. 分解结果满足A=Q*Rn); else flag_right_2=0; fprintf(2. 分解结果不满足A=Q*Rn); end if flag_right_1=0|flag_right_2=0 fprintf(所以，QR分解错误n); flag_right=0; else fprintf(所以，QR分解正确n); flag_right=1; B1=Q; B2=R; endend%-if mode=2%Householder reduction fprintf(Householder reductionn); I0=zeros(m,m); for i=1:m I0(i,i)=1; end P=I0;%将P初始化为单位矩阵 Ak=A; for k=1:n mm,nn=size(Ak);%Ak中保存的是第k行k列元素下方和右方的元素组成的矩阵 if mm=1 break; end e=zeros(mm,1); e(1)=1; u=Ak(:,1)-norm(Ak(:,1)*e; I=zeros(m-k+1,m-k+1); for i=1:m-k+1 I(i,i)=1; end R_hat=I-2*u*u/(u*u);%获得能使Ak矩阵第一列只有第一个元素不为0的矩阵 R=I0;%把R_hat矩阵放到R矩阵的右下方，此时的R矩阵能使P*A的第k列中k行以下的元素全为0 m_R,n_R=size(R_hat); for i=(m-m_R+1):m for j=(m-n_R+1):m R(i,j)=R_hat(i-m+m_R,j-m+n_R); end end P=R*P;%将R与P相乘 Ak=P*A;%P与A相乘，并取之k行k列以下的元素构成新的矩阵更新Ak Ak_hat=zeros(m-k,n-k); for i=k+1:m for j=k+1:n Ak_hat(i-k,j-k)=Ak(i,j); end end Ak=Ak_hat; end T=P*A; %验证Householder reduction结果是否正确 fprintf(验证结果是否正确n); I=zeros(m,m); for i=1:m I(i,i)=1; end eps=10(-10); if norm(P*P-I)eps & norm(P*P-I)eps%依据P的逆是否等于P的转置，即P的转置乘以P是否是单位阵，验证P是否是单位正交矩阵 flag_right_1=1; fprintf(1. P_T*P=I且P*P_T=I,P的确是单位正交矩阵(unitary matrix or orthogonal matrix) （P_T表示P的转置）n); else flag_right_1=0; fprintf(1. 不满足P_T*P=I且P*P_T=I，P不是是单位正交矩阵(unitary matrix or orthogonal matrix)（P_T表示P的转置）n); end if norm(A-P*T)eps flag_right_2=1; fprintf(2.分解结果满足A=P_T*T (P_T表示P的转置)n); else flag_right_2=0; fprintf(2.分解结果不满足A=P_T*T (P_T表示P的转置)n); end if flag_right_1=0|flag_right_2=0 fprintf(所以，Householder reduction错误n); flag_right=0; else fprintf(所以，Householder reduction正确n); flag_right=1; B1=P; B2=T; B3=P; endend%-if mode=3 %Givens reduction fprintf(Givens reductionn); T=zeros(m,n); I=zeros(m,m); for i=1:m I(i,i)=1; end P=I; Ak=A; for k=1:n for j=k+1:m%依次构造旋转矩阵，使得j行k列的元素等于0 P_hat=