14
$\begingroup$

I am building a matrix, but some of the lines are dependant. I would like to only keep independant lines (so erase all lines that depend from another).

I could compute the nullspace of the transpose matrix and use this information to detect which lines depends on which.

But I am wondering if mathematica already implements such a method.

I am looking, by order of priority :

  1. The fastest way to do it : I am doing this in a loop and I need to avoid to take too much time.
  2. The simplest way to write it.
$\endgroup$

3 Answers 3

13
$\begingroup$

there is of course no unique way to do this.

MatrixForm[m = RandomInteger[5, {20, 3}]]

enter image description here

While[Length[null = NullSpace[Transpose[m]]] > 0,
   m = Drop[m, 
   RandomChoice[Position[null[[1]], Except[0], {1}, Heads -> False]]]]

enter image description here

or maybe better to build up a basis like this:

Fold[ Function[{t},
      If[Length@NullSpace@Transpose@t == 0, t, #1]]@ 
     Append[#1, #2] &, { First@m }, Rest@m ] // MatrixForm

enter image description here

$\endgroup$
10
$\begingroup$

Interesting problem! Computing MatrixRank is just the wee teeniest bit faster than NullSpace, so use that. First compute the number of rows needed (using @george2079 random matrix).

MatrixForm[m = RandomInteger[5, {20, 3}]];

numRows = Length[m];
minRows=MatrixRank[m];

Then randomly sample without replacement for a subset of rows with the same rank. Reasonably, this ought to find a properly ranked submatrix pretty quick, while guarding against any odd structure in the original matrix (e.g., blocks of duplicate rows).

While[MatrixRank[m[[rows = RandomSample[Range[numRows], minRows]]]] < minRows];

ans=m[[rows]]

enter image description here

A nasty matrix could make the While loop take a...while.

badm = Join[Table[{1, 1, 1, 1}, {10}], {{1, 0, 0, 0}}]

enter image description here

cnt = 1; 
While[
 MatrixRank[
   badm[[rows = RandomSample[Range[numRows], minRows] // Sort]]] < minRows, cnt++];
badm[[rows]] // MatrixForm
cnt

enter image description here

   b (* cnt = 10 *)
$\endgroup$
7
$\begingroup$

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];
$\endgroup$

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