roman_pearce

Mr. Roman Pearce

1678 Reputation

19 Badges

20 years, 216 days
CECM/SFU
Research Associate
Abbotsford, British Columbia, Canada

I am a research associate at Simon Fraser University and a member of the Computer Algebra Group at the CECM.

MaplePrimes Activity


These are replies submitted by roman_pearce

The upshot is that most of the research done in Maplesoft-sponsored labs has become less about core Maple and more about specialized mathematics (there is great work being done, just not about core). I disagree ;)

I thought I would post a few more times for linear systems with smaller entries. The systems above have large solutions, ie: the numerators and denominators of the last solution have 20000 bits, or about 6000 decimal digits each. That's a lot of lifting steps.

These times are all from a Pentium 4 3.2 GHz, computer running 32-bit Maple 12. The matrix entries range from -5 to 5. The numerators and denominators of the last system have about 4000 decimal digits. On structured systems (where the solution is expected to be small) you could obviously get much better performance.

read "dpadic.mpl":
infolevel[solve] := 0:
with(LinearAlgebra):
for N from 8 to 11 do
  A := RandomMatrix(2^N, 2^N+1,generator=-5..5):
  TIMER := time(DensePadicLift(A,2^N));
  print(2^N = TIMER);
end do:
                    256 = 1.276

                    512 = 7.772

                   1024 = 54.499

                   2048 = 496.655

The dimension of the system really is 2, and if you substitute y2=1 you get a system of dimension 1.  I computed a lex Groebner basis and solve returned no solutions.  This ticked me off, so started manually eliminating linear equations.

sys := {x1^2-x2^2-x3+x4, (x2+1)^2-x5, 125+x1 - x2, 
1+x3 - x4, x8*x6+x7 - x8^2, x9*x6+x7 - x9^2-x3, 
y2*x6+x7 - x5-y2^2, y1*x6+x7+x5 - x5-y1^2-x4, 
x4+(x1-x2-1)*(x1+x2+1)}; 
tord := plex(x5, x3, x7, x6, x4, y2, y1, x9, x8, x1, x2):
R := Groebner[Basis](sys, tord);              
S, R := selectremove(a->degree(a)=1, R):
S := solve(S);
R := remove(`=`, subs(S, R), 0); 

This knocks out 5 variables and leaves us with an irreducible system of dimension 2 (prime ideal).  Solve does not return any solutions for this system.  It appears that the code that constructs solutions is buggy.

This is the bug.  In general, we know the dimension is two, but there exists a "maximal independent set" of variables V which, if we specialize them (put them in the coefficient ring) the system becomes zero dimensional (finite number of solutions in remaining variables).  If you pick any old set of variables, you might not get any solutions, or you might get something with the wrong dimension, etc.  It would appear that solve is not being careful about this.  The Groebner command to compute this set is Groebner:-MaximalIndependentSet.  It uses a Groebner basis, so you should call it with the Groebner bases already computed in solve.

In SolveTools:-UseGroebner:

elif dim <> 0 then ...
  gsol := Groebner:-Solve(eqs, 'plex'(orderG), notzero);
  sols := NULL;
  for i in gsol while nops([sols]) < maxsols do
     sols := sols, SolveTools:-PolynomialSystem(
     op(1, i), vars, op(3, i), maxsols - nops([sols]), 'avoidgroebner')
  end do

Change this to something like:
elif dim <> 0 then ...
  gsol := Groebner:-Solve(eqs, 'plex'(orderG), notzero);
  sols := NULL;
  for i in gsol while nops([sols]) < maxsols do
     vars1 := vars minus Groebner:-MaximalIndependentSet(op(1,i),op(2,i));
     sols := sols, SolveTools:-PolynomialSystem(
       op(1, i), vars1, op(3, i), maxsols - nops([sols]), 'avoidgroebner')
  end do

Note that in cases where solve does not use Groebner (ie: the old resultant code is faster) you may be hooped.

The dimension of the system really is 2, and if you substitute y2=1 you get a system of dimension 1.  I computed a lex Groebner basis and solve returned no solutions.  This ticked me off, so started manually eliminating linear equations.

sys := {x1^2-x2^2-x3+x4, (x2+1)^2-x5, 125+x1 - x2, 
1+x3 - x4, x8*x6+x7 - x8^2, x9*x6+x7 - x9^2-x3, 
y2*x6+x7 - x5-y2^2, y1*x6+x7+x5 - x5-y1^2-x4, 
x4+(x1-x2-1)*(x1+x2+1)}; 
tord := plex(x5, x3, x7, x6, x4, y2, y1, x9, x8, x1, x2):
R := Groebner[Basis](sys, tord);              
S, R := selectremove(a->degree(a)=1, R):
S := solve(S);
R := remove(`=`, subs(S, R), 0); 

This knocks out 5 variables and leaves us with an irreducible system of dimension 2 (prime ideal).  Solve does not return any solutions for this system.  It appears that the code that constructs solutions is buggy.

This is the bug.  In general, we know the dimension is two, but there exists a "maximal independent set" of variables V which, if we specialize them (put them in the coefficient ring) the system becomes zero dimensional (finite number of solutions in remaining variables).  If you pick any old set of variables, you might not get any solutions, or you might get something with the wrong dimension, etc.  It would appear that solve is not being careful about this.  The Groebner command to compute this set is Groebner:-MaximalIndependentSet.  It uses a Groebner basis, so you should call it with the Groebner bases already computed in solve.

In SolveTools:-UseGroebner:

