MathGroup Archive 2001

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

Search the Archive

Re: Echelon form of a matrix

  • To: mathgroup at smc.vnet.net
  • Subject: [mg28492] Re: Echelon form of a matrix
  • From: "Allan Hayes" <hay at haystack.demon.co.uk>
  • Date: Sun, 22 Apr 2001 01:30:24 -0400 (EDT)
  • References: <9bosuo$mjl@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

Rasmus,
Here are two different approaches, not properly packaged or provided with
options at the moment.
I would be interested to here how they perform.

In the following two methods I have tried to
-- maximize the use of Mathematicas matrix functions ( I have also used
Position directly to find pivots);
--  minimize the simplifcation used during reduction.
The final result can be simplified if desired.

If expression build-up is a problem during reduction then some more
simplification can be introduced into the reduction process.

*** RowReduce: no scaling; division used in reduction ***

pivotPosition[w_]:=
  Reverse/@
    Position[Transpose[w], x_/;!Simplify[x===0],{2},1, Heads\[Rule]False]

liftRow[w_,i_] := Block[{wi=w},wi[[{1,i}]]=w[[{i,1}]];wi]

clearCol[w_] :=
  Rest[w] - Simplify[Outer[Times,Rest[w[[All,1]]], w[[1]]]/w[[1,1]]]

pasteBottomRight[m_,w_]:=
  Block[{mi=m},
    mi[[Sequence@@(Range[-#,-1]&/@Dimensions[w])]]=w;mi
    ]

step[{w_,m_}]:=
  Block[{wi=w,pos, mi=m,i,j,s=t},
    If[Length[wi]>1,
      pos =pivotPosition[w];
        If [pos=!={},
        {i,j}= pos[[1]];
        If[j>1,wi= Drop[wi,{},j-1]];
        If[i>1,wi= liftRow[wi, i]];
        wi = clearCol[wi];
        mi = pasteBottomRight[mi,wi];
        wi= Drop[wi,{},{1}],
        Throw[{wi,mi}]
        ],
      Throw[{wi,mi}]
      ];
    {wi,mi}
    ]

RowEchelon[m_]:=Last[Catch[Nest[step,{m,m}, Length[m]]]]

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

Test

RowEchelon[m]//Timing

Simplify[%][[2]]//Timing

*** RowEchelon2: no scaling; notdivision in reduction; divide at end ***

Three modifications:

pivotPosition2[w_]:=
  Reverse/@Position[Transpose[w], x_/;!Expand[x]===0,{2},1,
Heads\[Rule]False]

clearCol2[w_] := w[[1,1]]Rest[w] -Outer[Times,Rest[w[[All,1]]], w[[1]]]

step2[{w_,m_, d_}]:=
  Block[{wi=w,pos, mi=m,i,j,s=t,di=d},
    If[Length[wi]>1,
      pos =pivotPosition2[w];
        If [pos=!={},
        {i,j}= pos[[1]];
        If[j>1,wi= Drop[wi,{},j-1]];
        If[i>1,wi= liftRow[wi, i]];
        k= 1-Length[wi];
        di[[Range[k,-1]]]=di[[Range[k,-1]]]/w[[1,1]];
        wi = clearCol2[wi];
        mi = pasteBottomRight[mi,wi];
        wi= Drop[wi,{},{1}],
        Throw[{wi,mi,di}]
        ],
      Throw[{wi,mi,di}]
      ];
    {wi,mi, di}
    ]

RowEchelon2[m_]:=(#[[2]] #[[3]])&[
    Catch[Nest[step2,{m,m, Table[1,{Length[m]}]}, Length[m]]]]

Test

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

RowEchelon2[m]//Timing

Simplify[%]//Timing

--
Allan
---------------------
Allan Hayes
Mathematica Training and Consulting
Leicester UK
www.haystack.demon.co.uk
hay at haystack.demon.co.uk
Voice: +44 (0)116 271 4198
Fax: +44 (0)870 164 0565

"Rasmus Debitsch" <debitsch at zeiss.de> wrote in message
news:9bosuo$mjl at smc.vnet.net...
> 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
>
>
>
>




  • Prev by Date: Re: PC - Mac Compatibility
  • Next by Date: Re: PC - Mac Compatibility
  • Previous by thread: Echelon form of a matrix
  • Next by thread: cylindrical vector plot