Re: [eigen] fast calculation of matrix rank

Hi Marco,

Thanks for the reply. A_i is the partioned matrix [A | B_i]. I have succeeded in using Eigen to make my rank calculations faster (factor of 2 to3), which was worth the effort for me.

You are right that in general rank(AB_i) <= min(rank(A), rank(B_i)). In my case however, I have used the transformation with Q and P such that A is orthogonal to

[ 0 ]

[ Y_i ]

Then rank(X_i) = rank( [A [ 0; Y_i ] ]) = rank( A) + rank [ 0; Y_i ] ) = rank( A) + rank [Y_i ] )

Taking the transpose of Q and P is indeed faster than the inverse. The reason for taking Q^{-1} and not Q^T in my code is that I only calculate Q^{-1} once (so computation overhead is small) and I was worried that Q^T might lead to larger rounding errors than Q^{-1}. I am not sure whether these rounding errors can occur (for P they cannot since it has integer entries).

With kind regards,

Pieter Eendebak

On Sat, Jan 30, 2016 at 2:39 AM, marco stamazza <giunghi@xxxxxxxxx> wrote:

Pieter,
probably you are better off asking this question at http://math.stackexchange..com/questions/tagged/linear-algebra

You mean that A_i is a partitioned matrix [A | B_i] or a product A * B_i?

what is the size of the matrices A and B_i?

Note that the rank of product and partitioned matrices is rather more complicated than the sum of the ranks you use (http://link.springer.com/chapter/10.1007%2F978-3-642-10473-2_6#page-1)
for example rank(AB_i) <= min(rank(A), rank(B_i)), not sure where the additive property you state comes from.

In any case, I don't understand why you are computing Q^(-1) A_i P^(-1) i, since Q iand (thin) P are orthogonal only need take the transpose. Furthermore, since AP = QR, Q^(-1) A_i P^(-1) = Q^(-1) [A B_i] P^(-1) = RB_iP' (or similar for partitioned matrix).. So the transformation doesn't help you, on the contrary, R is already available with its rank, so you still need determine the rank of (B_i)..

If A_i = [A | B_i] then you could compute the rank of A and the the ranks of only the B_i's with colPivQR. but you still have to deal with the formula in the reference above

If A_i is a product A*B_i, a strategy could be computing If A_iA_i' and then get the rank by computing only its eigenvalues or with LU decomposition, which could be faster then colPivQR, but I am not sure.

In my opinion with such small nrow and ncol there is no point in inventing algorithms to save few flops, unless you have some properties of the A and B_i matrices you can exploit. .

I hope this helps, Marco

On 30 January 2016 at 05:24, Pieter Eendebak <pieter.eendebak@xxxxxxxxx> wrote:
Hi,

Using Eigen I am trying to calculate the matrix rank of a large number of arrays. So far I have used FullPivLU but this is too slow for my application. The alternative ColPivHouseholderQR is faster, but still not fast enough. What are options for calculating the rank faster? My matrices have size 64x92 and integer elements (+1 or -1).

There is one additional fact that might used. All matrices are of the form

A_i = [A B_i]

with A a common submatrix of the A_i.

I can then calculate the rank of A_i in the following way. First perform ColPivHouseholderQR on A:

AP = QR

and determine the matrix rank of A as r0. Then for each matrix A_i we calculate

X_i = Q^(-1) A_i P^(-1)

The rank of X_i is equal to the rank of A_i. Next, take the submatrix

Y_i = X_i[r0:, r0:]

Then rank(A_i) = rank(X_i) = rank(A) + rank(Y_i).

Instead of calculating the rank of A_i we only have to calculate the rank of Y_i, which is much smaller in size than A_i. I have checked that calculation of the ColPivHouseholderQR decomposition of the Y_i matrices is indeed faster. However, the multiplication Q^(-1) A_i P^(-1) is almost as expensive as the decomposition on the original array X_i. Is there a way I can use the common submatrix in the rank calculations without performing the multiplication?

With kind regards,
Pieter

The matrices share a common block Most of matrices

--
"Anyone who lives within their means suffers from a lack of imagination."