Ask Your Question

Solve a "big" polynomial system numerically

asked 2013-10-16 00:25:28 -0500

I would like to solve numerically a polynomial system with around 10^5 variables and 10^7 equations. Each equation is of degree 4 with integer coefficients and contains around 10 variables, so that the jacobian of the system is a very "hollow" matrix with integer entries.

Let X=[x_0,x_1, x_2,...] and E=[eq_0,eq_1, eq_2...] the list of variables and equations :

Is there a SAGE function f(E,X) which find numerically solutions for such a system?

Is SAGE relevant for such intensive computation ?

I think the code needs to be 100% cython, for not losing time.

My computer has a hard disk of 980 Go, 4Go of RAM and its processor "Intel Core i5 CPU 760 @ 2.80GHz × 4" is around 10^9 FLOPS.

I guess my computer is not enough powerful for doing such computation in a reasonable time :
The jacobian J is a (10^5)x(10^7) matrix, so for example, the Gauss-Newton algorithm would require around 1000 Go of RAM and 10^16 operations.

I know very few things in numerical analysis, so some advices would be welcome.

edit retag flag offensive close merge delete

2 answers

Sort by » oldest newest most voted

answered 2013-10-16 07:04:24 -0500

Shashank gravatar image

If you think the problem is going to take a long time, then using cython is recommended. Start the notebook with %cython in the beginning and declare the matrix as a numpy matrix instead of a sage matrix, using command cdef np.matrix. Also the scipy algorithm for leastsq is very similar to newton-gauss method. You can use scipy.optimise.leastsq to solve the equation.

I hope it helps. I normally work with small matrices so I might be wrong, but I have to loop over a lot of matrices and using scipy within cython is much faster than default sage algorithms.

edit flag offensive delete link more


Thank you very much, I will try scipy algorithm and numpy matrix. Is it compatible with the "sparse mode" for matrix ?

Sébastien Palcoux gravatar imageSébastien Palcoux ( 2013-10-16 22:13:06 -0500 )edit

If you want sparse matrix you can use the scipy matrix instead to numpy.

Shashank gravatar imageShashank ( 2013-10-17 06:54:13 -0500 )edit

answered 2013-10-16 05:25:10 -0500

tmonteil gravatar image

If it can help, the mathematical english for "creuse" is "sparse", not "hollow". So you can have a look at:

sage: M = matrix(ZZ, 10^5, 10^7, sparse=True)
sage: M
100000 x 10000000 sparse matrix over Integer Ring (type 'print M.str()' to see all of the entries)
sage: M[3,4] = 10
sage: M[10,10]
sage: M[3,4]

And then see if Sage can survive your computations.

edit flag offensive delete link more


Thank you very much ! (I also do my PhD at the IML, Luminy, with Antony Wassermann in operators algebras, 2005-2009). So if I understand well, in mode sparse, the computer memorizes only the non-zero entries, and so in my case, this needs only 100 Mo of RAM (instead of 1000 Go).

Sébastien Palcoux gravatar imageSébastien Palcoux ( 2013-10-16 22:11:53 -0500 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower


Asked: 2013-10-16 00:25:28 -0500

Seen: 347 times

Last updated: Oct 16 '13