elif dim <> 0 then ...
  gsol := Groebner:-Solve(eqs, 'plex'(orderG), notzero);
  sols := NULL;
  for i in gsol while nops([sols]) < maxsols do
     sols := sols, SolveTools:-PolynomialSystem(
     op(1, i), vars, op(3, i), maxsols - nops([sols]), 'avoidgroebner')
  end do

Change this to something like:
elif dim <> 0 then ...
  gsol := Groebner:-Solve(eqs, 'plex'(orderG), notzero);
  sols := NULL;
  for i in gsol while nops([sols]) < maxsols do
     vars1 := vars minus Groebner:-MaximalIndependentSet(op(1,i),op(2,i));
     sols := sols, SolveTools:-PolynomialSystem(
       op(1, i), vars1, op(3, i), maxsols - nops([sols]), 'avoidgroebner')
  end do

Note that in cases where solve does not use Groebner (ie: the old resultant code is faster) you may be hooped.

Maple can also be very sensitive to hardware configurations [like cache size for numerics intensive computations]

Optimized, compiled routines can tell you quite a bit about the efficiency of the processor and memory architecture, however hardware review sites already report on this.  Synthetic benchmarks are a pretty good indicator of theoretical peak performance.  I'm just not sure what the performance of the Maple kernel (never mind the library) would be measuring.  Maplesoft could make a small tweak and radically change the timings.  It would be useful to have a measure of Maple performance over time, because some things have definitely slowed down.  Polynomial multiplication, for example, was much faster in Maple VR3-VR5, than Maple 7: http://www.cs.berkeley.edu/~fateman/papers/fastmult.pdf

And it can vary on whether you are running the computation using TTY, Classic or Standard

I believe this is because parts of the standard interface (like 2d math typesetting) is written in the Maple language.  This gets interpreted by the kernel, and the extra memory adds to the garbage collection time.

I have been thinking about creating a Maple benchmark, to show a performance ranking versus operating system and hardware.

Is there any interest in such a thing?

Yes, however Maple's performance tends to be fairly unaffected by a lot of things. I think it would be better to have benchmarks covering specific functionality: solving a linear system, factoring polynomials, etc. Then you could subtract algorithm time to measure things like the overhead of solve or simplify. I would also like to know the cost to evaluate functions at the default setting of Digits, and the time for the GUI to print X megabytes of 2d math. Designing benchmarks like this is quite a lot of work.

I'd hesitate to do this because it think it would be too hard to maintain. The Maple library often changes in a lot of small ways, and patches are very likely to introduce strange behavior in top level commands like solve, simplify, etc. It would be better, I think, to build a repository of good algorithms. Write standalone routines that perform a well-defined task and improve on Maple's. Maple's language can be very clear and concise, so if the code were documented it would be short and you could port to other languages easily.
On a 32-bit machine you can exhaust the 4GB address space. That's 32 billion bits, or log[10](2^32) decimal digits. In practice most operating systems will only let you use 2GB of RAM, so 10^9 digits is more realistic. Just don't expect to do anything with that number. For 64-bit machines, GMP still uses 32-bit (signed) integers for the number of words allocated, but words are now 64-bits. The largest integer you can represent is 16GB (2^31 words) and that gives you 10^11 decimal digits.
4800 seconds is an hour and twenty minutes! That is some obscene overhead for what boils down to a simple type check + software float. Then I looked at `evalf/erf`. What is the overhead of try/catch ? It seems to catch a numeric exception for the sole purpose of returning the same error. Here is something I wish Maplesoft would keep in mind (and I'm sure you agree): a factor of 300 slowdown makes otherwise useful software useless. And it shows. It shows badly whenever someone tries to do anything "modest", let alone big. The fact is, it took 512000 CPU cycles to evaluate each call to erf. Now Maple is interpreted and there are some underlying issues, but that number is just not defensible. I think in the long run, it would make sense to know the cost of software float arithmetic and functions. We should know, approximately, the number of cycles per bit of precision. And it better not be more than 10000. Nobody is going to take this stuff seriously otherwise. Serious projects do this. They estimate the cost. Then you get predictable performance. It also helps you argue that the software is fast. It doesn't matter that Maple is interpreted. We should know the cost of procedure calls and all that stuff in cycles (approximately). (edit: I calculated the CPU cycles wrong)
Try using the display3d command instead of display (I agree, this is silly).
The error message is telling you what to do - put the plots in a list or a set. display([CPA, CPB, CPC], title="...");
The problem certainly does exist and I think it's in Java/PowerPC. I use command line Maple a lot and it has no problems. The issue may be fixed in 10.4, I don't know. 10.3 is getting pretty old so I'll install 10.5 soon (iBook, 1GHz) and report back on this. As for new releases, Maplesoft tends to release every 12-18 months.
I don't agree. Providing a function for something this simple is too much, and I think it would be better to give an example of the add command on the help page for CoefficientList/Vector. I can see having a command if it could do more things, some of them hard. Horner form is good. Maybe it could efficiently convert a list of coefficients to a polynomial in different bases, the usual power basis being one of them.
I don't agree. Providing a function for something this simple is too much, and I think it would be better to give an example of the add command on the help page for CoefficientList/Vector. I can see having a command if it could do more things, some of them hard. Horner form is good. Maybe it could efficiently convert a list of coefficients to a polynomial in different bases, the usual power basis being one of them.
The internals of Maple can still be improved. I think the best way to scale the system is to develop special purpose high performance libraries and call them for large problems. It also helps to add "programmer packages" such as LinearAlgebra:-Modular that allow Maple developers to access high speed compiled code. With packages like that you can (carefully) code a high-level algorithm in Maple and get 90% of the speed and efficiency of C.
First 18 19 20 21 22 23 24 Last Page 20 of 39