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.

Debbee

For "all of the above", try scipy.optimize.linprog in a new enough SciPy -- that actually implements the simplex algorithm. – Orthotropous

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.

Isometry

But I thought the bases of column A would be a subset of the columns of A that would be linearly independent. q above seems to be something else. – Debbee

q is a set of orthogonal vectors which span the column space of A. There are potentially infinitely many bases of the column space, q is an especially nice one. But if you need the basis to consist of columns of A, then you can compute the QR decomposition and throw out the linearly dependent columns. For example, see here. – Isometry

I just tested and found that presently it is not correct to use np.linalg.qr(A) to find a basis for the column space of A. This is because we may need all of the orthonormal vectors it provides through q to be able to produce A—even for singular matrices. Instead, we should use scipy.linalg.qr(A, pivoting=True) so that the first n columns of q will give us a basis. scipy.linalg.qr – Fa

Answer 2

Try using this

scipy.linalg.orth(A)

this produces orthonormal basis for the matrix A

Rideout