Re: Speeding up Inverting matrices.
- To: mathgroup at smc.vnet.net
- Subject: [mg23042] Re: [mg23010] Speeding up Inverting matrices.
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Thu, 13 Apr 2000 02:43:27 -0400 (EDT)
- References: <200004120318.XAA04737@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
David McGloin wrote:
>
> 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
Below is a modification of a long response I made to a very similar
question that appeared on the news group sci.math.symbolic on 9/11/96.
Some of this also appeared in The Mathematica Journal Vol 7 (1), Winter
1997, pp26-27. It contains some ideas that should help. In short,
possibilities include the following.
(i) Use LinearSolve rather than Inverse, and try
Method->OneStepRowReduction.
(ii) Do several solutions for different values of your parameter in some
range on interest, and create a (vector-valued) interpolating function.
This is not mentioned explicitly below, but is hiddfen behind the
remarks about numeric solving.
(iii) Do a full-fledged polynomial interpolation. Details below.
Daniel Lichtblau
Wolfram Research
----------------------------------
The following problem, recast a bit for simplicity, appeared in the news
group sci.math.symbolic. A user has a sparse non-singular square matrix
mat.
The nonzero entries are low degree polynomials in a single variable. The
goal is to solve the linear system mat.sol == rhs, where rhs is, in this
case, a coordinate vector (that fact is irrelevant; we'd go through the
same steps for an arbitrary vector). The code below will generate a
matrix
and vector that suit these requirements. We show explicitly the matrix
and
vector we obtained.
Clear[randomPoly, randomSparsePoly, randomSparsePolyMatrix]
randomPoly[deg_, x_] := Table[Random[Integer, {-10,10}], {deg+1}] .
Table[x^j, {j, 0, deg}]
randomSparsePoly[x_] := Module[{r=Random[], deg},
Which [r < 0.6, Return[0],
r < 0.88, deg = 0,
r < 0.98, deg = 1,
r < 0.995, deg = 2,
True, deg = 3];
randomPoly[deg, x]
]
randomSparsePolyMatrix[dim_, x_] := Table[randomSparsePoly[x], {dim},
{dim}]
dim = 15;
mat = randomSparsePolyMatrix[dim, x]
rhs = Module[{r=Random[Integer, dim]}, Table[If [j==r, 1, 0], {j, dim}]]
mat = {{0, 0, -1, 0, 0, -6, -6 - 6*x - 3*x^2 - 6*x^3, 0, 0, 0, 0, 0, 0,
0, 1},
{0, -1, 8, 0, 6 + 6*x, 0, 0, 0, 0, 0, 9 - 7*x, 0, 0, 0, 0},
{0, 0, 0, 7 - 9*x, 0, 0, 0, 0, 0, 8, 0, 0, 8, 9, -1},
{4, 8, 0, 0, -6, 0, 0, 0, 0, 0, -5, -2 - x, -5, 0, -3},
{-5, 7, 4, -5 - 2*x, 0, 4, 1, 0, -4, -7 - x, 0, -6, 0, 0, -4},
{0, -4 + 3*x, 0, 7, 0, 0, 5, 5, -10, 0, 0, 7 - 6*x, 0, 2, 0},
{0, -4, 6 + 7*x, 0, 2 - 7*x, 9 + 9*x, -8, 0, 0, -4, 2, 5 + 5*x, 0, 0,
-9},
{0, 2 - 4*x, -5, 0, 0, 0, 0, 2, -5, -3 - 6*x, 0, 7 + x, 0, 6*x, 0},
{0, 5, 9, 6, -1, -7 - 9*x, 0, 0, 0, 7 - 10*x, 0, -2, 0, -10, 0},
{0, 3, 0, 0, 0, 6, 2, 0, 0, 0, 0, 0, -9, 10, 0},
{0, 0, -7, 0, 0, -5, 0, 0, 0, 0, -4, 0, -3 - 8*x, 0, 8},
{0, -3, -3, 0, 0, 0, 3, -2, -4 - 2*x, 5, 0, -3, -1, 0, 0},
{0, 0, 0, 0, -10 - 8*x, 4 + 10*x + 5*x^2, 0, 0, 0, 0, 0, 0, -2, 0,
0},
{0, 0, 0, 5, 0, 0, 4 + 2*x, 0, 0, -1, -3, 0, 0, 2, -2},
{7 - 6*x, 0, -9, 6 + 2*x + 3*x^2, 0, 0, 3 + 3*x, 0, 4 + 2*x, -6 -
6*x,
4 - 2*x, 5, 0, -2, 10}};
Apply[Plus, Last[RowReduce[mat /. x->2]]] (* shows mat is non-singular
*)
rhs = {0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
We now want to solve the system. Several comments can be made about this
process. First I should point out that, if you only want to later plug
in
numerical values for the variable (and especially if those values are
inexact), then it is often faster (and numerically MUCH more stable, if
using inexact values) to simply solve the system repeatedly with the
values
first plugged in. One can even reverse the process and use numerical
values
to obtain a symbolic solution. More on this later.
In regard to the problem as stated, the default linear algebra
technology
will suffer from severe intermediate swell. You can get an answer in
reasonable time if instead you do
Timing[ls1 = LinearSolve[mat, rhs, Method->OneStepRowReduction];]
This takes about 1 second on a Pentium Pro under Linux with 64 Meg RAM.
The
LeafCount is a bit over .7 Meg. Since we have in fact solved for one
column
of the inverse matrix, quite likely the full inverse computation
involves a
matrix around 15 times this big. Thus one sees why it is often better to
use
LinearSolve rather than multiplying left and right hand sides by
Inverse[mat]. (For inexact numeric matrices, it is preferable to use
LinearSolve for slightly diferent reasons. One gains a speed factor
around
2, and numeric stability is better).
I checked the solution by plugging x->1/3 into mat and into ls1, and
verified the identity mat.ls1==rhs.
In[4]:= Timing[ls1 = LinearSolve[mat, rhs,
Method->OneStepRowReduction];]
Out[4]= {0.91 Second, Null}
In[5]:= LeafCount[ls1]
Out[5]= 691818
In[6]:= m1 = mat /. x->1/3;
In[7]:= l1 = ls1 /. x->1/3;
In[8]:= m1 . l1 == rhs
Out[8]= True
That the one-step row reduction method was successful is good news. But
there is another method, not built in, that is typically better for
matrices of univariate polynomials. We know that the result can be
represented by rational functions in the variable x. We find the largest
exponents of each entry in each row, and sum these degrees, to get a
bound on the degree of Det[mat]. We can likewise bound it by summing
over Transpose[mat] (that is, sum the largest degrees that show up in
each column of mat).
In[12]:= Apply[Plus, Map[Max[Exponent[#, x]]&, mat]]
Out[12]= 18
In[13]:= Apply[Plus, Map[Max[Exponent[#, x]]&, Transpose[mat]]]
Out[13]= 17
We conclude that Det[mat] is a polynomial of degree no larger than 17.
We
can find it simply by interpolating from a table of 18 values. In the
problem as originally posed in the news group (unlike this random
example),
the matrix was stochastic; the variable was meant to take on values
between
zero and one. This is not relevant to the issue of interpolation, at
least
for exact matrices. We will do the simplest thing, that is, interpolate
at
small integer values. Provided they all give nonzero determinant, we can
reuse these values to interpolate cofactors later. Note that if our
polynomial entries contained rational coefficients, we would most likely
be better off multiplying our interpolation points by the least common
multiple of all the denominators. This would insure that we avoid
rational
arithmetic in the ensuing internal computations.
detlist = Table[{j, Det[mat /. x->j]}, {j, -8, 9}];
det = Expand[InterpolatingPolynomial[detlist, x]];
These take a fraction of a second on the Pentium Pro. Note that we can
actually get the result faster by direct means.
In[13]:= Timing[dd = Expand[Det[mat]];]
Out[13]= {0.79 Second, Null}
In[14]:= dd == det
Out[14]= True
While Det for symbolic matrices is far better than in previous versions,
I would still expect it to choke on matrices that are more dense. In
such
cases interpolation will typically work faster than the row reduction
used
internally. The method of interpolation is also desirable because, as we
will see, it lends itself to a MUCH more concise final result.
The entries in the result all arise as sums of integer products of
elements
in the inverse matrix (actually, the entries are just a particular
column of
the inverse). So they are all rational functions in x, with denominator
det
and numerators of degree no larger than 17. If we work instead with
(1/det)*mat, we have a matrix whose inverse contains only the polynomial
cofactors, because the denominator (that is, Det[(1/det)*mat]) is now 1.
We can thus again interpolate with 18 evaluation points, this time to
solve
the system over the integers. We can even reuse the values from detlist
provided none are zero, and we choose the same interpolation points. We
first check that nore of the determinants was zero.
In[17]:= Select[detlist, #[[2]]==0&]
Out[17]= {}
So we will reuse these points. We get rid of the first coordinates.
newdetlist = Map[#[[2]]&, detlist];
In[19]:= Timing[solnlist = Table[
LinearSolve[(1/newdetlist[[j+9]])*mat /. x->j, rhs], {j, -8, 9}];]
Out[19]= {1.13 Second, Null}
We use solnlist to interpolate the vector of polynomials vec that solves
(1/Det[mat]*mat).vec == rhs. Then we simply note that the solution sol
to
mat.sol == rhs is (obviously) sol==1/Det[mat]*vec.
abcissas = Table[j, {j, -8, 9}];
ordinateslists = Transpose[solnlist];
In[24]:= Timing[ls2 = (1/det)*Table[Expand[InterpolatingPolynomial[
Transpose[{abcissas,ordinateslists[[k]]}], x]], {k, 1, dim}];]
Out[24]= {0.18 Second, Null}
In[26]:= LeafCount[ls2]
Out[26]= 2041
l2 = ls2 /. x->1/3;
In[28]:= l2==l1 (* verify soln at x->1/3 *)
Out[28]= True
At the cost of some programming and a small factor of run time, we have
reduced the size of the result by a factor of more than 30. We can print
it
out even more concisely, by removing the denominator common to all
entries.
So we do
ls2numerators = Table[Expand[InterpolatingPolynomial[
Transpose[{abcissas,ordinateslists[[k]]}], x]], {k, 1,
dim}];
and then the result is (1/det)*ls2numerators, where det is
In[30]:= det//InputForm
Out[30]//InputForm=
-107759679244830 - 208362055969602*x - 46531268980354*x^2 +
40007106260995*x^3 - 149435262008162*x^4 - 261723072465380*x^5 -
483027101696815*x^6 - 498458053530926*x^7 - 237657126195021*x^8 -
101477072853478*x^9 - 52625582899214*x^10 - 15963313573744*x^11 -
1422536423784*x^12 + 669326184576*x^13 + 93960518400*x^14
and ls2numerators is
In[31]:= ls2numerators // InputForm
Out[31]//InputForm=
{-23002184559950 - 41160789773856*x - 7330880812072*x^2 +
2698130484509*x^3 - 13837854357467*x^4 - 6308052006349*x^5 -
36164144063000*x^6 - 48938362409374*x^7 - 21729624839460*x^8 -
9022519455968*x^9 - 4866877301276*x^10 - 1046518468128*x^11 -
10153048320*x^12 + 32043247488*x^13,
-1225128464878 + 6345075096038*x + 3840573681244*x^2 -
11464299489021*x^3 - 14451650302968*x^4 - 26524953405181*x^5 -
53459175942975*x^6 - 45471988115225*x^7 - 28727636644662*x^8 -
20142181144038*x^9 - 8520323161364*x^10 - 2335324449648*x^11 -
652951266048*x^12 - 45637966080*x^13,
-16298917327814 - 199435673198*x + 17490247565024*x^2 -
7445688930707*x^3 - 2708575819656*x^4 + 20416177360623*x^5 +
1864805112381*x^6 + 5767096509489*x^7 + 13613767073874*x^8 +
2848922018194*x^9 + 1193299207032*x^10 + 1634444790312*x^11 +
425680728768*x^12 + 78773990400*x^13,
9664348949460 + 8457957756134*x - 15198670562924*x^2 -
5609551195968*x^3 + 6970114925463*x^4 - 6476031018256*x^5 +
2450450778082*x^6 + 5603184747422*x^7 - 2708303333636*x^8 -
801408326872*x^9 + 733761562944*x^10 + 81214413696*x^11 +
13835743488*x^12, 2635066221642 + 2432665519286*x -
3260610658984*x^2 -
10877576978307*x^3 - 10832110036238*x^4 - 11386363866294*x^5 -
13727575761743*x^6 - 4122882246920*x^7 - 3513783074830*x^8 -
4801055186688*x^9 - 450603176680*x^10 + 280567594464*x^11 -
111966160320*x^12 + 11659334400*x^13,
3344726000456 + 3274575263812*x - 12539205219896*x^2 -
12150937612654*x^3 - 5346097485964*x^4 - 19222016066626*x^5 -
6354271717048*x^6 - 1268996265812*x^7 - 7388818894560*x^8 -
1245309502480*x^9 + 566236685184*x^10 - 193137057792*x^11 +
18654935040*x^12, -2012532849102 - 1353860010976*x +
12972071617870*x^2 + 1739120483917*x^3 - 2460709186561*x^4 +
5716222565292*x^5 + 1146148392754*x^6 - 1029835614318*x^7 -
119872828022*x^8 - 66455461248*x^9 - 13345528896*x^10,
25837928430684 - 47087786062952*x - 9923657341942*x^2 +
52249052678936*x^3 - 24731514322960*x^4 - 17275012889349*x^5 +
57719775018819*x^6 + 16308304769992*x^7 + 10606140768985*x^8 +
25788175688534*x^9 + 2491408362450*x^10 - 1579661919372*x^11 +
673706727360*x^12 - 349167262464*x^13 - 85369842432*x^14,
12331734417106 + 510578584906*x - 8815694612450*x^2 +
6218051858851*x^3 +
5121784131192*x^4 - 3054037533779*x^5 + 23647136406357*x^6 +
21689537131148*x^7 + 6575214568109*x^8 + 9157320946002*x^9 +
4841797800720*x^10 + 505083830112*x^11 + 269664318720*x^12 +
85369842432*x^13, 4949695872562 + 6279578576294*x +
6770618746088*x^2 -
1993010863407*x^3 - 12177381317020*x^4 - 15650253772834*x^5 -
14406578658385*x^6 - 4661835023730*x^7 - 89299023562*x^8 -
1092894877808*x^9 - 99237076944*x^10 + 311437164672*x^11 +
35222947200*x^12, 12595090314198 + 7299978026160*x -
8890432201960*x^2 +
7855256929629*x^3 + 21384652983217*x^4 + 10349893556394*x^5 +
17195030646773*x^6 + 15095466932028*x^7 + 1538943024338*x^8 +
1969785938876*x^9 + 3025748326080*x^10 + 754396555968*x^11 +
23418910080*x^12 + 9993715200*x^13,
-7744632409638 + 25266665161728*x + 18879036236978*x^2 -
16125171845563*x^3 + 1625156872353*x^4 - 3654020421635*x^5 -
36576941003048*x^6 - 22371120339943*x^7 - 14469621637748*x^8 -
16864443927672*x^9 - 6370305059972*x^10 - 1538664190104*x^11 -
834362919936*x^12 - 93960518400*x^13,
-6485879107298 + 569188046906*x + 6228672098184*x^2 -
11381135637787*x^3 - 5124037990520*x^4 + 4708395881715*x^5 -
8000533208245*x^6 - 6839677001857*x^7 - 2947854048970*x^8 -
4546795909530*x^9 - 2108884748280*x^10 - 268789711464*x^11 -
75223282176*x^12, -7074081687558 - 3085226442688*x +
9382741592356*x^2 +
139006243507*x^3 + 3423661228313*x^4 + 22585007442015*x^5 +
12420606247150*x^6 + 8453252015247*x^7 + 10422488251666*x^8 +
2711014818372*x^9 + 321027769626*x^10 + 574568829252*x^11 +
116991464832*x^12 + 13691389824*x^13,
-8305758419690 - 750341250794*x + 5926687339706*x^2 -
8220939123449*x^3 -
8321638163550*x^4 + 7666952392224*x^5 + 7965995902188*x^6 +
1235408161743*x^7 + 118383622488*x^8 - 1953521218736*x^9 -
2426718792702*x^10 - 523054081692*x^11 + 98840984832*x^12 -
1299182976*x^13}
- References:
- Speeding up Inverting matrices.
- From: David McGloin <dm11@st-andrews.ac.uk>
- Speeding up Inverting matrices.