Generating Random Non-Singular Integer Matrices

问题:

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.

回答1:

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

回答2:

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.

回答3:

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

回答4:

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]