Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

Mathematica 45x45 Symbolic Matrix Inversion (long)

825 views
Skip to first unread message

Michael R. Henry

unread,
Sep 10, 1996, 3:00:00 AM9/10/96
to

Hi,
I'm using Mathematica to compute the inverse of
a 45x45 symbolic matrix, and it's taking too long (4 days so
far). Anybody have any ideas? Here are the details:

I want to symbolically compute the stationary probability vector
for a 45-state Markov Chain, where the state transition
matrix has symbolic entries of the form {0, 1, 1-p, p/2}.
(Mostly zeros, as you can see from the matrix stm below).

I need the answer in terms of the variable 'p'. Basically,
I need to solve the following, where pi is 1x45 vector,
and the state transition matrix stm is a 45x45 matrix:

pi = pi . stm, subject to Sum[pi[[i]],{i,1,45}] = 1.

I've made the following definitions:

id=IdentityMatrix[45]; u=Table[1,{45}]; z=Table[0,{45}];

So my equation becomes
pi . (stm - id + 1) = u
pi = u . Inverse[stm-id+1]

For better efficiency, the manual says to compute this using
the LinearSolve[m,x] function, which finds the column-vector x
which satisfies m.x=b, for matrix m and column-vector b.
I converted my equation to use this form, getting
LinearSolve[Transpose[stm-id+1],u].

I try this out by setting p=1/2 first, and get a valid answer
for pi (in perhaps 20 seconds), since I write:
In[19]:= pi==pi.stm
Out[19]= True
In[20]:= pi.u==1
Out[20]= True

However, if I leave p undefined to get a solution in terms
of p, my Pentium 133 running the Linux version of Mma 2.2.4 is still
trying to compute the answer (Since Friday morning, at least
4 days so far).

My question is: Is this operation inherently so hard that a 45x45
symbolic matrix inversion takes many days to solve? I tried this problem
using MuPAD on a 486-66 (roughly 4x slower than my Pentium 133),
and MuPAD got an answer in eight days. I began looking for another
computer algebra package that could solve this *significantly* faster,
and my first attempt at using Mma has been disappointing. I would
have expected Mma to take roughly 2 days or less on the Pentium were it
faster than MuPAD.

Also, am I attempting this calculation in the best manner?
I tried the direct matrix inversion for a while, but only gave
Mma a few hours before going on to the LinearSolve[] attempt.
Perhaps there is some better Mma idiom to use.
Or, maybe I should look for some other package to use.

I freely offer this problem as a benchmark for anyone to show
off his computer algebra system and/or his barn-burning fire-breathing
monster computer ;-)

I eventually did get an answer out of MuPAD as stated earlier, but
I would really like to be able to repeat this calculation with
variations on stm, and cannot afford to wait 8 days or more between
iterations.
Anyone have any ideas?

TIA,
Michael Henry
(I apologize for the length of this post, but such is the nature
of 45x45 matrices :)

Mathematica-compatible definition of state transition matrix:
stm = {
{1-p, p/2, p/2, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 1-p, p/2, p/2,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
1-p, p/2, p/2, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 1-p, p/2, p/2, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1-p, p/2,
p/2, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 1-p, p/2, p/2, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1-p, p/2, p/2, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1-p,
p/2, p/2, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1-p, p/2, p/2, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
1, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 1-p, p/2, p/2, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1-p, p/2, p/2},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, p/2, p/2, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1-p, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, p/2,
p/2, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1-p, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, p/2, p/2, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1-p,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, p/2, p/2, 0,
0, 0, 0, 0, 0, 0, 0, 0,
1-p, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
p/2, p/2, 0, 0, 0, 0, 0, 0,
0, 1-p, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, p/2, p/2, 0, 0, 0,
0, 0, 1-p, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0}};


Carlos A. Felippa

unread,
Sep 11, 1996, 3:00:00 AM9/11/96
to

After 0.05 seconds, Mma says that matrix stm is singular for any
value of p.


Michael R. Henry

unread,
Sep 11, 1996, 3:00:00 AM9/11/96
to

In article <515agi$l...@lace.colorado.edu>,

Carlos A. Felippa <car...@mars.Colorado.EDU> wrote:

>After 0.05 seconds, Mma says that matrix stm is singular for any
>value of p.

Yes, matrix stm is a stochastic matrix, and will always be singular.
I didn't mean to be unclear, but the matrix I need to invert is
not stm, but rather (stm-IdentityMatrix[45] + 1), where the
value '1' stands for a 45x45 matrix of all ones in Mma notation.
This matrix does have an inverse. Actually I don't need the
inverse directly, but the stationary probability vector pi
defined below (as quoted from my original article):

