matlab - 将矩阵因式分解为初等矩阵

Question

MATLAB、Maple 或 Mathematica 中是否有执行此操作的程序包？

Answer 1

我想您所说的“基本”矩阵仅指那些执行行交换、行乘法和行加法等基本操作的矩阵。

您可能有兴趣知道这是 PLU 分解(因式分解)结果的一部分。 PLU分解得到的U是高斯消去的结果，PLU分解只是变相的GE。 PLU 分解的 P 和 L 对完成 GE 所采取的基本操作进行编码。 Maple、Matlab 和 Mathematica 都有很好的 PLU 分解例程。所以你可以获得基本因素。

现在假设我们不需要进行任何行交换​​。所以给定矩阵M我们可以得到下三角L和上三角M。位于主对角线下方的L的元素是构造基本行加法矩阵的值。

最后是 Maple 代码，展示了它是如何完成的。那里产生了三组基本矩阵。第一个集合，在表 T1 中，是由于将 M 变为行阶梯形式的 GE 步骤，并且来自使用 l1 M 的 PLU 分解的 L。这是完成的下三角。接下来将M的PLU分解的U转置u1，以处理M的上三角。

表T2中的第二组基本行加法矩阵是由于将u1^%T(M的PLU分解得到的U转置)到行阶梯形式的GE步骤。它们是使用 l2 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);