# Reduce a huge low-rank matrix

I have a huge square matrix, with a lot of zero columns and zero rows: $m=\left( \begin{array}{ccccccccc} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 3 & 7 & 9 & 0 & 0 & 0 \\ 0 & 0 & 0 & 4 & 8 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 2 & 6 & 3 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \end{array} \right)$

In general, I don't know where the zero columns (and rows), but when $i$-th row is empty, then $i$-th column is is also zero.

My question is: Is there any internal command to delete all zero columns and zero rows?

My motivation for this question is: I want to evaluate the eigenvalues of this huge matrix, which is slow. If I can reduce the matrix efficiently, Eigenvalues[] can be sped up.

• You can use DeleteCases first on the matrix, then on the transpose of the result, then transpose back. Nov 21 '13 at 1:16
• Transpose@DeleteCases[Transpose@DeleteCases[m, {0 ..}], {0 ..}] Nov 21 '13 at 1:17
• Slow, but... ImageData@ImageCrop[Image@m] :)
– rm -rf
Nov 21 '13 at 2:43

I might have found a way which is superior to the approach of deleting the rows and columns. It works as follows: You construct from the ArrayRules a new matrix, simply by re-numbering the matrix entries.

The implementation is straight forward

transformRules[m_] :=
Module[{idX, idY, r = ArrayRules[m]},
idX = DeleteDuplicates@Cases[r, ({x_Integer, _} -> _) :> x];
idY = DeleteDuplicates@Cases[r, ({_, y_Integer} -> _) :> y];
{idX, idY} = Thread[# -> Range[Length[#]]] & /@ {idX, idY};
SparseArray[r /. ({x_, y_} -> n_) :> ({x /. idX, y /. idY} -> n)]
]


## Details

Let's use a small example

m = ArrayPad[Partition[Range[9], 3], 2];
ArrayRules[m]

(* {{3, 3} -> 1, {3, 4} -> 2, {3, 5} -> 3, {4, 3} ->
4, {4, 4} -> 5, {4, 5} -> 6, {5, 3} -> 7, {5, 4} -> 8, {5, 5} ->
9, {_, _} -> 0} *)


What you see are the indexes for the non-zero elements. When we look only at the x-positions, we find that we have entries in

{3,3,3,4,4,4,5,5,5}


which means that there are only 3 distinct columns having non-zero entries: {3,4,5}. Now it doesn't matter where exactly those columns are, so we will simply rename them to {1,2,3} which puts them at the very beginning of our matrix (no matter how big the matrix is). Note that missing columns would simply be erased by this. So when we had elements at {4,23,45}, the result after renaming would still be {1,2,3}.

The two lines in the algorithm starting with idX = .. and idY = .. create exactly this list of duplicate-free, original column and row numbers. The next line starting with {idX,idY} = .. computes the renaming to id's starting with 1. After that we take the original array-rules and replace the old positions with our new ones and then we create a SparseArray from them.

Taking the initial m we get the following new rules:

transformRules[m] // ArrayRules

(* {{1, 1} -> 1, {1, 2} -> 2, {1, 3} -> 3, {2, 1} ->
4, {2, 2} -> 5, {2, 3} -> 6, {3, 1} -> 7, {3, 2} -> 8, {3, 3} ->
9, {_, _} -> 0} *)


The computation of Eigenvalues is invariant under this transformation and since the computation of ArrayRules is so incredible fast, the whole algorithm is fast, because it basically replaces only a few positions in a matrix of some probably some million entries.

## Example

Using a big example:

m = ArrayPad[Partition[Range[9], 3], 5000];
m[[1, 1]] = -3;
m[[45, 23]] = 18;


we can now test this:

Eigenvalues[transformRules[m]]

(* {18, 3/2 (5 + Sqrt[33]), -3, 3/2 (5 - Sqrt[33]), 0} *)


This takes only 0.2 seconds here. If we compare this to the already made suggestions, we see:

Using Nassers method, you see that N@Eigenvalues[m] takes forever. Only the last line, where you restrict the number of Eigenvalues and use N inside is faster. It takes about 1.1 seconds here and has the disadvantages which comes from numeric, rather than exact calculations

Eigenvalues[N@m, 5, Method -> Arnoldi]

(* {16.1168, -3., -1.11684, 1.66277*10^-16 + 7.08214*10^-10 I,
1.66277*10^-16 - 7.08214*10^-10 I} *)


The method suggested by Ray gives us exact Eigenvalues and needs about 1.23 seconds here

Eigenvalues[Nest[Transpose@DeleteCases[#, {0 ..}] &, m, 2]]

(* {18, 3/2 (5 + Sqrt[33]), -3, 3/2 (5 - Sqrt[33]), 0} *)


You do not really need to delete rows and columns. Simply convert to Sparse matrix. Sparse matrix supports eigenvalue as is.

Since your matrix is full of zeros any way, this will be more efficient as well.

m0 = {{3, 7, 9}, {4, 8, 1}, {2, 6, 3}};
m = Table[0, {10}, {10}];
m[[4 ;; 6, 4 ;; 6]] = m0;
MatrixForm[m]


convert to sparse and find eigenvalues

msp = SparseArray[m];
N@Eigenvalues[msp]


It is the same eigenvalues as the dense matrix ofcourse

N@Eigenvalues[m]


You can also take advantage of this method:

For many large sparse systems that occur in practical computations the Arnoldi algorithm is able to converge quite quickly.

 Eigenvalues[N@msp, 3, Method -> Arnoldi]


LinearAlgebraMatrixComputations.html

WorkWithSparseMatrices.html

• I am using a $1024\times 1024$ real matrix as a SparseArray, with 999 zero columns/rows. With the reduced matrix, the Eigenvalues[] computes faster by 3 to 4 times.
– 9527
Nov 21 '13 at 2:35

See Deleting entire row and/or column. The fastest solutions there are

Nest[Transpose@DeleteCases[#, {0 ..}]&, m, 2]     (* integers *)
Nest[Transpose@DeleteCases[#, {0. ..}]&, m, 2]    (* reals *)
Nest[Transpose@DeleteCases[#, {(0.|0)..}]&, m, 2] (* mixeds *)

• This answers my question. Originally I also expected other solutions like: is there any command to reduce/decompose a $10 \times 10$ square matrix with matrix rank equal to 3, to a $3\times 3$ matrix, such that the nonzero eigenvalues are not changed during the reduction/decomposition?
– 9527
Nov 21 '13 at 2:40

Rather than deleting cases, you can also directly select the rows and columns that are nonzero. For instance, with m as defined in Nasser's answer:

n = Select[m, # != ConstantArray[0, Length[m]] &];
mat = Select[Transpose[n], # != ConstantArray[0, Length[n]] &];


Then the eigenvalues of mat are as above.

• This method is faster than DeleteCases method, and is as fast as halirutan's solution.
– 9527
Nov 21 '13 at 3:02

Sparse matrices

Should your matrices be SparseArrraya, then this method works quite well:

Creating a random SparseArray.

n = 100000;
m = 200000;
k = 1000000;
A = SparseArray[
Rule[
Transpose[{RandomInteger[{1, m}, {k}],RandomInteger[{1, n}, {k}]}],
RandomReal[{-1, 1}, {k}]
], {m, n}, 0.];


This looks up the rows and columns with at least one nonzero entry and extracts them corresponing matrix. B = SparseArray[A] in the beginning guarantees that no "NonzeroPositions" with (numerically) zero "NonzeroValues" persist. (This is superfluous here since A was created by SparseArray, but in general, this might be necessary.)

Anew = With[{B = SparseArray[A]},
A[[Sequence @@ Map[Sort@*DeleteDuplicates, Transpose[B["NonzeroPositions"]]]]]
]; // RepeatedTiming


{0.073, Null}

If you would like to have a square matrix as output (which is desirable for computing Eigensystem), you can achieve that with

n = m = 100000;
k = 1000000;
A = SparseArray[
Rule[
Transpose[{RandomInteger[{1, m}, {k}],RandomInteger[{1, n}, {k}]}],
RandomReal[{-1, 1}, {k}]], {m, n},
0.];
Anew = With[{B = SparseArray[A]},
With[{ilist = Sort@DeleteDuplicates@Flatten[B["NonzeroPositions"]]},
A[[ilist, ilist]]
]
]; // RepeatedTiming


Dense matrices

For dense matrices that have few nonzero entries, it may be worth to convert them to SparseArrays. Otherwise, the following works (don't try that with the matrix from above!

n = 1000;
m = 2000;
k = 10000;
A = Normal[SparseArray[
Rule[
Transpose[{RandomInteger[{1, m}, {k}], RandomInteger[{1, n}, {k}]}],
RandomReal[{-1, 1}, {k}]
], {m, n}, 0.]];

Anew = With[{B = Abs[A]},
A[[
Flatten[Position[B.ConstantArray[1., Dimensions[B][[2]]], _?Positive]],
Flatten[Position[ConstantArray[1., Dimensions[B][[1]]].B, _?Positive]]
]]
]; // RepeatedTiming


{0.011, Null}