資源簡介
對于矩陣A,根據不同的參數number實現不同方法的矩陣分解:高斯消元法的LU分解,基于Gram-Schmidt、Householder、Givens 的QR分解

代碼片段和文件信息
function?[QR]=Reduction(Anumber)
if?number==1
?[QR]=LUcomposition(A);
elseif?number==2
??????[QR]?=?GramSchmidt(A);
elseif?number==3
??????[QR]?=?HouseHolder(A);
elseif?number==4
??????[QR]=Givens(A);
elseif?error(‘not?valid?number‘);
end
end
function?[LU]=LUcomposition(A)?%A=[2?2?2;4?7?7;6?18?22]??A=[10?20?30;20?45?80;30?80?171]
B=A;
[mn]=size(A);
if?m~=n
????error(‘不是方陣!‘);
end
if?rank(A)~=n
???error(‘奇異矩陣無法進行LU分解!‘);?%檢測參數異常停止程序
end
L=eye(n);
for?i=1:n-1
????if?abs(A(ii)) ????????disp(‘zero?pivot?encountered--use?PA=LU?‘);
????????[LU]=partialLUDecomposition(Bn);
????????return;
????end
????for?j=i+1:n????????????
????????????L(ji)=A(ji)/A(ii);
????????????A(j:)=A(j:)-(A(ji)/A(ii))*A(i:);
????end
end
U=A;
disp(‘-----LU?Reduction-----?‘);
end
function?[LU]=partialLUDecomposition(Am)??%?A=[1?2?-3?4;4?8?12?-8;2?3?2?1;-3?-1?1?-4]
L=zeros(m);
P=eye(m);
for?i=1:m-1
????index=find(abs(A(i:mi))==max(abs(A(i:mi))));
????k=index+i-1;
????if?k~=i
?????L([i?k]:)=L([k?i]:);
?????A([i?k]:)=A([k?i]:);?
?????P([i?k]:)=P([k?i]:);?
????end
????for?j=i+1:m??
????????L(ji)=A(ji)/A(ii);
????????A(j:)=A(j:)-(A(ji)/A(ii))*A(i:);
????end
end
U=A;
L=L+eye(m);?
disp(‘--partial?LU?Reduction---‘);?
end
%?對m*n維的矩陣A具有線性無關列,進行基于Gram-Schmidt的QR分解:A=QR,其中Q為m*n維正交矩陣,R為n*n維上三角矩陣且對角線元素大于0
function?[QR]?=GramSchmidt(A)??%A=[0?-20?-14;3?27?-4;4?11?-2]
[m?n]?=?size(A);
if?rank(A)~=n
???error(‘該矩陣無法進行進行基于Gram-Schmidt的QR分解!‘);?%檢測參數異常停止程序
end
Q?=?zeros(?m?n?);
R?=?zeros(?n?n?);
R(11)=?norm(?A(?:?1)?);
Q(:?1)?=?A(?:?1)/R(?1?1?);
?for?k?=?2:?n
?????for?j?=?1:?k?-?1
?????????R(?j?k?)?=?Q(?:?j?)‘?*?A(?:?k);
?????end
?????q?=?zeros(?m?1?);????
?????for?j?=?1:?k?-?1
?????????q?=?Q(?:?j)?*?R(?j?k?)?+?q;
?????end???
?????R(?k?k?)?=?norm(?A(?:?k?)?-?q?);
?????Q(?:?k?)?=?(?A(?:?k)?-?q?)?/?R(?k?k?);
?end
disp(‘---Gram-Schmidt?Reduction---‘);?
end
%%houseHolder進行QR分解?
function?[QR]=HouseHolder(A)?
[m?n]?=?size(A);?%?m-?number?of?rows?n-?number?of?columns
if?m?>?n
????num?=?n;%the?number?of?housholder?matrix
else
????num?=?m?-?1;
end
a=A;
for?i=1:num
????if?i==1?
???????u=a(:1)-norm(a(:1))*eye(m1);?
???????R=eye(m)-2*(u*u‘)/(u‘*u);?
???????Q=R;?
???????R=R*A;
???????fR=R;
????else
????????a=R(i:mi:n);?
????????u=a(:1)-norm(a(:1))*eye(m-i+11);?
????????R1=eye(m-i+1)-2*(u*u‘)/(u‘*u);
????????R=[eye(i-1)?zeros(i-1m-i+1);zeros(m-i+1i-1)R1];?
????????Q=Q*R;?
????????R=R*fR;?
????????fR=R;
????end
end
disp(‘------Householder?Reduction------‘);?
end
function?[QR]=Givens(A)%A是
[mn]=size(A);?%m行數n列數?A=[1?19?-34;-2?-5?20;2?8?37]
R=A;
Q=eye(m);
for?i=1:n-1?
????for?j=i+1:m?
????????x=R(:i);
????????G=givensmatrix(xij);
????????%Q=Q*G‘;
????????Q=G*Q;
????????R=G*R;
????end
end
Q=Q‘;
disp(‘----Givens?Reduction----‘);
end
function?G=givensmatrix(xij)
xi=x(i);??????????
xj=x(j);
r=sqrt(xi^2+xj^2);
c=x
?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????文件????????3187??2018-12-26?00:16??矩陣分解\Reduction.m
?????文件???????43759??2018-12-26?00:36??矩陣分解\說明文檔.docx
評論
共有 條評論