Scipy - find bases of column space of matrix

Question

I'm trying to code up a simple Simplex algorithm, the first step of which is to find a basic feasible solution:

1. Choose a set B of linearly independent columns of A
2. Set all components of x corresponding to the columns not in B to zero.
3. Solve the m resulting equations to determine the components of x. These are the basic variables.

I know the solution will involve using scipy.linalg.svd (or scipy.linalg.lu) and some numpy.argwhere / numpy.where magic, but I'm not sure exactly how.

Does anyone have a pure-Numpy/Scipy implementation of finding a basis (step 1) or, even better, all of the above?

Example:

>>> A
array([[1, 1, 1, 1, 0, 0, 0],
       [1, 0, 0, 0, 1, 0, 0],
       [0, 0, 1, 0, 0, 1, 0],
       [0, 3, 1, 0, 0, 0, 1]])

>>> u, s, v = scipy.linalg.svd(A)
>>> non_zero, = numpy.where(s > 1e-7)
>>> rank = len(non_zero)
>>> rank
4
>>> for basis in some_unknown_function(A):
...     print(basis)
{3, 4, 5, 6}
{1, 4, 5, 6}

and so on.

Answer 1

A QR decomposition provides an orthogonal basis for the column space of A:

q,r = np.linalg.qr(A)

If the rank of A is n, then the first n columns of q form a basis for the column space of A.