matlab自定义矩阵分解,matlab – 将矩阵分解为基本矩阵

我认为“基本”矩阵只表示那些进行行交换,行乘法和行加法的基本操作的矩阵.

您可能有兴趣知道这是PLU分解(分解)结果的一部分.来自PLU分解的U是高斯消除的结果,并且PLU分解仅仅是伪装的GE.并且PLU分解的P和L编码用于完成GE的元素操作.所有Maple,Matlab和Mathematica都有一个很好的PLU分解程序.所以你可以得到基本因素.

我们现在假设我们不需要进行任何行交换??.因此,给定矩阵M,我们可以得到下三角形L和上三角形M.位于主对角线下方的L的条目是用于构造基本行加法矩阵的值.

最后是Maple代码,展示如何完成它.那里有三套基本矩阵.表T1中的第一组是由于GE的步骤是将M形成行梯形形式,并且使用l1来表示M的PLU分解的L.这是下三角形.接下来我们将转换u,即M的PLU分解的U,以便处理M的上三角形.

表T2中的第二组基本行加法矩阵是由于获得u1 ^％T(从M的PLU分解将U的转置)转换为行梯形形式的GE步骤.它们是使用l1中的条目构建的,即u1 ^％T的PLU分解的L.

这只留下u2 U1 ^％T的PLU分解的U.它是一个对角矩阵(如果没有执行行交换).因此,我们为u2的每一行构造基本行缩放矩阵.

最后,我们可以按顺序放置所有内容,然后将它们相乘.注意,T2矩阵以相反的顺序出现,转置,因为它们必须相乘以形成u1 ^％T.类似地,T3出现在T1和T2集之间,因为T3构造了u2.

作为后来的编辑,这里它是一个Maple过程.现在它从排列结果生成行交换矩阵.它并没有返回一些恰好只是身份的不必要的因素.

请注意,这是针对精确矩阵,而不是浮点矩阵(由于如何根据大小选择枢轴以及如何进行比较,您的里程可能会有所不同).

ElemDecomp:=proc(M::Matrix(square))

local p1,u1,i,j,T1,T2,T3,p2,m,n,lu1,lu2,P1,P2;

uses LinearAlgebra;

(m,n):=Dimensions(M);

p1,lu1:=LUDecomposition(M,output=[':-NAG']);

for i from 1 to m-1 do

for j from 1 to i do

if lu1[i+1,j]<>0 then

T1[i*j]:=IdentityMatrix(m,compact=false);

T1[i*j][i+1,j]:=lu1[i+1,j];

end if;

end do; end do;

for i from 1 to m do

if p1[i]<>i then

P1[i]:=IdentityMatrix(m,compact=false);

P1[i][p1[i],i],P1[i][i,p1[i]]:=1,1;

P1[i][p1[i],p1[i]],P1[i][i,i]:=0,0;

end if;

end do;

u1:=Matrix(lu1,shape=triangular[upper]);

p2,lu2:=LUDecomposition(u1^%T,output=[':-NAG']);

for i from 1 to m-1 do

for j from 1 to i do

if lu2[i+1,j]<>0 then

T2[i*j]:=IdentityMatrix(m,compact=false);

T2[i*j][i+1,j]:=lu2[i+1,j];

end if;

end do; end do;

for i from 1 to m do

if lu2[i,i]<>1 then

T3[i]:=IdentityMatrix(m,compact=false);

T3[i][i,i]:=lu2[i,i];

end if;

end do;

for i from 1 to m do

if p2[i]<>i then

P2[i]:=IdentityMatrix(m,compact=false);

P2[i][p2[i],i],P2[i][i,p2[i]]:=1,1;

P2[i][p2[i],p2[i]],P2[i][i,i]:=0,0;

end if;

end do;

`if`(type(P1,table),entries(P1,':-nolist'),NULL),

seq(seq(`if`(assigned(T1[i*j]),T1[i*j],NULL),j=1..i),i=1..m-1),

seq(`if`(assigned(T3[i]),T3[i],NULL),i=1..min(m,n)),

seq(seq(`if`(assigned(T2[i*j]),T2[i*j]^%T,NULL),j=i..1,-1),i=m-1..1,-1),

`if`(type(P2,table),entries(P2,':-nolist'),NULL);

end proc:

A:=LinearAlgebra:-RandomMatrix(3,generator=1..4);

ElemDecomp(A);

LinearAlgebra:-Norm( `.`(%) - A);