Row reducing a matrix (not using sympy.rref)

Question

I have a matrix, G (size PxN), and I would like to make the first N/2 columns of G form the N/2×N/2 identity matrix by row operations. Any idea how I could write code for this? (not using sympy.rref).

Note P is a prime number slightly larger than N.

The matrix itself is a vandermonde matrix over the field of P i.e. all values are modulo P.

def ReedSolomon(k,p):
    rr = np.arange(k).reshape((-1, 1))
    cc = np.arange(p)
    return cc**rr % p

An example would be:

ReedSolomon(4,5) 
=          [[1, 1, 1, 1, 1],
           [0, 1, 2, 3, 4],
           [0, 1, 4, 4, 1],
           [0, 1, 3, 2, 4]]

And i would like to produce the following:

=          [[1, 0, -1, -2, -3],
           [0, 1, 2, 3, 4],
           [0, 1, 4, 4, 1],
           [0, 1, 3, 2, 4]]

In this case the N/2 x N/2 submatrix is the identity.

Using sympy.rref would lead to some rows being swapped around and for more rows to be reduced. I would like to stop once the N/2 columns have been turned into the identity. Hence writing custom code is preferred.

Thanks!

Answer 1

The following code works (adapted from: python built-in function to do matrix reduction)

def rref(V, tol = 1.0e-12):
    m, n = V.shape
    i, j = 0, 0
    jb = []
    while i < (n//2) and j < (n//2):
        # Find value and index of largest element in the remainder of column j
        k = np.argmax(np.abs(V[i:n, j])) + i
        p = np.abs(V[k, j])
        if p <= tol:
            # The column is negligible, zero it out
            V[i:n//2, j] = 0.0
            j += 1
        else:
            # Remember the column index
            jb.append(j)
            if i != k:
                # Do not swap the i-th and k-th rows
                pass
            # Divide the pivot row i by the pivot element A[i, j]
            V[i, j:n] = V[i, j:n] / V[i, j]
            # Subtract multiples of the pivot row from all the other rows
            for k in range(n//2):
                if k != i:
                    V[k, j:n] -= V[k, j] * V[i, j:n]
            i += 1
            j += 1
    return V

The adaptations include removing any row switching capability and limiting to the first n//2 columns and rows.

If anybody knows a more efficient way please let me know. Thanks