Factoring matrices into elementary matrices

Question

Is there a package in MATLAB, Maple or Mathematica which does this?

Answer 1

There are a number of factorization functions present in MATLAB, like LU factorization, orthogonal-triangular decomposition, and block LDL' factorization to name just a few.

Answer 2

The Matrix Decompositions page in the Mathematica documentation lists all the built-in matrix decomposition functions, such as SingularValueDecomposition, LUDecomposition, CholeskyDecomposition, SchurDecomposition, etc.

HTH!

Answer 3

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