I have a matrix A (shape (P,Q)) and I need to select a set of rows I such that A[I, :] is square and invertible. Is there an algorithm that I don't know to compute it?
More in general, I need to split A in as many non-singular matrices with shape (*, Q) as possible. Is there an algorithm to compute this?
A (pointer to a) numpy implementation would be appreciated.
My current (greedy and brute-force) implementation is the following:
def independent_columns_index(A, tol=1e-05):
"""
Found here: http://stackoverflow.com/questions/13312498/how-to-find-degenerate-rows-columns-in-a-covariance-matrix
"""
A = np.array(A)
Q, R = np.linalg.qr(A)
independent = np.where(np.abs(R.diagonal()) > tol)[0]
return independent
def independent_rows_index(A, tol=1e-05):
return independent_columns_index(A.T)
def n_independent_rows_indices(A, n, tol=1e-05):
if n > A.shape[1]:
return
independent = np.empty(0, dtype=int)
next_unchecked_row = 0
need_more = n
while A.shape[0] - next_unchecked_row > need_more:
need_more = n - len(independent)
first_row = next_unchecked_row
last_row = next_unchecked_row + need_more
next_unchecked_row = last_row
indices = np.append(independent, np.arange(first_row, last_row))
subA = A[indices, :]
ret = independent_rows_index(subA)
independent = indices[ret]
if len(independent) == n:
yield independent
independent = np.empty(0, dtype=int)
def test():
A = np.random.randint(0, 100, [43, 5])
k = 20
A[-k:, :] = A[:k, :]
np.random.shuffle(A)
A[1, :] = A[0, :]
for ind in n_independent_rows_indices(A, 5):
R = A[ind, :]
print(ind)
print(np.linalg.matrix_rank(R))
if __name__ == '__main__':
test()
EDIT: question corrected after Amit's comment. Added my naive algorithm after John's comment.
EDIT: improved solution:
def n_independent_rows_indices(A, n):
if n > A.shape[1]:
return
not_used = np.arange(A.shape[0])
while True:
try:
A_ = np.array(A[not_used, :][:3*n, :]).T
_, _, P = sp.linalg.qr(A_, pivoting=True)
now_used = not_used[P[:n]]
if len(now_used) < n:
return
not_used = np.delete(not_used, P[:n], axis=0)
yield now_used
except:
return