Re: Re: Complexity explosion in linear solve

• To: mathgroup at smc.vnet.net
• Subject: [mg80303] Re: [mg80264] Re: Complexity explosion in linear solve
• From: danl at wolfram.com
• Date: Fri, 17 Aug 2007 01:46:01 -0400 (EDT)
• References: <f9s3f8\$a8r\$1@smc.vnet.net><f9udg9\$ajv\$1@smc.vnet.net>

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

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

```

• Prev by Date: Re: Questions about InputField
• Next by Date: Re: Completely frivolous use of new v6 Mathematica
• Previous by thread: Re: Complexity explosion in linear solve
• Next by thread: Re: Complexity explosion in linear solve