9
$\begingroup$

This is related to Noah's recent question about solving quadratics in a number field, but about an even earlier and easier step.

Suppose I have a huge system of linear equations, say ~10^6 equations in ~10^4 variables, and I have some external knowledge that suggests there's a small solution space, ~100 dimensional. Moreover, the equations are sparse; in fact, the way I produce the equations gives me an upper bound on the number of variables appearing in each equation, ~10. (These numbers all come form the latest instance of our problem, but we expect to want to try even bigger things later.) Finally, all the coefficients are in some number field.

Which computer algebra system should I be using to solve such a system? Everyone knows their favourite CAS, but it's often hard to get useful comparisons. One significant difficulty here is that even writing down all the equations occupies a big fraction of a typical computer's available RAM.

I'll admit that so far I've only tried Mathematica; it's great for most of our purposes, but I'm well aware of its shortcomings, hence this question. A previous slightly smaller instance of our problem was within Mathematica's range, but now I'm having trouble.

(For background, this problem is simply finding the "low weight spaces" in a graph planar algebra. See for example Emily Peter's thesis for an explanation, or our follow-up paper, with Noah Snyder and Stephen Bigelow.)

Improve this question
$\endgroup$
1
  • $\begingroup$ If you feel up to it, a minimal description of the syntax required to import a large system of sparse linear equations into your preferred CAS would be great! $\endgroup$
    – Kim Morrison
    Commented Oct 25, 2009 at 20:16

8 Answers 8

Reset to default
3
$\begingroup$

It probably goes without saying that solving linear systems over number fields is probably far from the being among the most important user-level functionality of the main commercial computer algebra systems. That said, I do know that this functionality in Maple was written by someone with a specific interest in this sort of thing. If you have access to a recent version of Maple, take a look at the help page: ?SolveTools,Linear. 10^6 is pretty big, but it might still be within reach of the solver.

While, I do not know much about how Mathematica does these things, I do know that in Maple sparse linear systems are more efficiently solved as polynomials (rather than sparse-matrices) since underlying polynomial data-structure turns out to be well suited to sparse system solving.

If Maple does not work for you (or you do not have access to it), this strikes me as exactly the sort of problem that MAGMA might be targeting.

Improve this answer
$\endgroup$
4
$\begingroup$

LinBox

May be quite helpful, as it contains routines for finite fields and has an interface to Maple (which is significantly easier to use than LinBox as a C++ library).

Hope that helps!

Improve this answer
$\endgroup$
4
$\begingroup$

If Magma can do this, you may well look at Sage, which is open source, remarkably powerful, and with support for sparse linear algebra.

Improve this answer
$\endgroup$
3
$\begingroup$

This sounds like a problem you'll really need to run through a strong solver in C, and if you want it to pull in objects easily it'd probably be best to use a Perl C-wrapper with MPI and Parallel LAPACK. I'll look into easily accessible sparse matrix solvers, that are modular enough to change the inputs (oy that's gonna be interesting, but I'll take a look). Honestly at numbers that high I don't see any reasonably sized computer that matlab or mathematica can deal with solving it.

Will be back with results of search soon.

Improve this answer
$\endgroup$
3
$\begingroup$

For a problem of this size, you should consider which algorithm before you consider which computer algebra system.

In floating point arithmetic, there is an excellent algorithm for solving a sparse system of linear equations. Even for dense matrices it is more robust than Gaussian elimination; for sparse matrices it is better still. This is the conjugate gradient method. It is available in certain packages and certain computer algebra systems, but it is also not very hard to implement from scratch. For instance, to me it feels simpler to code conjugate gradient than to code Gaussian elimination.

If you find floating point solutions to high precision, there are algorithms to convert them to elements of a number field with bounded complexity. (See the Inverse Symbolic Calculator, etc.) Some of these also exist in computer algebra systems. If you use one of these solvers, you can then check that it is an exact solution to the original system of equations.

If the solution space is low-dimensional but not 1-dimensional, then you can run conjugate gradient repeatedly after adding constraints to eliminate the kernel found so far. Conjugate gradient is an algorithm to minimize a positive-definite quadratic form. (If the system is not positive definite symmetric, you have to square or make some other change to fix that.) Thus it still works if you add constraints.

Improve this answer
$\endgroup$
2
$\begingroup$

Matlab might be a better choice. Be sure to use the sparse matrix functionality rather than just regular matrices. I haven't used Matlab for this purpose but I've seen people using it for large systems.

Also, have you been using sparse matrices in Mathematica?

http://reference.wolfram.com/mathematica/howto/WorkWithSparseMatrices.html

This might improve the performance.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Embarrassingly, it occurs to me that I might not have tried using Mathematica's sparse matrices. I've actually been manipulating the equations themselves, using an ad-hoc mixture of Mathematica's built-in "Solve" command and a few hand-written tricks. $\endgroup$
    – Kim Morrison
    Commented Oct 26, 2009 at 3:51
  • $\begingroup$ How does Matlab cope with coefficients in number fields? My minimal experience with Matlab didn't get that far. $\endgroup$
    – Kim Morrison
    Commented Oct 26, 2009 at 4:08
  • $\begingroup$ My Matlab knowledge is not huge but my best guess would be that being aware of number fields is not a Matlab sort of functionality as I would think it requires some sense of symbolic manipulation. Noah's question seemed to suggest that involving number fields was an optional idea directed toward trying to speed things up. However, I think porting from Mathematica to Matlab might entail some speed up. The black box nature of Mathematica makes it harder (at least for me) to do optimizations in Mathematica than in Matlab. $\endgroup$
    – Kim Greene
    Commented Oct 26, 2009 at 4:24
  • $\begingroup$ It's absolutely essential that we get exact answers. If matlab can't cope with number fields, it's useless to us. To clarify, in Noah's question it was also essential to end up with exact answers (although I admit that our current technique is to solve things numerically then approximate by algebraic integers). Noah's point was that we have external knowledge that the solution still lies in the same number field, even though we're solving quadratics and the field isn't algebraically closed! He was wondering if it's possible to translate this knowledge into some sort of speedup. $\endgroup$
    – Kim Morrison
    Commented Oct 26, 2009 at 16:08
  • $\begingroup$ One idea is that it's easier to write tighter code in Matlab and therefore, you could get a faster but still numerical answer. The second idea is to optimize the Mathematica code. $\endgroup$
    – Kim Greene
    Commented Oct 27, 2009 at 3:39
2
$\begingroup$

A system that large is large enough where the speedup from skipping over all the layers of abstraction in most CAS packages is worth the trouble to write your own custom code to solve the problem.

If you are ok with a floating point solution, ScaLAPACK or another linear algebra package that has algorithms and data structures for sparse matrices would be a lot better than (P)LAPACK, which as far as I know uses only dense matrix data structures.

Python's scipy package also has sparse data structures and wrappers to call UMFPACK, and the syntax is easy enough that it wouldn't be significantly harder to use than a CAS program. It would be easier than writing a custom Fortran/C program straight up.

For solutions over the reals, there is an algorithm called the complete orthogonal decomposition (COD) that uses rank revealing QR factorization (see, e.g. Golub's Matrix Computations). This lets you separate out and project away the kernel of the problem, leaving behind only the part of the problem that lies in the range. I don't know if there is an analogue for arbitrary fields, but since you appear to have a problem with a small rank, it may be worth your trouble to look into this.

Improve this answer
$\endgroup$
0
$\begingroup$

I think here is someone who does related research:

http://www4.ncsu.edu/~kaltofen/bibliography/index.html

Improve this answer
$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .