7
$\begingroup$

I have a fairly large singular square matrix, of size 37 and rank 35.

Is there a way to get a vector with the (groups of) indexes of the columns in the original matrix that are linearly dependent?

RowReduce et al won't help me, because they perform an arbitrary number of row swaps.

$\endgroup$
1
  • $\begingroup$ But it doesn't perform column swap. $\endgroup$
    – user202729
    Commented Jan 3, 2020 at 8:46

2 Answers 2

7
$\begingroup$

One way is to start with empty matrix. Add the first column. Then loop, each time adding the next column, and checking if the rank of this matrix has increased from before, if so, keep it, else skip over to the next column. Keep doing this until you reach the last column in the original matrix, or have collected m columns, where m is the rank of the original matrix. (no need to keep trying if found m columns).

function

indepCols[mat_?(MatrixQ[#, NumericQ] &)] := Module[{nRows, nCols, m, idx, 
    i, vecs, candidate},

  {nRows, nCols} = Dimensions[mat];
  m = MatrixRank[mat];
  idx = Table[1, {m}]; (*arrary to collect the index of columns*)
  vecs = {mat[[All, 1]]}; (*first column is always in*)
  Do[
   candidate = Join[vecs, {mat[[All, i]]}];
   If[MatrixRank[candidate] > Length[vecs],(*did the rank increase?*)
    idx[[Length[vecs] + 1]] = i;
    vecs = candidate;
    If[Length[vecs] == m, Break[]](*bail out if got the rank*)
    ],
   {i, 2, nCols}
   ];

  {idx, vecs}
  ]

To use the above function:

First example

(mat = {{1, 0, -2, 1, 0}, {0, -1, -3, 1, 3}, {-2, -1, 1, -1, 3}, 
     {-2, -1, 1, -1, 3}, {0, 3, 9, 0, -12}}) // MatrixForm

Mathematica graphics

The above has rank 3. So there will be 3 L.I. columns

 {idx, out} = indepCols[mat];
 Transpose[out] // MatrixForm

Mathematica graphics

idx
(*{1, 2, 4}*)

Verified using "AdjacencyLists" thanks to Michael E2 above.

sa = Transpose@Unitize@SparseArray@RowReduce[mat];
row2col = Rule @@@ Reverse /@ First /@ GatherBy[sa["NonzeroPositions"], First];
Thread[Range[First@Dimensions@sa] -> (sa["AdjacencyLists"] /. row2col)]

Mathematica graphics

second example

Using example given by Michael E2, which is much larger

SeedRandom[1];
mat = RandomSample[#~Join~Accumulate@RandomSample[#, 2] &@
    RandomInteger[{-5, 5}, {35, 37}]];
{idx, out} = indepCols[mat];
idx

Mathematica graphics

Verified:

sa = Transpose@Unitize@SparseArray@RowReduce[mat];
row2col = Rule @@@ Reverse /@ First /@ GatherBy[sa["NonzeroPositions"], First];
Thread[Range[First@Dimensions@sa] -> (sa["AdjacencyLists"] /. row2col)]

Mathematica graphics

The above gives the index of the L.I. columns. To find the index of the L.D. columns, simply take the complement:

Complement[Range[1, Length[mat]], idx]
(*{36, 37}*)
$\endgroup$
2
  • $\begingroup$ Also very sweet. Thanks to the both of you. +1 for now until I can get my hands on my matrix again. $\endgroup$ Commented Nov 21, 2014 at 10:21
  • $\begingroup$ I ended up going with this solution, mostly because I understand it, given my relatively short knowledge of Mathematica advanced programming. $\endgroup$ Commented Nov 22, 2014 at 19:09
6
$\begingroup$

If RowReduce won't help, then perhaps I don't know what you're looking for. Here's my understanding of the question in which I use RowReduce to get the answer.

Example

A random matrix:

SeedRandom[1];
mat = RandomSample[#~Join~Accumulate@RandomSample[#, 2] &@
    RandomInteger[{-5, 5}, {35, 37}]];
MatrixRank[mat]
(*  35  *)

We can use the SparseArray property "AdjacencyLists" to find for each column, which columns it depends on. I'm not sure what form the final result should be in. The output below indicates for each column, which columns it is a linear of. An entry of the form i -> {i} indicates an independent vector. The entry 21 -> {3, 8} indicates that column 21 is a linear combinations of columns 3 and 8. The leading 1 in each nonzero row of the reduced matrix indicates the column with which the columns with nonzero entries in the row have a relationship. The set of rules row2col basically indicates for each row, the column in which the leading 1 is located.

sa = Transpose @ Unitize @ SparseArray @ RowReduce[mat2];
row2col = 
  Rule @@@ Reverse /@ First /@ GatherBy[sa["NonzeroPositions"], First];
Thread[Range[First @ Dimensions @ sa] -> (sa["AdjacencyLists"] /. row2col)]
(*
{1 -> {1}, 2 -> {2}, 3 -> {3}, 4 -> {4}, 5 -> {5}, 6 -> {6}, 7 -> {7},
 8 -> {8}, 9 -> {9}, 10 -> {10}, 11 -> {11}, 12 -> {12}, 13 -> {13}, 
 14 -> {14}, 15 -> {15}, 16 -> {16}, 17 -> {17}, 18 -> {18}, 
 19 -> {19}, 20 -> {20},
 21 -> {3, 8},
 22 -> {22}, 23 -> {23}, 24 -> {24},
 25 -> {3},
 26 -> {26}, 27 -> {27}, 28 -> {28}, 
 29 -> {29}, 30 -> {30}, 31 -> {31}, 32 -> {32}, 33 -> {33}, 
 34 -> {34}, 35 -> {35}, 36 -> {36}, 37 -> {37}}
*)
$\endgroup$
5
  • $\begingroup$ nice use of "AdjacencyLists", I had no idea such a thing exist. $\endgroup$
    – Nasser
    Commented Nov 21, 2014 at 3:23
  • $\begingroup$ @Nasser Mr.Wizard's answers make a nice course of study. $\endgroup$
    – Michael E2
    Commented Nov 21, 2014 at 3:35
  • $\begingroup$ @MichaelE2, are there more strings like "AdjacencyLists" and "NonzeroPositions", and if so, where can I read about them? I could not even find these two in the help page for SparseArray or the Tutorials "Sparse Arrays: Manipulating Lists/Linear Algebra". $\endgroup$ Commented Nov 21, 2014 at 7:14
  • $\begingroup$ This is good stuff. I'll test it on my matrix tomorrow, but it looks like exactly what I was looking for. +1 for now. $\endgroup$ Commented Nov 21, 2014 at 10:18
  • 1
    $\begingroup$ @MariusLadegårdMeyer Many data structures like SparseArray take an argument "Properties" and return a list of possible inputs. Try sa["Properties"]. I cannot always find out how to use them or even if they're implemented. Searching for the string + data structure on this site sometimes yields results. Or you can ask question here. Often they're not documented. $\endgroup$
    – Michael E2
    Commented Nov 21, 2014 at 11:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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