My motivation is solving arbitrary Lights Out puzzles with $n^2$ cells with a method given in MathWorld: basically imagine the game as a square grid graph and construct an $n^2 \times n^2$ "adjacency matrix" $A$. Then for initial configuration $b$, solving $Ax = b$ gives you button presses to solve the game.
For a matrix of size $n^2 \times n^2$, normally Gaussian elimination takes $O((n^2)^3) = O(n^6)$ time, but maybe having a sparse binary matrix (by binary I mean all operations are mod 2, over field $\mathbb Z / 2 \mathbb Z$) can make the operation faster. I know the matrix is very sparse because every cell only has at most 4 neighbors, so almost all the matrix entries are zero. Any references to algorithms or software implementations that can solve binary sparse matrices efficiently are appreciated.
Lights Out solvers are already known, see, for instance, here or here, so I assume you plan to deal with big puzzle boards.
But the dimension of the matrix to deal can be reduced to $O(n)$ by the following row-by-row solving algorithm. Let the board has $n_1$ rows and $n_2$ columns. Suppose that $n_2\le n_1\ge 2$. Then $n_2\le n$.
Given a lit instance, going down from the second to the bottom rows and switching cells just below lit cells in the previous row, we switch off all lights but, possible, some of the bottom row. Let $b’$ be the binary vector corresponding to the lit instance of the bottom row.
Now we shall look for a binary vector $x_1$ with the following property. Assume that all lights at the board are switched off. Now switch the lights in the first row according to the vector $x_1$. Again go down from the second to the bottom rows and switch cells just below lit cells in the previous row. For each natural $i$ from $1$ to $n_1$ let $x_i$ be a vector corresponding to lights switched in $i$-th row. Put $x_0=0_{n_1}$. It is easy to see that for each natural $i$ from $1$ to $n_1-1$, $x_{i+1}=Bx_i+x_{i-1}$, where $B=\|b_{ij}\|$ is an $n_2\times n_2$-matrix such that $b_{ij}$ equals $1$ if $|i-j|\le 1$, and equals $0$, otherwise. Then $x_i=B_ix_1$ for each $i$, where $B_0=0$, $B_1=I$, and $B_{i+1}=BB_{i}+B_{i-1}$ for each $1\le n_1-1$. It follows that $B_i=f_i(B)$, where $f_i$ is the $i$-th Fibonacci polynomial defined over $\Bbb Z_2$ by $f_0=0$, $f_1=1$, and $f_k=xf_{k-1}+f_{k-2}$ for $k\ge 2$, see Theorem 2 in [GCT]. Our aim is to obtain in the bottom row a lit instance corresponding to $b$, that is to achieve $Bx_{n_1}+x_{n_1-1}=b$. This means that $x_1$ is a solution of an equation $Ax=b$, where $A=f_{n_1+1}(B)=BB_{n_1}+B_{n_1-1}$.
Remark that since the matrix $A$ does not depend on $b$, it can be calculated once for a board, at a preprocessing step. Moreover, a calculation of the polynomial $f_{n_1+1}$ can be shortcut using the equalities from Lemma 4 of [GCT].
References
[GCT] J. Goldwasser, W. Klostermeyer, and G. Trapp, Characterizing Switch Setting Problems, Linear and Multilinear Algebra, 43:1–3 (1997) 121–136.