Python - Get steps needed to convert matrix to RREF

Question

I have the following problem:

I have two binary matrices, that might look like this:

    | 1 | 0 | 1 | 0 | 0 | 1 |       | 1 | 0 | 0 | 0 | 0 | 0 |
    | 0 | 1 | 1 | 1 | 0 | 1 |       | 0 | 1 | 0 | 0 | 0 | 0 |
a = | 0 | 0 | 0 | 1 | 1 | 0 |   b = | 0 | 0 | 1 | 0 | 0 | 0 |
    | 1 | 0 | 0 | 0 | 0 | 1 |       | 0 | 0 | 0 | 1 | 0 | 0 |
    | 0 | 0 | 0 | 1 | 0 | 0 |       | 0 | 0 | 0 | 0 | 1 | 0 |
    | 1 | 0 | 1 | 0 | 0 | 1 |       | 0 | 0 | 0 | 0 | 0 | 1 |

Now i would like to find the row echelon form (not necesarily reduced row echelon form) of the matrix a, and then apply the same matrix operations on matrix b, which would result in something like this:

    | 1 | 0 | 1 | 0 | 0 | 1 |       | 1 | 0 | 0 | 0 | 0 | 0 |
    | 0 | 1 | 1 | 1 | 0 | 1 |       | 0 | 1 | 0 | 0 | 0 | 0 |
a = | 0 | 0 | 1 | 0 | 0 | 0 |   b = | 1 | 0 | 0 | 1 | 0 | 0 |
    | 0 | 0 | 0 | 1 | 1 | 0 |       | 0 | 0 | 1 | 0 | 0 | 0 |
    | 0 | 0 | 0 | 0 | 1 | 0 |       | 0 | 0 | 1 | 0 | 1 | 0 |
    | 0 | 0 | 0 | 0 | 0 | 0 |       | 1 | 0 | 1 | 0 | 1 | 1 |

Using numpy to convert the first matrix to rref works great, except I have no way of knowing what row operations were performed, so I can't apply the same operations on the second matrix.

Now this is just an example, but the actual matrices are going to be 50.000x50.000 or larger, and not necessarily square. I tried implementing my own solution, but it's just way too slow. Is there something out there, that can do what I want, or do I have to try to optimize my own solution?

Thanks for the help

/Morten

I'm rusty on linear algebra, have you checked this page to see if anything works? http://docs.scipy.org/doc/numpy/reference/routines.linalg.html — steven, Feb 18 '16 at 14:59

score 1 · Answer 1 · answered Feb 29 '16 at 14:20

Horizontally concatenate your matrices c = np.c_[a,b] and use rref on that matrix such that you don't need to store the intermediate matrix.

Note that rref is computationally useless except very very very special occasions. So prefer LU decomposition if you have to.

score 0 · Answer 2 · answered Feb 20 '16 at 01:13

Turns out that it's a bit messy, but there is a solution. scipy.linalg.qr https://en.wikipedia.org/wiki/QR_decomposition http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.qr.html#scipy.linalg.qr

import scipy.linalg as la
matrix = [[randint(2) for k in range(4)] for j in range(5)]
(q, r) = la.qr(matrix)

matrix:

[[1, 1, 1, 1], [0, 1, 1, 1], [0, 0, 0, 1], [1, 1, 1, 1], [1, 0, 0, 1]]

r:

array([[-1.73205081, -1.15470054, -1.15470054, -1.73205081],
       [ 0.        , -1.29099445, -1.29099445, -0.77459667],
       [ 0.        ,  0.        ,  0.        ,  1.        ],
       [ 0.        ,  0.        ,  0.        ,  0.63245553],
       [ 0.        ,  0.        ,  0.        ,  0.        ]])

numpy.dot(q,r):

array([[  1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
          1.00000000e+00],
       [  0.00000000e+00,   1.00000000e+00,   1.00000000e+00,
          1.00000000e+00],
       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          1.00000000e+00],
       [  1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
          1.00000000e+00],
       [  1.00000000e+00,   0.00000000e+00,   1.11022302e-16,
          1.00000000e+00]])

So matrix = q*r, and r is the row echelon form of matrix. All that is left to do is solve matrix2 = q*x for x. The rounding issue where 0 is not always exactly 0 is a known issue solving matrices numerically.

Python - Get steps needed to convert matrix to RREF

2 Answers2