How can I solve an underdetermined linear system?

[abdullahkhalids]

Suppose I have $m$ equations in $n$ variables for $m < n$. How can I solve such a system using galois? Usually there will be many solutions. I am fine with any valid solution to the system.

Here is an example, where we need a $x$ such that $Ax = b$.

import galois
GF = galois.GF(2)
A = GF([[1, 1, 1, 0],[0,0,1,1]])
b = GF([0,1])

#solve?

[mhostetter]

Following these links, I did the following. I did this in $\mathrm{GF}(7)$ so it was more clear -- things get obfuscated in $\mathrm{GF}(2)$.

• https://en.wikipedia.org/wiki/Underdetermined_system
• https://math.stackexchange.com/questions/1128047/solving-an-underdetermined-system-of-equations

import numpy as np
import galois

GF = galois.GF(7)
A = GF([[3, 0, 5, 2], [4, 1, 3, 0]])
print(f"A = \n{A}")

b = GF([6, 1])
print(f"b = \n{b}")

Ab = np.hstack((A, np.atleast_2d(b).T))
print(f"Ab = \n{Ab}")

# Must be true for solutions to be consistent
print(np.linalg.matrix_rank(A) == np.linalg.matrix_rank(Ab))

Abrr = Ab.row_reduce()
print(f"Abrr = \n{Abrr}")

for _ in range(10):
    x3 = GF.Random()  # Independent variable
    x4 = GF.Random()  # Independent variable
    x1 = -GF(4) * x3 - GF(3) * x4 + GF(2)  # Dependent variable
    x2 = -GF(1) * x3 - GF(2) * x4 + GF(0)  # Dependent variable
    x = GF([x1, x2, x3, x4])
    print(f"x = {x}")

    print(np.array_equal(A @ x, b))

From the row-reduced augmented [A | b] matrix, you can read it the relationships between the variables. Below are 10 valid solutions to the underdetermined equation.

A = 
[[3 0 5 2]
 [4 1 3 0]]
b = 
[6 1]
Ab = 
[[3 0 5 2 6]
 [4 1 3 0 1]]
True
Abrr = 
[[1 0 4 3 2]
 [0 1 1 2 0]]
x = [2 3 6 6]
True
x = [4 4 3 0]
True
x = [6 0 4 5]
True
x = [5 6 1 0]
True
x = [0 2 2 5]
True
x = [6 1 6 0]
True
x = [5 3 2 1]
True
x = [2 4 1 1]
True
x = [6 1 6 0]
True
x = [3 1 3 5]
True

[abdullahkhalids]

Thank you for your quick answer. It gave me the necessary structure to be able to write a function to do the necessary. This seems to work.

I would appreciate if you glance at this and see if there are any obvious mistakes.

import numpy as np
import galois

def solve_underdetermined_system(A, b):
    n_vars = A.shape[1]
    A_rank = np.linalg.matrix_rank(A)
    
    # create augmented matrix
    Ab = np.hstack((A, np.atleast_2d(b).T))

    # Must be true for solutions to be consistent
    if A_rank != np.linalg.matrix_rank(Ab):
        return False, None

    # reduce the system
    Abrr = Ab.row_reduce()
    
    # additionally need reduced row echelon form (RREF)
    swaps = []
    for i in range(Abrr.shape[0]):
        if Abrr[i,i] == 0:
            for j in range(i+1,n_vars):
                if Abrr[i,j] == 1:
                    Abrr[:, [i,j]] = Abrr[:, [j,i]]
                    swaps.append((i,j))
                    break

    # randomly generate some independent variables
    n_ind = n_vars - A_rank
    ind_vars = A.Random(n_ind)

    # compute dependent variables using reduced system and dep vars
    dep_vars = -Abrr[:A_rank,A_rank:n_vars]@ind_vars + Abrr[:A_rank,-1]
    
    # x is just concatenation of the two
    x = np.hstack((dep_vars, ind_vars))
    
    # Need to swap the variables according to column swaps for RREF
    for s in reversed(swaps):
        x[s[0]], x[s[1]] = x[s[1]], x[s[0]]
    
    return True, x

p = 5
K = 3
n = 5

GF = galois.GF(p)
A = GF.Random((K,n))
b = GF.Random(K)


print(f"A = \n{A}")
print(f"b = \n{b}")

solved, x = solve_underdetermined_system(A, b)

if solved:
    print(f"x = {x}")
    print("Solution correct:", np.array_equal(A @ x, b))
else:
    print("solution does not exist.")
    
### Now check lots of cases
noexist = 0
wrong = 0

GF = galois.GF(p)
for i in range(10000):
    A = GF.Random((K,n))
    b = GF.Random(K)
    solved, x = solve_underdetermined_system(A, b)
    
    if solved:
        sol_check = np.array_equal(A @ x, b)
        if not sol_check:
            wrong += 1
    else:
        noexist += 1

print(f'{noexist=}')
print(f'{wrong=}')

A = 
[[1 0 3 1 1]
 [0 3 1 2 2]
 [4 3 4 2 3]]
b = 
[2 1 1]
x = [1 3 2 0 0]
Solution correct: True
noexist=88
wrong=0

[mhostetter]

Looks pretty good to me.

One thing to note: Ab.row_reduce() already returns the reduced row echelon form. So no need to recompute it.

[abdullahkhalids]

I seem to recall that RREF meant that the columns must be swapped as well to put the identity matrix on the left (or right). Wikipedia disagrees with me on these semantics.

The Ab.row_reduce() operation does not put the identity matrix on the left. While not necessary, moving the identity to the left helps with indexing when we compute the dependent variables. That's why I did it.

P.S. It would be nice to have a cleaned up version of this in the library.