Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2001
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2001

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

Search the Archive

Re: Echelon form of a matrix lCorrection]

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

I noticed some bugs in my pevious posting, and Rasmus Debitsch kindly sent
me some useful examples that brought up another one. Apologies for these.
Here is a corrected and improved versions of  RowEchelon (RowEchelon2 seems
too slow for larger matrices)

Also, in the original posting DixonGE used the options

    GESimplify->Together,
    GEZeroTest->(Simplify[#===0]&)

Here Simplify can do nothing since each instance of #===0 will evaluate to
True or false before Simplify acts.

We are really working with

    GESimplify->Together,
    GEZeroTest->(#===0]&)

If are not concerned with big final entries then we can use

    GESimplify->Identity
    GEZeroTest->Together[#]===0&)

which seems quicker.

I suspect that the method used in DixonGE will be quicker for big matrices,
particularly sparse ones.

*** RowEchelon ***

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

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

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

pasteBottomRight[m_, w_] := PadLeft[w, Dimensions[m], m]

step[{w_, m_}] :=
  Block[{wi, pos, mi, i, j},
    If[MatchQ[w, {_} | {{}, ___}], Throw[{w, m}] ];
    pos = pivotPosition[w];
    If [pos === {}, Throw[{w, m}]];
    {i, j} = pos[[1]];
    wi = Drop[w, {}, j - 1];
    If[i > 1, wi[[1]] = -wi[[1]]];  (*additional step to previous code*)
    wi = liftRow[wi, i];
    wi = clearCol[wi];
    mi = pasteBottomRight[m, PadLeft[wi, Dimensions[w]]];
    wi = Drop[wi, {1}, {1}];
    {wi, mi}
    ]

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



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

"Allan Hayes" <hay at haystack.demon.co.uk> wrote in message
news:9btr8n$1uh at smc.vnet.net...
> 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: RealTime3D without RealTime3D
  • Next by Date: nxn matrix into point coordinate
  • Previous by thread: Re: RealTime3D without RealTime3D
  • Next by thread: nxn matrix into point coordinate