Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Echelon form of a matrix

• To: mathgroup at smc.vnet.net
• Subject: [mg28476] Echelon form of a matrix
• From: "Rasmus Debitsch" <debitsch at zeiss.de>
• Date: Fri, 20 Apr 2001 04:24:24 -0400 (EDT)
• Sender: owner-wri-mathgroup at wolfram.com

Hello Mathgroup.

To solve a problem in robotics, I'm want to use resultants of polynomial
systems. In "Mathematica in Education and Research", Vol.6 No.3 1997 one can
found an instructive article "Elimination with the Dixon Resultant" by Nakos
and Williams. Working with their implementation of the algorithms, it turns
out, that the computation of the echolon form of a matrix by using
elementary row operations without scaling is the main performance gap of the
algorithm. The elemination step for symbolic matrices generates very long
expressions and is very time comsuming.
I have the strong feeling, that somewhere is a build in Mathematica
function, which does the job much faster. But I haven't found this function
yet. One candidate is RowReduce -- but I haven't found a way to use it for
this special elemination process.

I have attached the implementation of the elimination procedure an an
example.

<< LinearAlgebra`MatrixManipulation`

Clear[DixonGE,DixonPivot];

DixonPivot[
    matrix_?MatrixQ,
    options___?OptionQ
    ]:=Module[{
      i=1,j=1,
      n,m,zerotest
      },
    zerotest=GEZeroTest/.Flatten[{options}]/.Options[Dixon];

    {n,m}=Dimensions[matrix];
    While[
      j<m&&VectorQ[matrix[[All,j]], zerotest[#]&],
      j++];
    While[
      i<n&&zerotest[matrix[[i,j]]],
      i++];
    {i,j}
    ];

DixonGE[
      matrix_?MatrixQ,
      options___?OptionQ
      ]:=Module[{
        n,m,i,j,matrixGE,matS,q,qr,
        verbose,print,simplify,zerotest
        },
      verbose=Verbose/.Flatten[{options}]/.Options[Dixon];
      print=VerbosePrint/.Flatten[{options}]/.Options[Dixon];
      simplify=GESimplify/.Flatten[{options}]/.Options[Dixon];
      zerotest=GEZeroTest/.Flatten[{options}]/.Options[Dixon];

      {n,m}=Dimensions[matrix];
      If[n<2,
        Return[matrix]
        ];

      If[verbose,print["DixonGE ",n," ",m]];

      Fold[
        Function[{mat,row},
          {i,j}=DixonPivot[TakeMatrix[mat,{row,row},{n,m}]]+{row-1,row-1};
          If[!VectorQ[{i,j},IntegerQ],
            mat,
            If[i!=row,

              matS=mat[[Range[n]/.{row->i,i->row}]];
              matS[[i]]*=-1,
              matS=mat;
              ];

            If[verbose,
              print["DixonGE:"," row ",row," pivot ",{i,j}];
              ];

            q=matS[[row,j]];
            If[Not[zerotest[q]],
              qr=matS[[row]]/q;
              MapAt[simplify[#-#[[j]] qr]&,matS,
Transpose[{Range[row+1,n]}]], matS ]
            ]
          ],
        matrix,
        Range[1,n-1]
        ]
      ];

Clear[m, a, b, c];

m={{-270*a*b, 36*b^2, 21*a^2*c}, {36*b^2, 108*a*b, 0}, {0, 9*b^2*c,
27*a*b*c},
 {9*b^2*c, 27*a*b*c, 0}};

DixonGE[m,
    Verbose->True,
    VerbosePrint->Print,
    GESimplify->Together,
    GEZeroTest->(Simplify[#===0]&)
]

The typical matrix size for my special problems ist 50x60. The resulting
expressions have leaf counts up to 1000000.

Regards
Rasmus

--

        Rasmus Debitsch
        Carl Zeiss Lithos GmbH
        Carl Zeiss Strasse
        D-73447 Oberkochen

        eMail : debitsch at zeiss.de