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