2
$\begingroup$

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.

$\endgroup$
9
  • 2
    $\begingroup$ Can you use Chop to remove the negligible values prior to RowReduce? $\endgroup$
    – Bob Hanlon
    Feb 16, 2023 at 16:15
  • 2
    $\begingroup$ RowReduce has an option "ZeroTest". Look it up in the help. $\endgroup$ Feb 16, 2023 at 16:20
  • 2
    $\begingroup$ Is there an implementation of partial pivoting in Mathematica? Could you make a MWE? I would be very surprised if Mathematica does not do this automatically. Yes, it does. Looking at help/implementation notes, it says RowReduce use Gaussian elimination with partial pivoting So I think you have something else going on. $\endgroup$
    – Nasser
    Feb 16, 2023 at 16:20
  • $\begingroup$ Ate these computations done with exact integer arithmetic? Or with approximate values. I would expect different behavior for these cases. $\endgroup$ Feb 17, 2023 at 3:22
  • 1
    $\begingroup$ Also do you know in advance how the rows were orthonormalized? If by Gram-Schmidt, there might be a possibility of undoing the operations. $\endgroup$ Feb 17, 2023 at 18:56

1 Answer 1

2
$\begingroup$

Maybe this?:

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

Mathematica graphics

  • 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.
$\endgroup$
2
  • $\begingroup$ Thanks for the answer! I'll try QR and see how it goes. Unfortunately in my actual application, I don't have a way to determine precisely how large the errors are. In other words, I don't have access to intmat or realmat, just approxmat, since the data comes with some amount of uncertainty/error. $\endgroup$ Feb 19, 2023 at 16:21
  • $\begingroup$ @DanielLongenecker "I don't have a way to determine precisely how large the errors are" -- maybe you can guess or approximate the tolerance. As I understand the problem, you're going to have to do it no matter what approach you take, no? You picked 2*10^-10 in the question -- maybe that was a reasonable choice. If you lower the tolerance too much, you get large-term fractions for numbers close integers from Rationalize and small nonzero numbers from Chop. $\endgroup$
    – Michael E2
    Feb 19, 2023 at 17:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.