# Suggestions for solving a large linear system

I have a system of $\approx 200,000$ linear equations in $\approx 40,000$ variables (with rational coefficients) and I would like to determine the dimension of the solution space, which I know to be positive (by other means).

The naive Solve does not work: the kernel runs out of memory and quits. One thing I tried is to break the system into 40 or so systems of $\approx 5,000$ equations each, to solve each subsystem and then to collate the solutions, but I wonder whether there is a more efficient way. In particular, applying the solutions already found to each subsystem in turn is quite slow.

Edit (based on the comments) The kernel crashes (a few seconds after evaluation) whenever I run CoefficientArrays even if the kernel appears to be "idle". The suggestion of applying N to the equations does speed things up considerably, but now the solutions look like

β2[0, 1, 1] ->
0. - 0.166667 β1[0, 1, 9] - 0.166667 β1[0, 2, 9] -
0.166667 β1[0, 3, 9] - 0.166667 β1[0, 4, 9] -
0.166667 β1[0, 5, 9] - 0.166667 β1[0, 6, 9] -
0.166667 β1[0, 7, 9] - 0.166667 β1[0, 8, 9] +
1.38778*10^-17 β1[0, 9, 0] - 0.166667 β1[0, 9, 9] -
0.166667 β1[0, 10, 9] + 1.38778*10^-17 β1[1, 9, 0] +
1.38778*10^-17 β1[2, 9, 0] + 1.38778*10^-17 β1[3, 9, 0] +
1.38778*10^-17 β1[4, 9, 0] + 1.38778*10^-17 β1[5, 9, 0] +
1.38778*10^-17 β1[6, 9, 0] + 1.38778*10^-17 β1[7, 9, 0] +
1.38778*10^-17 β1[8, 9, 0] - 1.38778*10^-17 β1[9, 10, 0]


whereas the exact solution is clearly

β2[0, 1, 1] -> - 1/6(β1[0, 1, 9] + β1[0, 2, 9] +
β1[0, 3, 9] + β1[0, 4, 9] + β1[0, 5, 9] + β1[0, 6, 9]
+ β1[0, 7, 9] + β1[0, 8, 9] + β1[0, 9, 9] + β1[0, 10, 9])


This has the unfortunate consequence that expressions which ought to be zero are just very small (but nonzero) and this seems to spoil my computation of the dimension of the solution space: an equation of the form $\epsilon \beta = 0$ where $\epsilon$ is very small but nonzero, would imply $\beta = 0$, whereas the exact result would have $\epsilon = 0$ and hence $\beta$ remains free.

Question Is there a way to force Mathematica to treat such small numbers as if they were zero? I don't expect very large denominators in the system, so it should be safe to treat as zero any number which is $\lt 10^{-4}$, say.

• Is the system sparse, by any chance? Commented Aug 3, 2015 at 14:37
• In that case I would recommend throwing a bit of N over them as this may have a dramatic impact on speed and memory use. Have a look at LinearSolve as well. Commented Aug 3, 2015 at 14:56
• Additionally, if your system is not already in matrix-vector format, you can use CoefficientArrays[] for the conversion. Commented Aug 3, 2015 at 15:00
• See Chop. It is provided for this specific purpose. Commented Aug 5, 2015 at 10:11
• I'll assume you can find a way to get at the underlying matrix. You want the dimension and not necessarily a solution. For the system mat.x==b (mat is the matrix, b the right side vector, x the vector of unknowns) this is equivalent to finding the dimension of the null space of the augmented matrix mat|b (using common linear algebra terminology, not Mathematica, to denote augmentation). Now check this prior MSE thread for possibilities. Commented Aug 5, 2015 at 15:54