Symbolic inverse of a 12x12 matrix
Madhusudan Singh
I am solving a problem using f95 that involves solving a 12x12 coupled
linear system. I have written out the matrix A (A x = f), and instead of
calculating the inverse repeatedly for different values of the input
parameters, I am trying to calculate the symbolic inverse of the matrix,
post multiply it with the RHS and code the resulting 12 equations for x_i
into fortran.
I believe this is the sensible way to do it, since this system needs to be
solved lots of times (its part of a half-infinite integral) and inverting
such a system numerically would not only be time consuming, but also
obscure possible stability problems.
Now, I am using Mathematica 5.1 (multiprocessor Sun machine) to do it, but
the inverse has been running for the past 40 minutes without any result. A
prior attempt to run it with Mathematica 4.1 (single processor Pentium M
Linux machine) resulted in an aborted process.
I am wondering if such long inversion times for this medium sized system are
typical for Mathematica. Could I make it faster by compiling the notebook
in some fashion ? I have sparingly used maxima (formerly known as macsyma)
in the past, and have heard that it generally gives faster results than
Mathematica does. Is it reliable enough to make it worthwhile to try that ?
Mathematica makes it possible to render the result in a human readable
format using MatrixForm. Is such a device available in maxima (if it is
suitable for this task) ?
Though I have not heard good things about Matlab's symbolic computation
package, is that also a tool worth trying ?
Thanks.
Gordon Sande
A 12x12 matrix would have 144 entries so at the naive level an
expression for the inverse would involve the 144 unknowns.
To cover all the cases the will be 12! (12 factorial) terms in
the expression. Each term will be the product of 12 variables.
Remember this will just be a determinant expansion.
The long execution time for the symbolic inversion is just
the computers subtle way of trying to point out that it
just might no be the most sensible thing to do.
When (and even if) you get the expression you will then face the
problem of evaluating it. The expression may not be a very stable
representation of what you are trying to do.
Sometimes the "simplicity" of a symbolic solution requires a rather
special view to convince yourself that it is in fact simple.
Perhaps you should try 3x3, then 4x4 and so on to see just what
you are getting into.
We all await your report.
Madhusudan Singh
> A 12x12 matrix would have 144 entries so at the naive level an
> expression for the inverse would involve the 144 unknowns.
>
No. Nearly 2/3 of the entries are zero, and the terms that remain are all
functions of 4-6 variables.
Unfortunately, the matrix is not tridiagonal, or UT, or LT or any simple
form like that.
> To cover all the cases the will be 12! (12 factorial) terms in
> the expression. Each term will be the product of 12 variables.
> Remember this will just be a determinant expansion.
Please read above. I should have made this clearer earlier.
>
> The long execution time for the symbolic inversion is just
> the computers subtle way of trying to point out that it
> just might no be the most sensible thing to do.
>
I think it is, in light of the above. Awaiting further suggestions.
Alec Mihailovs
news:425db863$0$79461$1472...@news.sunsite.dk...
> No. Nearly 2/3 of the entries are zero, and the terms that remain are all
> functions of 4-6 variables.
>
> Unfortunately, the matrix is not tridiagonal, or UT, or LT or any simple
> form like that.
Sometimes LU-decomposition works pretty good for such problems. Did you try
it?
Alec
Julian V. Noble
The symbolic inverse of a 12x12 matrix is a fairly tall order. It might
be worth your while to express the result as the inverse of --say-- 6x6
matrices, using identities like those of Strassen's method.
Also you should first try your Mathematica program on 4x4 or 6x6 matrices
to see how long they take.
(Just my 2¢ worth.)
--
Julian V. Noble
Professor Emeritus of Physics
j...@lessspamformother.virginia.edu
^^^^^^^^^^^^^^^^^^
http://galileo.phys.virginia.edu/~jvn/
"For there was never yet philosopher that could endure the
toothache patiently."
-- Wm. Shakespeare, Much Ado about Nothing. Act v. Sc. 1.
Jon Harrop
> I think it is, in light of the above. Awaiting further suggestions.
Plot graphs the time taken to analyse similar but smaller matrices and guess
how long this is going to take?
--
Dr Jon D Harrop, Flying Frog Consultancy
http://www.ffconsultancy.com
Madhusudan Singh
Hi
Out of a hunch, I just exported the notebook to text, and used "math" (the
text version of Mathematica) to run it through. Less that two minutes
later, I had an answer. The problem is that the answer (in nohup.out) is so
ill formatted (I used MatrixForm) that it is impossible to make sense of.
OTOH, the graphical Notebook is still chugging away. Has anyone seen this
odd behaviour before ?
Thanks.
Daniel Lichtblau
Unless the structure is in some way VERY nice, a 12x12 symbolic
inversion is not to be in your near future. In Mathematica or
otherwise.
Things you might try:
(1) Inverse[matrix, Method->OneStepRowReduction]
(2) Have a look at
The method indicated there, using interpolation, will be difficult to
adapt because you have multivariate entities. In general this will
involve need for a sparse polynomial interpolator, which is not
terribly easy to implement. And i'm not certain offhand if you can
compute the needed information for a determinant computation in the way
it is found for e.g. Zippel's sparse multivariate gcd algorithm.
(3) Work with the numeric matrices. I put this last because I figured
only after you saw (1) and (2) would you change your mind about the
numeric approach. It is almost certainly the appropriate method. Note
that you can check conditioning for each system you solve, hence get a
fair assessment of numerical problems. In contrast, a huge symbolic
result involving rational functions might have horrible numeric
instabilities in evaluating at finite precision, either from
cancellation error or from a denominator expression with vanishing set
that comes close to many of the points you might find of interest.
Daniel Lichtblau
Wolfram Research
Madhusudan Singh
>> I think it is, in light of the above. Awaiting further suggestions.
>
>
> Unless the structure is in some way VERY nice, a 12x12 symbolic
> inversion is not to be in your near future. In Mathematica or
> otherwise.
>
I know it is. Precisely 1.5 minutes in my future if I start it now. Please
see my latest post on this thread (interesting developments) which
describes a rather unexpected occurence. Not that I can read the answer
that Mathematica is giving me.
Jeremy Watts
"Madhusudan Singh" <spammers...@spam.invalid> wrote in message
> Hi
>
> I am solving a problem using f95 that involves solving a 12x12 coupled
> linear system. I have written out the matrix A (A x = f), and instead of
> calculating the inverse repeatedly for different values of the input
> parameters, I am trying to calculate the symbolic inverse of the matrix,
> post multiply it with the RHS and code the resulting 12 equations for x_i
> into fortran.
is there any particular reason why you want to calculate the symbolic
inverse? I should imagine its a mammoth task even for relatively small
matrices such as yours. As bourne out by your execution times. Also the
results would be likely prone to ill conditioning I should imagine.
By the way,a better way of solving linear systems is via the LU
decomposition or do it iteratively using Jacobi or Gauss-Seidel, rather than
finding the inverse.
Marcus Stollsteimer
> Out of a hunch, I just exported the notebook to text, and used
> "math" (the text version of Mathematica) to run it through. Less
> that two minutes later, I had an answer. The problem is that the
> answer (in nohup.out) is so ill formatted (I used MatrixForm) that
> it is impossible to make sense of.
>
> OTOH, the graphical Notebook is still chugging away. Has anyone seen
> this odd behaviour before ?
Hi,
generally: you should first try with smaller matrices to get
a feeling for the behavior.
in this case: try using ";" to suppress output
(see Mathematica book 1.4.8).
the problem is not the calculation of the inverse, but
rather the rendering of the output, which takes much more time:
here the time (sec) used for calculating and displaying the inverse
for matrices of various sizes
3: 0 0
4: 0.01 0.02
5: 0.02 0.18
6: 0.05 1.59
7: 0.11 14.55
8: 0.36 -
(for 8x8-matrix, rendering needs too much memory and crashes kernel.)
the times were obtained with
m = Table[a[i, j], {i, 3}, {j, 3}]
t0 = TimeUsed[];
inv = Inverse[m];
t1 = TimeUsed[];
inv
t2 = TimeUsed[];
t1 - t0
t2 - t1
Regards,
Marcus
--
Is there anything out there? And why?
And must they be so noisy? -- W. Allen
John Lansberry
news:425d988c$0$79462$1472...@news.sunsite.dk...
>
> I am solving a problem using f95 that involves solving a 12x12 coupled
> linear system. I have written out the matrix A (A x = f), and instead of
> calculating the inverse repeatedly for different values of the input
> parameters, I am trying to calculate the symbolic inverse of the matrix,
> post multiply it with the RHS and code the resulting 12 equations for x_i
> into fortran.
>
> I believe this is the sensible way to do it, since this system needs to be
> solved lots of times (its part of a half-infinite integral) and inverting
> such a system numerically would not only be time consuming, but also
> obscure possible stability problems.
It's hard to tell from your description whether the matrix A is always fixed
or if it will change from call to call. It's also hard to tell if your A
matrix has any special structure.
In general, I would quibble with your notion that the execution speed of
code that implements the symbolic inverse of a 12 x 12 matrix will be faster
than just solving the system of equations numerically. Even if Mathematica
could return the symbolic inverse of a general 12 x 12 matrix it would be so
complicated that the number of multiplications, divisions and additions
would be huge and very difficult to code in f95 (especially without making a
mistake). I would also quibble with your notion that computing the symbolic
inverse is better from a numerical stability point of view.
I would solve the system of equations using LAPACK in f95 (get it from
netlib.org). Assuming A is a general matrix (no special structure) and that
the A matrix does not change, then compute its LU factorization using XGETRF
(where X depends on the precision you want, single, double, etc.). Then,
for different values of F, use XGETRS to solve for x. If the A matrix does
change from call to call, then use XGESV.
I bet the speed of LAPACK will be far faster than anything you code up from
a symbolic inverse.
John
Albert Reiner
> I believe this is the sensible way to do it, since this system needs to be
> solved lots of times (its part of a half-infinite integral) and inverting
But if it is to be evaluated as some parameter varies along that
integration range, it seems likely that the coefficients of your
matrix and right hand side also vary sufficiently slowly. In that
case you might try to just keep the first order changes and check
whether that simplifies things sufficiently. Or you could try to get
the inverse of your matrix only at the beginning of the integration,
and to get an expression for its change along the integration path.
Albert.
Jimbo
Based on my past experience with Macsyma I believe that your run time
problem arises because you are resubstituting intermediate results into
down stream computations.
Herein lies the problem. You must compute your factors for the matrix
and write out the results for L and U, i.e., the full (sparse) symbolic
expressions. At this point these expressions will be numbers in your
real inverse matrix calculation. You can generate backsubstitution
symbolic equations, but only use single variable names for each
non-zero matrix location. (If you were to substitute your previous
symbolic solutions here you would generate an exponential explosion in
the number of terms-----remember in your Fortran code its just numbers
at this point of the calculation). This same strategy applies for both
the L and U backsubstitutions.
For the record, I tried a Masyma LU symbolic demo case for inverting a
4x4 matrix. The result was humongous!----400k operations (I checked
with utilities). This was obviously nonsense. I got into the code and
preformed the rearrangement described above and the operation count was
reduced to ~4k (a 1000-to-1 reduction!!!!!). The numerical answers
were the same. As for myself, I preferred the second approach!!! In
fact, applying this strategy I have been able to symbolically invert
50x50 systems (file size~2MB).
Good luck in you endeavours, trust the numbers, and for Heavens sake
don't resubstitute intermediate symbolic results.
Jim Turner
Manager, Business Development
Dynacs, Military and Defense
Houston, TX 77058
281-226-5213
jtu...@dynacs.com
also
Adjunct Associate Professor, TAMU AERO
Albert Reiner
> don't resubstitute intermediate symbolic results.
In Maxima, there is an optimize command that goes part of the way to
doing that:
,----
| - Function: OPTIMIZE (exp)
| returns an expression that produces the same value and side
| effects as exp but does so more efficiently by avoiding the
| recomputation of common subexpressions. OPTIMIZE also has the side
| effect of "collapsing" its argument so that all common
| subexpressions are shared. Do EXAMPLE(OPTIMIZE); for examples.
`----
But it doesn't go all the way:
,----
| (%i2) optimize(1+a+b+sin(a+b)+cos(a+b));
|
| (%o2) BLOCK([%1], %1 : b + a, SIN(%1) + COS(%1) + b + a + 1)
`----
Albert.
Richard Fateman
.. <snip>...
> Even if Mathematica
> could return the symbolic inverse of a general 12 x 12 matrix it would be so
> complicated that the number of multiplications, divisions and additions
> would be huge and very difficult to code in f95 (especially without making a
> mistake).
Most CAS have facilities to dump formulas out into Fortran, C, or other
forms. This should decrease the probability of copying errors.
> I would also quibble with your notion that computing the symbolic
> inverse is better from a numerical stability point of view.
I agree, though perhaps SEEING the 12x12 sparse matrix would make
this discussion more concrete. It could be that there is a nice
solution, and that one of the CAS would find it. There are a number
of papers over the last 30 years which discuss various techniques including
graph-theoretic decompositions, hash-coding of common expressions,
etc.
Someone also asked about Macsyma/Maxima. It is
plausible to do a symbolic matrix inverse in this program. How much
time is taken depends on which of several possible algorithms is
chosen.
RJF
Madhusudan Singh
>
> "Madhusudan Singh" <spammers...@spam.invalid> wrote in message
> news:425d988c$0$79462$1472...@news.sunsite.dk...
>> Hi
>>
>> I am solving a problem using f95 that involves solving a 12x12 coupled
>> linear system. I have written out the matrix A (A x = f), and instead of
>> calculating the inverse repeatedly for different values of the input
>> parameters, I am trying to calculate the symbolic inverse of the matrix,
>> post multiply it with the RHS and code the resulting 12 equations for x_i
>> into fortran.
>
> is there any particular reason why you want to calculate the symbolic
> inverse? I should imagine its a mammoth task even for relatively small
> matrices such as yours. As bourne out by your execution times. Also the
> results would be likely prone to ill conditioning I should imagine.
>
The long execution time seems to be related to the time Mathematica spends
on rendering the result as Marcus pointed out.
Let me describe what I did next :
"
Out of a hunch, I just exported the notebook to text, and used "math" (the
text version of Mathematica) to run it through. Less that two minutes
later, I had an answer. The problem is that the answer (in nohup.out) is so
ill formatted (I used MatrixForm) that it is impossible to make sense of.
"
So, if math (text version of Mathematica) can churn out the answer in 1.5
minutes, the problem cannot lie in calculating the inverse but in rendering
the result. Do you know if there is a way to make Mathematica render it any
faster ?
Albert Reiner
> So, if math (text version of Mathematica) can churn out the answer in 1.5
> minutes, the problem cannot lie in calculating the inverse but in rendering
> the result. Do you know if there is a way to make Mathematica render it any
> faster ?
You might try InputForm[] instead of MatrixForm[].
Albert.
Madhusudan Singh
> Hi,
>
> generally: you should first try with smaller matrices to get
> a feeling for the behavior.
>
> in this case: try using ";" to suppress output
> (see Mathematica book 1.4.8).
>
> the problem is not the calculation of the inverse, but
> rather the rendering of the output, which takes much more time:
>
Aha !
The only place where I am o/p anything is the following :
Print[MatrixForm[Ainv]];
Print[MatrixForm[Ainv.f]];
Let me try commenting out the first one and get back to you.
This was an illuminating response. Thank you.
car...@colorado.edu
element research since 1982. The biggest symbolic matrix I have
inverted was order 48 but that had special structure. It took about
30 hours. If the matrix entries are general and complicated
expressions it is usually hopeless to go beyond, say, 8 or 10.
In the worst case the number of terms increases factorially. For
example the adjoint of a 48 x 48 will have 48! =
124139155925360726708622890473 733750385214863546777600000000 00
terms, or O(10^61), which is greater than the number of atoms in the
universe.
Surprisingly for symbolic matrices the deprecated Cramer's
rule often works much better than LU factorization. The
reason is that taking determinants involves only products. No
intermediate divisions and no pivot searches. The first thing I do
when confronted with a symbolic A to invert is to take
Det[A] and print it. Then I try to simplify it. If it
doesnt look too bad and takes reasonable time, say minutes, I run my
own inverse program and hope it will finish within a few hours:
MyInverse[A_]:=Module[{i,j,B,C ,detA,detC,n=Length[A],z},
detA=Simplify[Det[A]]; z=Table[0,{n}]; B=Table[0,{n},{n}];
For [i=1,i<=n,i++, For[j=1,j<=n,j++, C=Transpose[A];
C[[j]]=z; C[[j,i]]=1; C=Transpose[C];
detC=Det[C]; B[[j,i]]=Simplify[detC]/detA ]];
Return[B]];
To test:
Ainv=MyInverse[{{a,b},{c,d}}]; Print[Ainv//MatrixForm];
Peter Spellucci
In article <425d988c$0$79462$1472...@news.sunsite.dk>,
a 12 by 12 inverse has 144 elements and these will in most cases be nonzero
even if the original matrix is sparse (for example the inverse of an
irreducible tridiagonal matrix is completely dense. )
hence you force Mathematica to compute these 144 formulas involving
coefficients of A which itself are functions of several variables as you
write. but the formal inverse elements are quotients of determinants,
all 11 by 11 minors in the numerator (the elements of the adjoint matrix)
and the determinant itself in the denominator. but a n by n determinant
has factorial(n) summands each on consisting of n factors.
factorial(12).=.10^8 hence you cannot hope to get a result even with a much
larger and faster machine
hth
peter
car...@colorado.edu
the fastest method for *symbolic* matrices. The main problem is
that each step introduces one level of symbolic division.
Such "nested division" expressions can be slow to
tunravel for simplification at the end.
Here is a timing example. The 8 x 8 matrix A below is a typical
connection matrix of a two-dimensional stress-hybrid finite element:
A={{(-2*a*C11-b*C13)/4,(-2*a*C12-b*C23)/4,(-2*a*C13-b*C33)/4,(a*C11)/4,0,a,0,b/2},
{(-2*b*C12-a*C13)/4,(-2*b*C22-a*C23)/4,(-2*b*C23-a*C33)/4,0,(b*C22)/4,0,b,-a/2},
{(2*a*C11-b*C13)/4,(2*a*C12-b*C23)/4,(2*a*C13-b*C33)/4,-(a*C11)/4,0,a,0,b/2},
{(-2*b*C12+a*C13)/4,(-2*b*C22+a*C23)/4,(-2*b*C23+a*C33)/4,0,-(b*C22)/4,0,b,
a/2},{(2*a*C11+b*C13)/4,(2*a*C12+b*C23)/4,(2*a*C13+b*C33)/4,(a*C11)/4,0,a,0,-b/2},
{(2*b*C12+a*C13)/4,(2*b*C22+a*C23)/4,(2*b*C23+a*C33)/4,0,(b*C22)/4,0,b,a/2},
{(-2*a*C11+b*C13)/4,(-2*a*C12+b*C23)/4,(-2*a*C13+b*C33)/4,-(a*C11)/4,0,a,0,-b/2},
{(2*b*C12-a*C13)/4,(2*b*C22-a*C23)/4,(2*b*C23-a*C33)/4,0,-(b*C22)/4,0,b,-a/2}};
Using LUDecomposition, followed by solving for 8 unit RHS and
simplifying, takes 22 sec to get Ainv on my laptop.
With the brute-force Cramer's rule (which BTW I teach undergraduate
students NOT to use for computations):
MyInverse[A_]:=Module[{i,j,B,C,detA,detC,n=Length[A],z},
detA=Simplify[Det[A]]; z=Table[0,{n}]; B=Table[0,{n},{n}];
For [i=1,i<=n,i++, For [j=1,j<=n,j++, C=Transpose[A];
C[[j]]=z; C[[j,i]]=1; C=Transpose[C]; detC=Simplify[Det[C]];
B[[j,i]]=detC/detA; ]];
Return[B]];
Print[Timing[Simplify[MyInverse[A]]]];
produces Ainv in 7 sec.
Leo Harten
of 10^80.
<car...@colorado.edu> wrote in message
news:1113506921.6...@o13g2000cwo.googlegroups.com...
>I had similar experiences since I use CAS to do finite
> element research since 1982. The biggest symbolic matrix I have
> inverted was order 48 but that had special structure. That took about
> 30 hours. If the matrix entries are general and complicated
> expressions it is usually hopeless to go beyond, say, 8 or 10.
>
> In the worst case the number of terms increases factorially. For
> example the adjoint of a 48 x 48 will have 48! =
> 12413915592536072670862289047373375038521486354677760000000000
> terms, or O(10^61), which is greater than the number of atoms in the
> universe.
>
> Surprisingly for symbolic matrices the deprecated Cramer's
> rule often works much better than LU factorization. The reason is that
> taking determinants involves only products. No intermediate
> divisions and no pivot searches. Also it takes no extra memory. The
> first thing I do when confronted with a symbolic A to invert is to take
> Det[A] and print it. Then I try to simplify it. If it doesnt look too
> bad and takes reasonable time, say minutes, I run my own inverse
> program and hope it will finish within a few hours:
>
> MyInverse[A_]:=Module[{i,j,B,C,detA,detC,n=Length[A],z},
> detA=Simplify[Det[A]]; z=Table[0,{n}]; B=Table[0,{n},{n}];
> For [i=1,i<=n,i++, For[j=1,j<=n,j++, C=Transpose[A];
> C[[j]]=z; C[[j,i]]=1; C=Transpose[C];
Leo Harten
held at General Electric, Schenectady, NY, Ed V. Ellen Lewis and M.A.
Hussain,
article is "Using MACSYMA to Generate (Somewhat) Optimized FORTRAN Code"
and includes the source code, and examples with analysis. A focus was on the
computation of a matrix with symbolic entries and its derivative.
"Albert Reiner" <are...@tph.tuwien.ac.at> wrote in message
news:vw8d5sx...@berry.phys.ntnu.no...
Jimbo
Optimize is not the issue. You only use optimize before you print out
the intermediate results. I wrote my own LU macsyma code so that I
could generate the solution in the most efficient way.
Jim Turner
Manager, Business Development
Dynacs, Military and Defense
Houston, TX 77058
281-226-5213
car...@colorado.edu
hope to invert a full 48x48.
Dave Rusin
Madhusudan Singh <spammers...@spam.invalid> wrote:
>I am solving a problem using f95 that involves solving a 12x12 coupled
>linear system. I have written out the matrix A (A x = f), and instead of
>calculating the inverse repeatedly for different values of the input
>parameters, I am trying to calculate the symbolic inverse of the matrix,
It might be of interest to see an example. I don't think I saw in this
thread a specification of what the matrix looked like -- how many
"parameters"? how do they appear? (i.e. is A_ij a linear combination
of them? polynomial? rational function?)
Seems to me I heard a tale of such a symbolic matrix inversion problem
which was attacked head-on; only when the (remarkable) answer was finally
known did someone stop to ask the origin of the matrix and was told that
A was known from the outset to be orthogonal...
dave
spasmous
Madhusudan Singh wrote:
>... Awaiting further suggestions.
A simple test using a diagonal matrix in MATLAB:
syms A a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12
A = diag([a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12]);
inv(A)
ans = ... (the right answer)
The time taken was a fraction of a second. I checked up to 5x5 for a
full matrix and inv() still produced answers within seconds -
presumably the correct answers, I didn't check. The equations aren't
pretty although having lots of zeros in your matrix will help.
Unless your matrix is singular:
A=[a1 0;a2 0];
inv(A)
??? Error using ==> sym.inv
Error, (in inverse) singular matrix
As long as it's not singular, I guess a 12x12 symbolic inverse will
work fine in MATLAB.
spasmous
for a 10x10 symbolic inverse.
syms A a b c d e f g h i j k l m n o p q r s t u v w x y z;
syms aa bb cc dd ee ff gg hh ii jj kk ll mm nn oo pp qq rr ss tt uu vv;
A = [a 0 b 0 c 0 d 0 e 0 ;
f 0 0 g h 0 0 i j pp;
0 k l 0 0 m n o 0 0 ;
p 0 0 q r s 0 0 t 0 ;
0 u 0 v w 0 x y 0 qq;
aa 0 0 bb 0 0 cc 0 0 rr;
0 dd ee 0 0 0 ff 0 gg 0 ;
hh ii 0 0 0 0 0 0 jj 0 ;
0 0 0 0 0 kk ll mm nn ss;
tt 0 0 0 uu 0 0 0 vv 0];
B=inv(A);
??? Error using ==> sym.inv
Error, (in normal/nlinear) integer too large in context
Madhusudan Singh
Thanks for the two posts though. I could have wasted some time trying out
matlab for this.
Mathematica is happily handling the problem (the raw inverse of the 12x12
matrix takes about 11 seconds), except that a subsequent FullSimplify step
to clean up the expressions is causing some memory problems on some
machines.
Matlab/octave are good for numerical problems involving matrices (though I
would always code anything that needed a lot of speed in f95) but as far as
symbolic stuff goes, Mathematica is very nearly the king (my experience of
maxima is limited, and of maple, almost non-existent).
car...@colorado.edu
and labs.
Bob Lewis and Susan Lewis
>
> > Actually I retract what I said - MATLAB bails out
> after a few minutes
> > for a 10x10 symbolic inverse.
>
> Thanks for the two posts though. I could have wasted
> some time trying out
> matlab for this.
>
> Mathematica is happily handling the problem (the raw
> inverse of the 12x12
> matrix takes about 11 seconds), except that a
> subsequent FullSimplify step
> to clean up the expressions is causing some memory
> problems on some
> machines.
>
> Matlab/octave are good for numerical problems
> involving matrices (though I
> would always code anything that needed a lot of speed
> in f95) but as far as
> symbolic stuff goes, Mathematica is very nearly the
> king (my experience of
> maxima is limited, and of maple, almost
> non-existent).
I just noticed this thread. Sorry to be jumping in so late.
(1) Can we see the actual matrix you want to invert? In my experience the number of symbolics and their placement will make a huge difference. If it is a secret, can you post an "equivalent" or slightly reduced case to test? I get the impression that the matrix contains only polynomials. Is that right?
(2) In the experience of me and quite a few other people Mathematica does poorly at this kind of problem. The two best systems for polynomials and matrices are Fermat (written by me) and Magma.
(3) You say Mathematica "finished" in 11 seconds but wouldn't or couldn't simplify. That is a common behavior with several large systems; they really are no way near completion -- at least as most people would define completion. When Michael Wester and I tested CA systems in 1999 [ref below] we met this behavior often.
(4) I took an example from my web site [below] and reduced it to 12x12 to invert. The matrix has 12 different symbolic vars and about 98 entries that are 0 or 1. It takes Fermat about 12 minutes to invert, and about 400 meg of RAM. The entries in the inverse have around 10000 terms.
(5) Someone mentioned LU decomposition. I used Fermat to produce the fraction-free LU decomposition [ref below] of the above matrix. The L part has several entries with more than 1000000 terms.
Hope this helps.
Robert H. Lewis
Fordham University
R. Lewis and M. Wester, Comparison of Polynomial-oriented Computer Algebra Systems, SIGSAM Bulletin, Jan 2000.
Nakos, Turner, and Williams, Fraction-free Algorithms for Linear and Polynomial Equations, SIGSAM Bulletin Sept 1997.
Corless and Jeffrey, The Turing Factorization of a Rectangular Matrix, SIGSAM Bulletin Sept 1997.
Madhusudan Singh
Thanks for your post.
I managed to work out the memory problems and finish the simplification in
two steps (memory ran out halfway through simplification).
Finally, taking the inverse took 11 seconds. Taking the dot product with RHS
took 1 second, and simplification of the resulting 12 components took about
2 hours in all.