Services & Resources / Wolfram Forums
MathGroup Archive
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2000

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Speeding up Inverting matrices.

  • To: mathgroup at
  • Subject: [mg23035] Re: [mg23010] Speeding up Inverting matrices.
  • From: "Mark Harder" <harderm at>
  • Date: Thu, 13 Apr 2000 02:43:23 -0400 (EDT)
  • Sender: owner-wri-mathgroup at

    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:
Aprime = Table[Random[Real, {-10., 10.}], {i, 4}, {j, 4} ];
MatrixForm[Aprime ]
Det[Aprime ]

    RowBox[{"(", "\[NoBreak]", GridBox[{
            "8.079691866377622`", \(-1.426734910675819`\), \
\(-7.824867955066297`\), "3.727533573934952`"},
          {"1.0413965630255468`", \(-4.560122761941616`\),
            "5.829457548181766`", \(-7.928158183700707`\)},
            "4.943212521615639`", \(-1.4136769328239183`\), \
          {\(-4.4084483689155585`\), \(-5.980249047434835`\), \
\(-4.59332986628173`\), \(-6.125015284805416`\)}
          }], "\[NoBreak]", ")"}],
    (MatrixForm[ #]&)]\)

b = Table[Random[Real, {-2., 2.} ], {i, 4} ]

{1.84429, 0.102742, -0.97328, -1.07906}

ApInv = Inverse[Aprime];

    RowBox[{"(", "\[NoBreak]", GridBox[{
          {"0.06684139203781155`", "0.06031430797989468`",
            "0.021701857028159632`", \(-0.04399962320922808`\)},
          {\(-0.10579506911857166`\), \(-0.09527370237604071`\),
            "0.11340322473344215`", "0.02440994290443818`"},
            "0.07609778337870071`", \(-0.04306852356860005`\), \
            "0.06247971263886825`", \(-0.007457004396521146`\), \
\(-0.0940443542616278`\), \(-0.08695597076225173`\)}
          }], "\[NoBreak]", ")"}],
    (MatrixForm[ #]&)]\)

Needs["LinearAlgebra`MatrixManipulation`" ];

P = ZeroMatrix[4 ];
x0 = ApInv.b
xList = {x0};
dList = {1., 2., 3.};
      P[[2, 3]] = dList[[j]];
      AppendTo[xList, FixedPoint[(x0 - ApInv.P.#1) &, xList[[j]] ] ]),
    {j, 3} ];

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.30073}, {0.142184, -0.320066, 0.113109, 0.301514},
{0.13663, -0.311292,
               0.106102, 0.302201}}

(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

-----Original Message-----
From: David McGloin <dm11 at>
To: mathgroup at
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 McGloin
>Dept. of Physics
>Univ. of St. Andrews
>dm11 at

  • Prev by Date: Re: Re: making a column into a list
  • Next by Date: Re: Speeding up Inverting matrices.
  • Previous by thread: Re: Speeding up Inverting matrices.
  • Next by thread: Mathematica Wavelet-Explorer Question: Inverse Wavelet Transform on a subset of the coefficients