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