Re: Speeding up Inverting matrices.
- To: mathgroup at smc.vnet.net
- Subject: [mg23035] Re: [mg23010] Speeding up Inverting matrices.
- From: "Mark Harder" <harderm at ucs.orst.edu>
- Date: Thu, 13 Apr 2000 02:43:23 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
David,
2 points, one practical, one algebraic:
Practical: 64Mb is not very much memory. On my PII 266MHz WinNT machine,
Mathematica can eat up 64Mb pretty easily. When this happens, the program
goes to virtual memory, i.e. putting & getting data to & from disk with an
access time in ms, instead of RAM in ns. I'd bet a 6-pack that your 16x16
problem is running in RAM +cache, but 24x24 gets bogged down with disk
access. Suggestion: buy yourself at least another 64Mb of RAM, better yet,
add 128Mb.
Algebraic: Since you will be solving your problem repeatedly with a
coefficient matrix that is part constant (& numeric) and part variable,
decompose the matrix--and therefore the problem -- into numeric and variable
parts; solve the numeric part once; use that solution to solve the whole
problem for succesive values of d, starting with d=0, of course. What I
mean is this:
Let A=A' + P, where A' is the numeric part of the coefficient matrix,
and P is a zero-matrix except for the one element that =d (& therefore when
d=0, P= all zeros).
So we have (A'+P).x=b.
=> A'.x +P.x=b
=> x = Inverse[A'].b - Inverse[A'].P.x ( eqn
* )
Begin with d=0, solve for x(0) = Inverse[A'].b
Now let d= d1 where d1 is the smallest value of interest to you,
& solve eqn * for x(d1) by a fixed-point method
(see Mathematica's FixedPoint[]function.), starting with x=x(o) as an
initial value.
Now, repeat that last step for d=d2 and initial x=x(d1), the idea being
that dn is always a step larger than dn-1.
Here's how I implemented this:
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
SeedRandom[];
Aprime = Table[Random[Real, {-10., 10.}], {i, 4}, {j, 4} ];
MatrixForm[Aprime ]
Det[Aprime ]
\!\(\*
TagBox[
RowBox[{"(", "\[NoBreak]", GridBox[{
{
"8.079691866377622`", \(-1.426734910675819`\), \
\(-7.824867955066297`\), "3.727533573934952`"},
{"1.0413965630255468`", \(-4.560122761941616`\),
"5.829457548181766`", \(-7.928158183700707`\)},
{"9.361455471297852`",
"4.943212521615639`", \(-1.4136769328239183`\), \
\(-1.864841482188929`\)},
{\(-4.4084483689155585`\), \(-5.980249047434835`\), \
\(-4.59332986628173`\), \(-6.125015284805416`\)}
}], "\[NoBreak]", ")"}],
(MatrixForm[ #]&)]\)
Out[14]=7302.79
In[17]:=
SeedRandom[];
b = Table[Random[Real, {-2., 2.} ], {i, 4} ]
Out[18]=
{1.84429, 0.102742, -0.97328, -1.07906}
In[19]:=
ApInv = Inverse[Aprime];
MatrixForm[ApInv]
Out[20]//MatrixForm=
\!\(\*
TagBox[
RowBox[{"(", "\[NoBreak]", GridBox[{
{"0.06684139203781155`", "0.06031430797989468`",
"0.021701857028159632`", \(-0.04399962320922808`\)},
{\(-0.10579506911857166`\), \(-0.09527370237604071`\),
"0.11340322473344215`", "0.02440994290443818`"},
{\(-0.009726094269161961`\),
"0.07609778337870071`", \(-0.04306852356860005`\), \
\(-0.09130648849920328`\)},
{
"0.06247971263886825`", \(-0.007457004396521146`\), \
\(-0.0940443542616278`\), \(-0.08695597076225173`\)}
}], "\[NoBreak]", ")"}],
(MatrixForm[ #]&)]\)
In[25]:=
Needs["LinearAlgebra`MatrixManipulation`" ];
In[36]:=
P = ZeroMatrix[4 ];
x0 = ApInv.b
xList = {x0};
dList = {1., 2., 3.};
Do[(
P[[2, 3]] = dList[[j]];
AppendTo[xList, FixedPoint[(x0 - ApInv.P.#1) &, xList[[j]] ] ]),
{j, 3} ];
xList
Out[37]={0.155828, -0.341618, 0.130324, 0.299827}
Out[41]={{0.155828, -0.341618, 0.130324, 0.299827}, {0.148524, -0.33008,
0.121108,
0.30073}, {0.142184, -0.320066, 0.113109, 0.301514},
{0.13663, -0.311292,
0.106102, 0.302201}}
In[35]:=
(Aprime + P).xList[[4]] - b
Out[35]=\!\({0.`, \(-4.440892098500626`*^-16\),
\(-2.220446049250313`*^-16\), 0.`}\)
<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
Sorry for the evil formatting produced by MatrixForm and whatnot, but the
last line shows that the solution for d=3. is within 10^-15 of zero. The
Do[] loop part took much less than one second on my system. The solution is
not symbolic, of course, which probably has saved lots of time, but you said
you wanted to plot elements of x vs. d, so this should work, just solve for
the x-vectors, then extract elements and plot them.
good luck,
mark harder
harderm at ucs.orst.edu
-----Original Message-----
From: David McGloin <dm11 at st-andrews.ac.uk>
To: mathgroup at smc.vnet.net
Subject: [mg23035] [mg23010] Speeding up Inverting matrices.
>I wish to solve the matrix equation Ax = b for x where A is a 24 x 24
>matrix and x and b are column matrices. Most of the values in the matrix
>are numbers (and many are equal to zero), but one remains unevaluated
>i.e element [1,10] may be 160 + d, where d is unevaluated. Currently
>we're using the command:
>
>x = {Inverse [A]. b}
>
>this works fine for the smaller matrices we use (8 x 8 and 16 x 16) but
>the calculation has now be runing for over 2 days (the smaller matrices
>may take many minutes if not seconds). The program is running on a PII
>350MHz with 64Mb of RAM. Does anyone have any ideas about how to
>optimise our calculation?
>
>Ultimately I want to extract arbitary elements of x and plot them
>against the unevaluated element.
>
>Thanks for any help!
>
>David
>
>--
>***************************************
>David McGloin
>Dept. of Physics
>Univ. of St. Andrews
>dm11 at st-and.ac.uk
>
>
>