Partial pivoting in row reduction

Question

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.

Answer 1

Maybe this?:

tol = Sqrt[Max[Abs[approxmat - realmat]]]
{q, r} = QRDecomposition[approxmat];
RowReduce[r, Tolerance -> tol] //
  Rationalize[#, tol] & //
   MatrixForm

• The tolerance should depend on the noise, imo; Sqrt[] may be overaggressive (not sure, though).
• r is equal to q . approxmat; multiplying by q, which is invertible (orthogonal), preserves row-equivalence.