2

Suppose I have a symmetric, positive semidefinite matrix A of dimension n x n. The rank of A is d < n, and I want to find a set of indices of length d such that A[indices, indices] has rank d. For example, if

A <- structure(c(3, 2, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), .Dim = c(4L, 4L))

then I would like to identify one of the sets of indices c(1, 2, 3) or c(1, 2, 4).

I can do this by brute force going through all combinations of rows and columns, but my guess is that there exist much more elegant and efficient methods for doing this.

Lars Lau Raket
  • 1,905
  • 20
  • 35
  • Here's a scratch that should work: do `eigen(A)` and look at the eigenvector which has zero eigenvalue (i.e. the last column). It has non-zero components (mind numerical issues here) in the last two rows, i.e. one of the indices 3 and 4 can be dropped to arrive at your set of indices. Now as `n` gets larger there needs to be some combinatorics done but it seems doable. – mts Aug 10 '15 at 08:05
  • @mts: I have looked at this, and due to numerical precision (and very large matrices) it became very messy to use the eigenvalues to identify the proper columns. – Lars Lau Raket Aug 10 '15 at 09:06
  • edit: I found out that the `trim.matrix` method in the `subselect` package solves the problem using the eigenvalue approach, but it is considerably slower than the straightforward QR decomposition. – Lars Lau Raket Aug 10 '15 at 09:14

1 Answers1

3

You can use the QR decomposition, qr

q <- qr(A)
q$pivot[seq(q$rank)]
# [1] 1 2 3

A[,q$pivot[seq(q$rank)]]
#      [,1] [,2] [,3]
# [1,]    3    2    1
# [2,]    2    2    1
# [3,]    1    1    1
# [4,]    1    1    1

These resulting columns should be linearly independent.

Rorschach
  • 31,301
  • 5
  • 78
  • 129