9

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.

MikeRand
  • 4,788
  • 9
  • 41
  • 70
  • For "all of the above", try `scipy.optimize.linprog` in a new enough SciPy -- that actually implements the simplex algorithm. – Fred Foo Nov 27 '14 at 18:29

2 Answers2

11

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.

SiHa
  • 7,830
  • 13
  • 34
  • 43
jme
  • 19,895
  • 6
  • 41
  • 39
  • 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. – MikeRand Nov 27 '14 at 18:18
  • 2
    `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](http://math.stackexchange.com/a/748538). – jme Nov 27 '14 at 18:27
  • 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](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.qr.html#scipy.linalg.qr) – Shahrokh Bah Jan 19 '21 at 17:09
6

Try using this 

scipy.linalg.orth(A)

this produces orthonormal basis for the matrix A

Jayanth Kumar
  • 120
  • 1
  • 6