将矩阵分解为初等矩阵

发布于 2024-09-12 12:04:42 字数 47 浏览 4 评论 0原文

MATLAB、Maple 或 Mathematica 中是否有包可以执行此操作？

原文

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

[旋木] 2024-09-19 12:04:43

我认为“基本”矩阵仅指那些执行行交换、行乘法和行加法基本运算的矩阵。

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

现在我们假设我们不需要进行任何行交换。因此，给定矩阵 M，我们可以得到下三角 L 和上三角 M。L 中位于主对角线下方的项是用来构造初等行加法矩阵的值。

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

表T2中的第二组基本行加法矩阵是由于得到u1^的GE步骤而产生的%T（M 的 PLU 分解中的 U 的转置）到行梯形形式。它们是使用 u1^%T 的 PLU 分解的 L 中的条目构建的。

只剩下 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);

I suppose that by "elementary" matrices you mean only those which do the elementary operations of row swapping, row multiplication, and row addition.

You might be interested to know that is part of the result of the PLU decomposition (factorization). The U coming from the PLU decomposition is the result of Gaussian Elimination, and that PLU decomposition is merely GE in disguise. And the P and L of the PLU decomposition encode the elemetary operations taken to accomplish GE. And all of Maple, Matlab, and Mathematica have a good PLU decomposition routine. So you can get the elementary factors.

Let's assume for now that we won't need to do any row swaps. So given matrix M we can get lower triangular L and upper triangular M. The entries of L that lie below the main diagonal are the values that with which to construct the elementary row addition matrices.

At the end is Maple code to show how it can be done. There are three sets of elementary matrices getting produced there. The first set, in table T1, are due to the GE steps of getting M to row echelon form and come from using l1 the L of the PLU decomposition of M. That's the lower triangle done. Next we'll transpose u1, the U of the PLU decomposition of M, so as to deal with the upper triangle of M.

The second set of elementary row addition matrices, in table T2, are due to the GE steps of getting u1^%T (the transposition of the U from PLU decomposition of M) to row echelon form. They are constructed using entries in l2 the L of the PLU decomposition of u1^%T.

That just leaves u2 the U of the PLU decomposition of u1^%T. It is a diagonal matrix (if there were no row swaps performed). So we construct elementary row scaling matrices for each row of u2.

Finally, we can put then all in the right order, and multiply them together. Note that the T2 matrices appear in reverse order, transposed, because they must multiply together to form u1^%T. Similarly, the T3 appear in between the T1 and T2 sets, since the T3 construct u2.

As a later edit, here it is as a Maple procedure. Now it generates the row swap matrices from the permutation results. And it doesn't return some unnecessary factors which happen to be just the identity.

Note that this is for exact Matrices, not floating-point Matrices (for which your mileage may vary, on account of how pivots are chosen by magnitude and how it does comparisons).

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);

回复收藏 0 原文

夜血缘 2024-09-19 12:04:43

Mathematica 文档中的矩阵分解页面列出了所有内置矩阵分解函数，例如 SingularValueDecomposition、LUDecomposition, CholeskyDecomposition，SchurDecomposition 等

HTH！

回复收藏 0 原文

恋你朝朝暮暮 2024-09-19 12:04:43

MATLAB 中存在许多因式分解/分解函数（请参阅此处的列表在“特征值和奇异值”和“矩阵分解”下），例如 LU 分解、正交三角分解和 < a href="https://www.mathworks.com/help/matlab/ref/ldl.html" rel="nofollow noreferrer">阻止 LDL 分解仅举几例。