Generating Random Non-Singular Integer Matrices

8
votes

As part of a synthetic noise generation algorithm, I have to construct on the fly a lot of large non-singular square matrices

a i,j (i,j:1..n) / ∀ (i,j) a i,j ∈ ℤ and 0 ≤ a i,j ≤ k and Det[a] ≠ 0

but the a i,j should also be random following a uniform distribution in [0, k].

In its current incarnation the problem has n ≅ 300, k≅ 100.

In Mathematica, I can generate random element matrices very fast, but the problem is that I must also check for singularity. I am currently using the Determinant value for this.

The problem is that this check, for the 300x300 matrices takes something near 2 seconds, and I can't afford that.

Of course I may construct the rows by selecting a random first row and then constructing successive orthogonal rows, but I'm not sure how to guarantee that those rows will have their elements following an uniform distribution in [0,k].

I am looking for a solution in Mathematica, but just a faster algorithm for generating the matrices is also welcome.

NB> The U[0,k] condition means that taken a set of matrices, each position (i , j) across the set should follow a uniform distribution.

4
algorithmmathrandomwolfram-mathematica
Maybe you should post this question in a Gap or Sage forum. You basically want random matrices in GL(n,ZZ). For example in Sage you could use random_matrix(ZZ, n, n, algorithm='subspaces', rank=n), but for large n this is really slow... I'm sure the experts would be able to show you a better way. - Simon
@Simon I've a few alternatives. I am trying to avoid the problem by using random but Hermitian operators, but I've to prove some equivalence relations that are not easy. Also, I may go with your answer and see if the overall timing is good enough. And the GF(p) suggestion is also a possible way. That is all before going to Sage, which is one of my many areas of ignorance. - Dr. belisarius
You might find the article of Dana Randal useful eecs.berkeley.edu/Pubs/TechRpts/1991/CSD-91-658.pdf - Sasha
@Sasha Thanks! It seems much more efficient than the method I finally used. I'll keep the pointer in case I need to redo the filtering stage. - Dr. belisarius
For the record, Randall's paper is for nonsingular matrices over finite fields. A matrix of zeros and ones that has determinant "2" is singular over F_2, but nonsingular over the reals. - Rick

4 Answers

5
votes

In both Matlab and Octave, the determinant and LU factorization of a 500x500 matrix are basically instantaneous. Does Mathematica have a mode where it can call out to LAPACK or some similar library? You might need to annotate that your arrays should be treated as floating point numbers and not symbolically; that might make it a lot faster. As a point of comparison, LU on a 5000x5000 matrix takes 8.66 seconds on my system using Octave; 500x500 should be around 1000 times faster than that.

5
votes

You could use MatrixRank instead. On my machine it's about n/10 times faster for large nxn integer matrices.

3
votes

If you use numeric approximate matrices in the singularity tests you will get much better speed.

k = 100; n = 500;
mat = RandomInteger[100, {n, n}];

AbsoluteTiming[Det[mat] == 0]

Out[57]= {6.8640000, False}

AbsoluteTiming[Det[N@mat] == 0.] (*warning light!!*)

Out[58]= {0.0230000, False}

AbsoluteTiming[MatrixRank[N@mat] != n]

Out[59]= {0.1710000, False}

Unfortunately that fastest test is not reliable. But the rank test should work well. Here is a quick example wherein we replace the last row by the sum of prior rows.

mat2 = mat;
mat2[[-1]] = Total[Most[mat]];

AbsoluteTiming[Det[mat2] == 0]

Out[70]= {9.4750000, True}

AbsoluteTiming[Det[N@mat2] == 0.]

Out[69]= {0.0470000, False}

AbsoluteTiming[MatrixRank[N@mat2] != n]

Out[71]= {0.1440000, True}

In principle I suppose there is a smallish chance that the rank test could give a false negative., say due to ill conditioning. As your usage will better tolerate false positives (that is, incorrect claims of singularity) you could instead test for singularity over a prime modulus. I think this was one of the recommendations others have made.

Continuing the above examples:

AbsoluteTiming[Det[mat, Modulus -> Prime[1000]]]

Out[77]= {0.6320000, 4054}

AbsoluteTiming[Det[mat2, Modulus -> Prime[1000]]]

Out[78]= {0.6470000, 0}

This is slow but faster than working over the rationals. For what it's worth, for most usage I'd be fairly confident in the outcomes of the faster testing via MatrixRank[N[matrix]].

Daniel Lichtblau Wolfram Research

1
votes

Here's an expansion of a comment that I made a bit. I agree with Dan that it is highly improbable that the numerical version would return a false positive. Nonetheless, you can avoid that scenario by examining the singular values and returning False reliably, if the smallest singular value is larger than some error tolerance. (Admitedly, it might be a bit tricky to find a provable tolerance.) If the smallest singular value is uncomfortably small, you can apply Det to the integer matrix.

Here's a function that should quickly return False for most non-singular matrices. If the matrix is close to singular, a more expensive integer computation is performed.

singularQ[M_?MatrixQ] := If[
  Last[SingularValueList[N[M], Tolerance -> 0]] > 1/10.0^8,
 False, Det[M] == 0];

Here are 200 matrices that meet your description. One in the middle has been rigged to be singular.

SeedRandom[1];
matrices = RandomInteger[{0, 100}, {200, 300, 300}];
matrices[[100, -1]] = Sum[RandomInteger[{0, 10}]*matrices[[100, k]],
 {k, 1, Length[matrices[[100]]] - 1}];

Now, let's find the indices of all the singular matrices, watching as we go.

Flatten@Monitor[Table[
  If[singularQ[matrices[[k]]], k, {}],
  {k, 1, Length[matrices]}], k]