I want to solve a linear matrix equation in mod 2, which takes the form AX = B but with additional restrictions described below. Here, columns of A are linearly independent to begin with. The RHS, B on the other hand, has some redundancy as some columns are all 0. X is unknown, but can be solved using known GF(2) solvers.
Clearly X is not a one-to-one transformation as it maps distinct columns in A to the same vector (the all 0 vector). So X is not invertible in this system.

Now, I place an additional condition: that A and B are also defined as binary matrices on a larger vector space (think more coordinates per column, and also more columns stacked horizontally).
In the bigger space, the equation now takes the form A' X' = B' .
Equivalently, you could think of A,X,B as A’,X’,B’ restricted to a subspace of coordinates.
Since A has linearly independent columns, A’ on the bigger space is guaranteed to be column-wise linearly independent as well.
This time, I also know that columns B’ are linearly independent.

On the LHS, both A and A’ are fully known. On the RHS, I only know B (on the restricted subspace) but not B’ on the bigger space. In fact B’ on the rest of the entries is unconstrained. I don’t really care as long as columns of B’ are some linear combinations of columns in A’, which is guaranteed by the equation A’X’ = B’.
I want to solve this system of equations and find X.

I am currently using Mathematica LinearSolve[A,B, Modulus->2] to solve just the AX = B equation.
Due to redundancies in B, there are multiple solutions on X and I want a way to write my additional constraint from the bigger subspace into this solver, and find X that preserves linear independence on the bigger space.

Answers in Sage, Python, or if needed Matlab also welcome!