QRDecomposition
for the transpose of m
could help. Try this:
m = RandomInteger[5, {20000, 150}];
R = QRDecomposition[Transpose[N[m]]][[2]];
indices =
Union @@ Map[
Function[row,FirstPosition[row, _?(Function[x, Abs[x] > 10^-12]), {1}]],
R
];
m[[indices]] // MatrixForm
MatrixRank[m[[indices]]]
The leading nonzero entries in each row of R
tell you where to find a useful row in m
. I converted m
to floating-point in order to guarantee that an efficient LAPACK-routine is used as the backend. The tolerance 10^-12
was chosen arbitrarily. Maybe there is an even better integer routine out there...
Edit: My old solution had some trouble with the matrix
m = {{0, 0, 0}, {0, 1, 1}, {0, 0, 0}, {0, 1, 1}, {1, 0, 1}, {1, 0, 0}};
I somewhat assumed that R
were in row echelon form. I have just learnt about the command RowReduce
. One can directly apply RowReduce
to Transpose[N[m]]
in order to get to reduced row echelon form. This can be done as follows:
On["Packing"]
U = RowReduce[Transpose[N[m]]];
indicesU = Union @@ Map[
Function[row, FirstPosition[row, 1, {1}]],
U
]; // AbsoluteTiming
m[[indicesU]] // MatrixRank
This would be even faster if Mathematica would not break a PackedArray
on its way from LAPACK to the math kernel by applying something like UpperTriangularize
in the background. However, the pivots are easier to detect this way...
Edit: Recently, I discovered that also FirstPosition
unpacks arrays. A faster way to find the leading nonzero row entries of a matrix u
in row echelon form is provided by the following function.
getIndices = Compile[{{u, _Real, 2}},
Block[{i, j, dims, imax, bag, jmax},
dims = Dimensions[u];
imax = dims[[1]];
jmax = dims[[2]];
i = 1;
j = 1;
bag = Internal`Bag[Most[{0}]];
While[(i < imax + 1) && (j < jmax),
While[(Abs[u[[i, j]]] < 1. 10^-12) && (j < jmax),
++j
];
Internal`StuffBag[bag, j];
++i; ++j;
];
Internal`BagPart[bag, All]
],
CompilationTarget -> "C"
];
Usage is:
u = Developer`ToPackedArray[N[U]];
indices = getIndices[u];