1

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
marcotama
  • 1,991
  • 2
  • 19
  • 24

0 Answers0