Re: Complexity explosion in linear solve
- To: mathgroup at smc.vnet.net
- Subject: [mg80363] Re: Complexity explosion in linear solve
- From: carlos at colorado.edu
- Date: Sun, 19 Aug 2007 01:30:47 -0400 (EDT)
- References: <f9s3f8$a8r$1@smc.vnet.net><fa3cq4$5g$1@smc.vnet.net>
On Aug 16, 11:47 pm, d... at wolfram.com wrote:
> > On Aug 15, 2:28 am, Jens-Peer Kuska <ku... at informatik.uni-leipzig.de>
> > wrote:
> >> Hi,
>
> >> what do you expect ? matrix inversion is an N^3 process
> >> when N is the dimension of the matrix. This mean for numerical
> >> data, that you have to do N^3 operations to get the inverse.
> >> For symbolic data a single operation gives not a single number
> >> it simply add another leaf to your expression because 1+2 is 3
> >> (leaf count 1) but a+b is a+b (leaf count 2) so every of your N^3
> >> operations add a new leaf to your expression and you end up (in the
> >> worst case) with N^3 more leafs in your expression ..
>
> >> Regards
> >> Jens
>
> >> car... at colorado.edu wrote:
> >> > LinearSolve (and cousins Inverse, LUDecomposition etc) works
> >> > efficiently with
> >> > numerical float data. But it is prone to complexity explosion in
> >> > solving symbolic
> >> > systems of moderate size. Here is an example from an actual
> >> > application.
>
> >> > Task: solve A x = b for x where A is 16 x 16 and b is 16 x 1:
>
> >> > A={{7/12,1/24,-1/2,0,-1/4,-5/24,-1/2,0,
> >> > 0,-1/24,-1/4,0,-1/3,5/24,-1/4,0},{1/24,91/144,-2/3,0,
> >> > 5/24,-3/16,-1/3,0,-1/24,0,-1/3,0,-5/24,-4/9,-2/3,
> >> > 0},{-1/2,-2/3,(-9+(1378*C11*hh)/225)/2,(-4*C13*hh)/15,
> >> > 1/2,-1/3,(-3-2*C11*hh)/2,4*C13*hh,1/4,
> >> > 1/3,(-3-(32*C11*hh)/25)/2,0,-1/4,
> >> > 2/3,(-3-(128*C11*hh)/45)/2,(-64*C13*hh)/15},{0,
> >> > 0,(-4*C13*hh)/15,(224*C22*hh)/5,0,0,-4*C13*hh,16*C23*hh,0,0,0,
> >> > 16*C23*hh,0,0,(64*C13*hh)/15,(64*C23*hh)/5},{-1/4,5/24,1/2,0,
> >> > 7/12,-1/24,1/2,0,-1/3,-5/24,1/4,0,0,1/24,1/4,
> >> > 0},{-5/24,-3/16,-1/3,0,-1/24,91/144,-2/3,0,5/24,-4/9,-2/3,
> >> > 0,1/24,0,-1/3,0},{-1/2,-1/3,(-3-2*C11*hh)/2,-4*C13*hh,
> >> > 1/2,-2/3,(-9+(314*C11*hh)/45)/2,(4*C13*hh)/15,1/4,
> >> > 2/3,(-3-(128*C11*hh)/45)/2,(64*C13*hh)/15,-1/4,
> >> > 1/3,(-3-(32*C11*hh)/15)/2,0},{0,0,4*C13*hh,16*C23*hh,0,
> >> > 0,(4*C13*hh)/15,(832*C22*hh)/15,0,
> >> > 0,(-64*C13*hh)/15,(64*C23*hh)/5,0,0,0,(80*C23*hh)/3},{0,-1/24,
> >> > 1/4,0,-1/3,5/24,1/4,0,7/12,1/24,1/2,0,-1/4,-5/24,1/2,
> >> > 0},{-1/24,0,1/3,0,-5/24,-4/9,2/3,0,1/24,91/144,2/3,0,
> >> > 5/24,-3/16,1/3,0},{-1/4,-1/3,(-3-(32*C11*hh)/25)/2,0,
> >> > 1/4,-2/3,(-3-(128*C11*hh)/45)/2,(-64*C13*hh)/15,1/2,
> >> > 2/3,(-9+(1378*C11*hh)/225)/2,(-4*C13*hh)/15,-1/2,
> >> > 1/3,(-3-2*C11*hh)/2,4*C13*hh},{0,0,0,16*C23*hh,0,
> >> > 0,(64*C13*hh)/15,(64*C23*hh)/5,0,0,(-4*C13*hh)/15,(224*C22*hh)/5,
> >> > 0,0,-4*C13*hh,16*C23*hh},{-1/3,-5/24,-1/4,0,0,1/24,-1/4,
> >> > 0,-1/4,5/24,-1/2,0,7/12,-1/24,-1/2,0},{5/24,-4/9,2/3,0,
> >> > 1/24,0,1/3,0,-5/24,-3/16,1/3,0,-1/24,91/144,2/3,
> >> > 0},{-1/4,-2/3,(-3-(128*C11*hh)/45)/2,(64*C13*hh)/15,
> >> > 1/4,-1/3,(-3-(32*C11*hh)/15)/2,0,1/2,
> >> > 1/3,(-3-2*C11*hh)/2,-4*C13*hh,-1/2,
> >> > 2/3,(-9+(314*C11*hh)/45)/2,(4*C13*hh)/15},{0,
> >> > 0,(-64*C13*hh)/15,(64*C23*hh)/5,0,0,0,(80*C23*hh)/3,0,0,
> >> > 4*C13*hh,16*C23*hh,0,0,(4*C13*hh)/15,(832*C22*hh)/15}}
>
> >> > b={0,0,0,0,q*b/2,0,0,0,q*b/2,0,0,0,0,0,0,0}
>
> >> > in which all symbolic variables are atoms. Matrix A has rank 13.
> >> > Three BCs: x1=x2=x13=0 are imposed by removing equations 1, 2 and
> >> > 13,
> >> > which gives the reduced 13-system Ahat xhat = bhat. Matrix Ahat is
> >> > symmetric but indefinite, so Cholesky is not an option. Instead
> >> > LinearSolve is used
>
> >> > xhat=LinearSolve[Ahat,bhat]
>
> >> > Solving takes about 3 mins on a dual G5 running 5.2. LeafCounts of
> >> > the xhat entries reach hundreds of millions:
>
> >> > {325197436,235675446,292306655,256982512,146153454,73076907,35324210,
> >> > 18877804,9441784,4745440,2429139,1354800,903023}
>
> >> > Obviously Simplify would take a long, long time so I didnt attempt it.
> >> > Another solution method, however, gave this solution in about 10 sec:
>
> >> > xhat={q/3, 0, 4*q, 0, q/3, 0, 4*q, 0, q/3, 0, 0, q/3, 0}
>
> >> > My question is: is there a way to tell LinearSolve to Simplify as it
> >> > goes
> >> > along? That would preclude or at least alleviate the leafcount
> >> > explosion.
>
> > Hi Jens,
>
> > I know *that* - I teach numerical analysis, among other stuff.
> > Also I teach the kids to avoid Cramer's rule for systems
> > of order >3. Then can you explain me why
>
> > MyLinearSolve[A_,b_]:=Module[{i,j,Ab,detA,detAb,n=Length[A],x},
> > detA=Simplify[Det[A]]; x=Table[0,{n}];
> > For [i=1,i<=n,i++, Ab=Transpose[A];
> > Ab[[i]]=b; detAb=Simplify[Det[Ab]];
> > x[[i]]=Simplify[detAb/detA] ];
> > ClearAll[Ab];
> > Return[x]];
>
> > beats LinearSolve for a symbolic system of order 13?
> > Here are my times on a dual G5 Mac:
>
> > LinearSolve 235 sec unsimplified x
> > MyLinraSolve 6 sec simplified x
>
> > For the simplification of the LinearSolve solution x
> > I estimate 10^8 years.
>
> Hard to say much without seeing the specific input you use. I'll venture a
> guess below.
>
> As for Cramer's rule...for symbolic systems it is often a "good" method
> but requires huge memory to do efficiently (one must memo-ize smaller
> cofactors in order to do the computations recursively). For sufficiently
> large dimension-- 12 and up if I recall correctly-- it is not used as a
> default method in LinearSolve et al. In any case it won't work for
> nonsquare or rank deficient matrices. Also for LinearSolve the built in
> cofactor method (Method ->CofactorExpansion) is not terribly efficient.
> Instead of using Cramer's rule is in effect inverting the matrix via
> cofactors and then multiplying the rhs by that inverse. This is probably a
> weakness in the implementation.
>
> I'll show what happens with what I think is the matrix and rhs you
> actually use (this is something you should have tried by now). I use mat
> as your A above.
>
> b = q/2*{0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0};
> rng = Complement[Range[16], {1, 2, 13}];
> mat2 = mat[[rng, rng]];
> b2 = b[[rng]];
>
> In[28]:= Timing[
> sol1 = LinearSolve[mat2, b2, Method -> "OneStepRowReduction"]]
>
> Out[28]= {0.271, {q/9, 0, (4 q)/3, 0, q/9, 0, (4 q)/3, 0, q/9, 0, 0,
> q/9, 0}}
>
> Your method also works well for this example.
>
> In[29]:= Timing[sol2 = myLinearSolve[mat2, b2]]
>
> Out[29]= {3.104, {q/9, 0, (4 q)/3, 0, q/9, 0, (4 q)/3, 0, q/9, 0, 0,
> q/9, 0}}
>
> In[36]:= Timing[sol3 = myLinearSolve[mat2, b2]]
>
> Out[36]= {1.882, {q/9, 0, (4 q)/3, 0, q/9, 0, (4 q)/3, 0, q/9, 0, 0,
> q/9, 0}}
>
> I'd define it a bit differently, though.
>
> myLinearSolve[A_, b_] :=
> Module[{i, detA = Simplify[Det[A]], atran = Transpose[A]},
> Table[Simplify[Det[ReplacePart[atran, i -> b]]/detA],
> {i,Length[A]}]]
>
> Daniel Lichtblau
> Wolfram Research
I reloaded 6.0 under Mac OS X 10.4.9 on a MacBookPro with the Intel
processor, and tried OneStepRowReduction. It works
fine and is 3-4 times faster than MyLinearSolve with your change.
The 13 x 13 system is solved in about 1.5 sec. Thanks for the tip!
There was an weird but unrelated font problem. Some matrices had
all greek letters in their symbols for example \[Beta]\[Pi]. When
printing the results some letters are replaced by black dots. (Does
not
happen on the G5 Macs). Also copy and paste of greeks
to the Find boxes doesnt work. Is a fix posted on the WRI website?