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$
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$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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