# Solve a "big" polynomial system numerically

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 close merge delete

Sort by » oldest newest most voted

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]
0
sage: M[3,4]
10


And then see if Sage can survive your computations.

more

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.

more

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

( 2013-10-17 05:13:06 +0200 )edit
1

If you want sparse matrix you can use the scipy matrix instead to numpy. http://docs.scipy.org/doc/scipy/reference/sparse.html

( 2013-10-17 13:54:13 +0200 )edit