> id=IdentityMatrix[45]; u=Table[1,{45}]; z=Table[0,{45}];
>
> So my equation becomes

> pi = u . Inverse[stm-id+1]

And Mma is still trying (5 days so far)

Thanks,
Michael Henry

Daniel Lichtblau

unread,
Sep 11, 1996, 3:00:00 AM9/11/96
to

I'll make several comments about this problem. First I should point out
that, if you only want to plug in numerical values for p (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 quite
quickly if instead you do

Timing[ls1 = LinearSolve[Transpose[stm-id+1], u,
Method->OneStepRowReduction];]

This takes a bit under 10 sec on a Pentium Pro under linux with 64Meg RAM.

SHAMELESS PLUG FOR NEXT RELEASE
Or wait a few weeks, get Mathematica version 3, and it will take about 2
sec.
END SHAMELESS PLUG


The result will be quite large. Quite likely the intermediate computation
involved a matrix around 45 times this big, so all that RAM probably
played a significant role in keeping down the time.

In[4]:= LeafCount[ls1]
Out[4]= 1209698

MORE ADVERTISEMENT HYPE: LeafCount in version 3 is under a meg.


I checked the solution by plugging p->1/3 into stm and it, and verified
the identities pi==pi.stm and pi.u==1.

In[15]:= l1 = ls1 /. p->1/3;

In[18]:= stm1 = stm /. p->1/3;

In[20]:= l1 == l1 . stm1
Out[20]= True

In[21]:= l1 . u
Out[21]= 1


Next, observe that while the entries are still simple, you do not exactly
have a sparse matrix once you go adding 1 to every entry. Observing that
row 39 of (stm-id) is all zeros except for 1 in the last column and -1 in
column 39, and also that there are other rows with nonzero elements in
these columns, we guess that row 39 is linearly dependent on the remaining
rows. Well, it turns out that this was just inspiration, because actually
it is column 39 that will be of interest. On the other hand, as the rank
is 44 (plug in any numerical value for p and take Eigenvalues to verify
this), we expect that almost any row will be dependent on the rest. Given
the sparsity of the matrix, maybe this is a bad guess. Who cares. We'll
just check as below. We remove row 39 of the transposed matrix to form a
new matrix, plug in a value for p, row reduce, and check that the last row
has a pivot. I tried the substitution p->1 first.

In[108]:= tmp = Drop[Transpose[(stm-id)/.p->1], {39}];

In[109]:= Last[RowReduce[tmp]]
Out[109]= {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0,
> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1}

Note that the choice p->0 is not so lucky.

In[112]:= Drop[Transpose[(stm-id)/.p->0], {39}];

In[113]:= Last[RowReduce[tmp]]
Out[113]= {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0,
> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}

No matter. The result from p->1 substitution suffices to show that row 39
of the transposed matrix depends linearly on the rest. Thus if we simply
add 1 to each entry of row 39 of Transpose[stm - id], we can then have a
right-hand-side vector of zeroes, except for a one in the 39th slot. This
recovers the sparsity.

stm2 = Transpose[stm - id];
stm2[[39]] += 1;
u = Table[0, {45}];
u[[39]] = 1;


Now we will see that things improve in terms of speed and size of result.

In[114]:= Timing[ls2 = LinearSolve[stm2, u, Method->OneStepRowReduction];]
Out[114]= {5.85 Second, Null}

In[115]:= LeafCount[ls2]
Out[115]= 285045

(Version 3 does it in just over a second, result has LeafCount around
155K).

Again we check the result at p->1/3.

In[116]:= l2 = ls2 /. p->1/3;

In[117]:= l2 == l1 (* verify soln at p->1/3 *)
Out[117]= True

That one-step row reduction metod was successful is good news. But there
is another method (not built in, I'm afraid) that might be better still
for e.g. denser matrices of univariate polynomials. We know that the
result can be represented by rational functions in the variable p. We find
the largest exponents of each entry in each row, and sum these degrees, to
get a bound or the degree of Det[stm2]. We can likewise bound it by
summing over Transpose[stm2] (that is, sum the largest degrees that show
up in each column of stm2).

In[130]:= Apply[Plus, Map[Max[Exponent[#, p]]&, stm2]]
Out[130]= 39

In[131]:= Apply[Plus, Map[Max[Exponent[#, p]]&, Transpose[stm2]]]
Out[131]= 17

We conclude that Det[stm2] is a polynomial of degree no larger than 17. We
can find it simply by interpolating from a table of 18 values. Notice that
even though p will take values between zero and one, there is no reason to
so restrict the interpolation values. If we use even integer values, we
can avoid all rational arithmetic (because the largest denominator in stm2
is 2).

In[135]:= Max[Denominator[Flatten[Together[stm2]]]]
Out[135]= 2

detlist = Table[{j, Det[stm2 /. p->j]}, {j, -16, 18, 2}];
det = Expand[InterpolatingPolynomial[detlist, p]];

These take a total of maybe 3 seconds on the pentium pro. Note that in
version 3 we can actually get the result faster by direct means.

In[31]:= Timing[dd = Expand[Det[stm2]];]
Out[31]= {1.39 Second, Null}

In[32]:= det == dd
Out[32]= 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 row reduction. 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 column 39
of the inverse). So they are all rational functions in p, with denominator
det and numerators of degree no larger than 17. If we work instead with
(1/det)*stm2, we have a matrix whose inverse contains only the polynomial
cofactors, because the denominator is now 1. We can thus (again)
interpolate with 18 evaluation points. We could even use the values from
detlist provided we choose the same interpolation points. Except this will
fail, because the value of the determinant when p->0 is zero and hence at
that point the system is singular. Not a surprise, given what happened in
Out[113]. So we compute at even integers from 2 to 36.


newdetlist = Table[Det[stm2 /. p->j], {j, 2, 36, 2}];

In[157]:= Timing[solnlist = Table[LinearSolve[(1/newdetlist[[j/2]])*stm2
/. p->j, u],
{j, 2, 36, 2}];]

Out[157]= {31.94 Second, Null}

(Takes 8.5 sec in version 3.)

We use solnlist to interpolate the vector of polynomials x that solves
1/Det[stm2]*stm2 . x == u. Then we simply note that the solution y to
stm2 . y == u is (obviously) y==1/Det[stm2]*x.

abcissas = Table[j, {j, 2, 36, 2}];
ordinateslists = Transpose[solnlist];

In[186]:= Timing[ls3 =
(1/det)*Table[Expand[InterpolatingPolynomial[
Transpose[{abcissas,ordinateslists[[k]]}], p]],
{k, 1, 45}];]
Out[186]= {0.65 Second, Null}

(* This one, sob, slowed to a crawling 1.2 sec in version 3. Something for
us to improve).


In[187]:= LeafCount[ls3]
Out[187]= 6498

In[188]:= l3 = ls3 /. p->1/3;

In[190]:= l3==l1 (* verify soln at p->1/3 *)
Out[190]= True


If we instead do

ls3numerators = Table[Expand[InterpolatingPolynomial[
Transpose[{abcissas,ordinateslists[[k]]}], p]], {k, 1,
45}];

then the result is (1/det)*ls3numerators, where det is

In[194]:= det//InputForm
Out[194]//InputForm=
(23*p^5)/8 + (209*p^6)/32 - (121*p^7)/64 - (165*p^8)/16 - (151*p^9)/128
+
(401*p^10)/64 + (107*p^11)/64 - (203*p^12)/128 - (19*p^13)/32 + p^14/8
+
p^15/16

and ls3numerators is

In[195]:= ls3numerators // InputForm
Out[195]//InputForm=
{p^5/16 + p^6/8 - p^7/4 - (21*p^8)/128 + (11*p^9)/32 - p^10/32 -
(13*p^11)/128 + p^12/32, p^6/32 + p^7/16 - p^8/8 - (21*p^9)/256 +
(11*p^10)/64 - p^11/64 - (13*p^12)/256 + p^13/64,
(5*p^5)/32 + (5*p^6)/64 - (35*p^7)/128 - (3*p^8)/16 + (13*p^9)/64 +
(39*p^10)/256 - (11*p^11)/128 - (9*p^12)/256 + p^13/64,
p^6/16 + (5*p^7)/32 - (11*p^8)/64 - (59*p^9)/256 + (7*p^10)/32 +
(5*p^11)/128 - (17*p^12)/256 + (3*p^13)/256,
p^7/64 + p^8/32 - p^9/16 - (21*p^10)/512 + (11*p^11)/128 - p^12/128 -
(13*p^13)/512 + p^14/128, p^6/4 + (5*p^7)/64 - (3*p^8)/8 - (3*p^9)/32
+
(111*p^10)/512 + (13*p^11)/256 - p^12/16 - (5*p^13)/512 + p^14/128,
(5*p^5)/32 - (3*p^6)/64 - (29*p^7)/128 + (13*p^8)/128 + (15*p^9)/128 -
(11*p^10)/128 - (3*p^11)/256 + (5*p^12)/256 - p^13/256,
(5*p^6)/64 + (5*p^7)/128 - (35*p^8)/256 - (3*p^9)/32 + (13*p^10)/128 +
(39*p^11)/512 - (11*p^12)/256 - (9*p^13)/512 + p^14/128,
(21*p^6)/64 + (13*p^7)/128 - (139*p^8)/256 - p^9/8 + (23*p^10)/64 +
(21*p^11)/512 - (25*p^12)/256 - p^13/512 + p^14/128,
p^6/16 + p^7/8 - p^8/4 - (21*p^9)/128 + (11*p^10)/32 - p^11/32 -
(13*p^12)/128 + p^13/32, p^6/32 + p^7/8 + p^8/64 - (35*p^9)/128 -
(9*p^10)/256 + (29*p^11)/128 - p^12/32 - (7*p^13)/128 + p^14/64,
p^6/32 + p^7/8 + p^8/64 - (35*p^9)/128 - (9*p^10)/256 + (29*p^11)/128 -
p^12/32 - (7*p^13)/128 + p^14/64,
p^7/32 + (5*p^8)/64 - (27*p^9)/256 - (27*p^10)/256 + (19*p^11)/128 +
p^12/256 - (3*p^13)/64 + (3*p^14)/256,
p^7/64 + p^8/16 - (3*p^9)/256 - (65*p^10)/512 + (11*p^11)/512 +
(25*p^12)/256 - (15*p^13)/512 - (11*p^14)/512 + p^15/128,
p^7/64 + p^8/16 - (3*p^9)/256 - (65*p^10)/512 + (11*p^11)/512 +
(25*p^12)/256 - (15*p^13)/512 - (11*p^14)/512 + p^15/128,
p^6/4 - (5*p^7)/32 - (25*p^8)/64 + (69*p^9)/256 + (47*p^10)/256 -
(37*p^11)/256 - p^12/64 + (3*p^13)/128 - p^14/256,
p^6/64 + (31*p^7)/128 + (11*p^8)/256 - (105*p^9)/256 - (13*p^10)/256 +
(133*p^11)/512 + (13*p^12)/512 - (37*p^13)/512 - (3*p^14)/512 +
p^15/128,
p^6/64 + (31*p^7)/128 + (11*p^8)/256 - (105*p^9)/256 - (13*p^10)/256 +
(133*p^11)/512 + (13*p^12)/512 - (37*p^13)/512 - (3*p^14)/512 +
p^15/128,
(5*p^5)/32 + (3*p^6)/64 - (43*p^7)/128 - p^8/16 + (73*p^9)/256 -
(5*p^10)/256 - (9*p^11)/128 + p^12/64,
p^6/4 + p^7/16 - (13*p^8)/32 - p^9/32 + (33*p^10)/128 - (9*p^11)/256 -
(7*p^12)/128 + p^13/64, p^6/4 + p^7/16 - (13*p^8)/32 - p^9/32 +
(33*p^10)/128 - (9*p^11)/256 - (7*p^12)/128 + p^13/64,
(3*p^6)/32 + (13*p^7)/64 - (17*p^8)/128 - (47*p^9)/128 + (37*p^10)/256
+
(15*p^11)/64 - (3*p^12)/32 - (3*p^13)/64 + (5*p^14)/256,
p^7/8 + p^8/16 - (65*p^9)/256 - (47*p^10)/512 + (99*p^11)/512 +
(13*p^12)/256 - (33*p^13)/512 - (5*p^14)/512 + p^15/128,
p^7/8 + p^8/16 - (65*p^9)/256 - (47*p^10)/512 + (99*p^11)/512 +
(13*p^12)/256 - (33*p^13)/512 - (5*p^14)/512 + p^15/128,
(11*p^6)/32 + p^7/64 - (77*p^8)/128 + p^9/128 + (111*p^10)/256 -
(15*p^11)/256 - (29*p^12)/256 + (3*p^13)/128 + p^14/256,
(11*p^6)/64 + (53*p^7)/128 - (69*p^8)/256 - (81*p^9)/128 +
(41*p^10)/256 + (197*p^11)/512 - (23*p^12)/512 - (51*p^13)/512 +
(3*p^14)/512 + p^15/128, (11*p^6)/64 + (53*p^7)/128 - (69*p^8)/256 -
(81*p^9)/128 + (41*p^10)/256 + (197*p^11)/512 - (23*p^12)/512 -
(51*p^13)/512 + (3*p^14)/512 + p^15/128,
p^5/16 + (7*p^6)/32 - (7*p^7)/64 - (49*p^8)/128 + (43*p^9)/256 +
(7*p^10)/32 - (13*p^11)/128 - (9*p^12)/256 + p^13/64,
p^6/32 + (13*p^7)/64 - (13*p^8)/128 - (25*p^9)/64 + (47*p^10)/256 +
(7*p^11)/32 - (7*p^12)/64 - p^13/32 + p^14/64,
p^5/32 + (13*p^6)/64 - (13*p^7)/128 - (25*p^8)/64 + (47*p^9)/256 +
(7*p^10)/32 - (7*p^11)/64 - p^12/32 + p^13/64,
(11*p^5)/32 + (11*p^6)/64 - (95*p^7)/128 - (23*p^8)/128 + (73*p^9)/128
+
(13*p^10)/256 - (11*p^11)/64 - p^12/256 + p^13/64,
(11*p^6)/32 + (21*p^7)/64 - (93*p^8)/128 - (45*p^9)/128 +
(137*p^10)/256 + p^11/8 - (5*p^12)/32 - p^13/64 + p^14/64,
(11*p^5)/32 + (21*p^6)/64 - (93*p^7)/128 - (45*p^8)/128 + (137*p^9)/256
+
p^10/8 - (5*p^11)/32 - p^12/64 + p^13/64,
p^5/16 + (7*p^6)/32 - (7*p^7)/64 - (49*p^8)/128 + (43*p^9)/256 +
(7*p^10)/32 - (13*p^11)/128 - (9*p^12)/256 + p^13/64,
p^6/32 + (13*p^7)/64 - (13*p^8)/128 - (25*p^9)/64 + (47*p^10)/256 +
(7*p^11)/32 - (7*p^12)/64 - p^13/32 + p^14/64,
p^5/32 + (13*p^6)/64 - (13*p^7)/128 - (25*p^8)/64 + (47*p^9)/256 +
(7*p^10)/32 - (7*p^11)/64 - p^12/32 + p^13/64,
(11*p^5)/32 + (11*p^6)/64 - (95*p^7)/128 - (23*p^8)/128 + (73*p^9)/128
+
(13*p^10)/256 - (11*p^11)/64 - p^12/256 + p^13/64,
(11*p^6)/32 + (21*p^7)/64 - (93*p^8)/128 - (45*p^9)/128 +
(137*p^10)/256 + p^11/8 - (5*p^12)/32 - p^13/64 + p^14/64,
(11*p^5)/32 + (21*p^6)/64 - (93*p^7)/128 - (45*p^8)/128 + (137*p^9)/256
+
p^10/8 - (5*p^11)/32 - p^12/64 + p^13/64,
p^5/16 + (3*p^6)/16 - p^7/8 - (3*p^8)/8 + (41*p^9)/256 + (15*p^10)/64 -
(13*p^11)/128 - (11*p^12)/256 + (5*p^13)/256,
p^6/32 + (7*p^7)/64 - (7*p^8)/128 - (49*p^9)/256 + (43*p^10)/512 +
(7*p^11)/64 - (13*p^12)/256 - (9*p^13)/512 + p^14/128,
p^5/32 + (15*p^6)/64 + p^7/128 - (57*p^8)/128 - p^9/128 +
(155*p^10)/512 - (21*p^12)/256 - p^13/512 + p^14/128,
(11*p^5)/32 + (11*p^6)/64 - (75*p^7)/128 - (21*p^8)/128 + (51*p^9)/128
+
p^10/64 - (25*p^11)/256 + (3*p^12)/256 + p^13/256,
(11*p^6)/64 + (11*p^7)/128 - (95*p^8)/256 - (23*p^9)/256 +
(73*p^10)/256 + (13*p^11)/512 - (11*p^12)/128 - p^13/512 + p^14/128,
(11*p^5)/32 + p^6/2 - (41*p^7)/64 - (185*p^8)/256 + (57*p^9)/128 +
(105*p^10)/256 - (67*p^11)/512 - (13*p^12)/128 + (7*p^13)/512 +
p^14/128}


Daniel Lichtblau
Wolfram Research, Inc.
da...@wolfram.com

In article <513ru0$2...@vdal4.eng.umd.edu> mrh...@Glue.umd.edu (Michael R.

> < matrix deleted for brevity (danl)>


0 new messages