Obtaining an invertible square matrix from a non-square matrix of full rank in numpy or matlab

Question

Assume you have an NxM matrix A of full rank, where M>N. If we denote the columns by C_i (with dimensions Nx1), then we can write the matrix as

A = [C_1, C_2, ..., C_M]

How can you obtain the first linearly independent columns of the original matrix A, so that you can construct a new NxN matrix B that is an invertible matrix with a non-zero determinant.

B = [C_i1, C_i2, ..., C_iN]

How can you find the indices {i1, i2, ..., iN} either in matlab or python numpy? Can this be done using singular value decomposition? Code snippets will be very welcome.

EDIT: To make this more concrete, consider the following python code

from numpy import *
from numpy.linalg.linalg import det

M = [[3, 0, 0, 0, 0],
     [0, 0, 1, 0, 0],
     [0, 0, 0, 0, 1], 
     [0, 2, 0, 0, 0]]
M = array(M)

I = [0,1,2,4]
assert(abs(det(M[:,I])) > 1e-8)

So given a matrix M, one would need to find the indices of a set of N linearly independent column vectors.

Answer 1

My first thought would be to try each possible combination of N out of the M columns. That could be done like this (in Python):

import itertools
import numpy.linalg

# 'singular' returns whether a matrix is singular.
# You could use something more efficient than the determinant
# (I'm not sure what options there are in NumPy)
singular = lambda m: numpy.linalg.det(m) == 0

def independent_square(A):
    N,M = A.shape
    for colset in itertools.combinations(xrange(M), N):
        B = A[:,colset]
        if not singular(B):
            return B

If you want the column indices instead of the resulting square matrix, just replace return B with return colset. Or you could get both with return colset,B.

I don't know of any way SVD would help here. In fact, I can't think of any purely mathematical operation that would convert A to B (or even any that would figure out the MxN column selection matrix Q such that B=A.Q), other than by trial and error. But if you want to find out whether one exists, math.stackexchange.com would be a good place to ask.

If all you need is a way to do it computationally, the above code should suffice.

Answer 2

Easy, peasy in MATLAB. Use QR, specifically, the pivoted QR.

M = [3 0 0 0 0;
     0 0 1 0 0;
     0 0 0 0 1; 
     0 2 0 0 0]

[Q,R,E] = qr(M)
Q =
     1     0     0     0
     0     0     1     0
     0     0     0     1
     0     1     0     0

R =
     3     0     0     0     0
     0     2     0     0     0
     0     0     1     0     0
     0     0     0     1     0

E =
     1     0     0     0     0
     0     1     0     0     0
     0     0     1     0     0
     0     0     0     0     1
     0     0     0     1     0

The first 4 columns of E designate the columns of M to be used, i.e., columns [1,2,3,5]. If you want the columns of M, just form the product M*E.

M*E
ans =
     3     0     0     0     0
     0     0     1     0     0
     0     0     0     1     0
     0     2     0     0     0

By the way, Using det to determine if a matrix is singular is the absolutely, positively, definitely worst way you can do that.

Use rank instead.

Essentially, you should virtually NEVER use det in MATLAB unless you understand why it is such a bad thing to do, and you choose to use it despite that fact.