Determine if an operator is in the stabilizer group

Question

I was working on error correction algebraically but I now want to use some computational method. I have a set of generators of a stabilizer group which I represent as a rectangular matrix over $\mathbb{F}_2$ in a standard way. More precisely, if $i$'th generator $S_i$ contains $Z$ acting on $j$'th qubit, $(i,j)$ entry is $1$. If it contains $X$ acting on $j$'th qubit, $(i,j)$ entry is $1$.

Now, I have an operator $L$ represented as a vector in the same way, I want to see if it's contained in the group.

I tried using Python numpy and galois to calculate the pseudo-inverse of the matrix at first, which didn't work because of this reason. https://math.stackexchange.com/questions/1612708/computing-one-sided-inverse-of-a-matrix-over-some-finite-field The main problem is, over a finite field, even though the rows of the matrix is linearly independent, $(A A^T)^{-1}$, which is used to calculate the pseudo-inverse, is not necessarily invertible.

Next thing I tried was to use numpy.linalg.lstsq. Although galois claims that it supports combined use with numpy.linalg, it doesn't give an approximated solution over $\mathbb{F}_2$.

Answer 1

Never mind, I found a solution. I just asked a chatGPT and this efficiently gave me a solution.

import numpy as np

def solve_binary_linear_system(A, b):
    """
    Solve the linear system Ax = b over GF(2).

    Parameters:
    A (numpy.ndarray): Coefficient matrix (m x n)
    b (numpy.ndarray): Right-hand side vector (m)

    Returns:
    numpy.ndarray: Solution vector (n) if a solution exists, otherwise None
    """
    m, n = A.shape
    augmented_matrix = np.hstack((A, b.reshape(-1, 1)))

    # Perform Gaussian elimination in GF(2)
    for col in range(min(m, n)):
        # Find the pivot row
        pivot_row = None
        for row in range(col, m):
            if augmented_matrix[row, col] == 1:
                pivot_row = row
                break
        
        if pivot_row is None:
            continue
        
        # Swap the current row with the pivot row
        if pivot_row != col:
            augmented_matrix[[col, pivot_row]] = augmented_matrix[[pivot_row, col]]

        # Eliminate other rows
        for row in range(m):
            if row != col and augmented_matrix[row, col] == 1:
                augmented_matrix[row] ^= augmented_matrix[col]
    
    # Check for consistency and extract solution
    for row in range(m):
        if np.all(augmented_matrix[row, :-1] == 0) and augmented_matrix[row, -1] == 1:
            return None  # No solution
    
    # Back-substitution to get the solution
    x = np.zeros(n, dtype=int)
    for row in range(min(m, n) - 1, -1, -1):
        if augmented_matrix[row, row] == 1:
            x[row] = augmented_matrix[row, -1] ^ np.dot(augmented_matrix[row, row + 1:-1], x[row + 1:])
    
    return x

# Example usage
A = np.array([[1, 0, 1], [0, 1, 1]], dtype=int)
b = np.array([1, 0], dtype=int)
solution = solve_binary_linear_system(A, b)

print("Solution:", solution)
```