I have a large $m \times 2m$ numerical matrix which I obtained from NullSpace
. I expect that the matrix is equivalent under row operations to a sparse matrix with small integer values.
However, when I RowReduce
, the matrix becomes ill-conditioned because it tries to pivot on very small values which should be zero up to numerical errors. Partial pivoting is the algorithm which pivots on the largest entry of a row at each step in row reduction, which would prevent the matrix from becoming ill-conditioned.
Is there an implementation of partial pivoting in Mathematica?
Edit 1: Here is some Mathematica code to generate examples of this kind.
(*Generate a random sparse integer matrix of size m x 2m*)
m = 20;
intmat =
SparseArray[
Flatten[Table[{j, RandomInteger[{1, 2 m}]} ->
RandomChoice[{-1, 1}], {i, 1, 5}, {j, 1, m}]]]
(*Multiply on the left by a random real unit determinant matrix.
This preserves the row span but obscures the sparse integer
representation*)
realmat =
Module[{randommat = RandomReal[{-1, 1}, {m, m}]},
randommat/Abs[Det[randommat]]^(1/m)] . intmat;
(*To recover an integer representation, row reduce. The row reduction
agrees with the row reduction of the original integer matrix, as
expected*)
Rationalize[Chop[RowReduce[realmat]]] // MatrixForm
RowReduce[intmat] // MatrixForm
(*Introduce small numerical errors, and then observe that the row
reduction is quite bad*)
approxmat = RandomReal[{-2*10^-10, 2*10^-10}, {m, 2 m}] + realmat;
Rationalize[Chop[RowReduce[approxmat]]] // MatrixForm
(*However, this is fixed by testing for zeroes with a high enough
tolerance*)
Rationalize[
Chop[RowReduce[approxmat, ZeroTest -> (Chop[#1, 10^-8] == 0 &)],
10^-8], 10^-7] // MatrixForm
So ZeroTest
fixes this problem perfectly for small matrices. However, try the same approximate example above with m=500. I find that the errors accumulate over the course of RowReduce
, requiring a larger and larger tolerance for ZeroTest
. In my case, my matrix has around $m=2000$. (Beware: running RowReduce
with ZeroTest
for $m=2000$ takes about 10 minutes on my 16 core machine.) For $m$ large enough, there is no value of the tolerance which gives the correct answer, once the accumulated errors become larger than the smallest nonzero values of the real matrix.
In my case, the solution will be to re-generate the original data at higher precision. This should allow me to get a clear separation between the size of the errors (even as they accumulate) and the smallest nonzero values. Even though this method is computationally costly, I'm lucky that I have this option at all.
I wonder if there would be any more stable way to analyze large real matrices that are hiding simple integer information. In the scenario where the uncertainties in the data are irreducible, such a stable analysis would be required.
Chop
to remove the negligible values prior toRowReduce
? $\endgroup$