We can see this problem as an instance XOR-SAT problem, though it is more general than the problem posed here, that’s one way to go.
Just to gain some intuition I provide a very simple example. Suppose you have systems of three switches and three bulbs like this:
S    B
1    1, 3    // toggles bulbs 1 and 3
2    1, 2
3    1, 2, 3
It is equivalent to have the following formula, which we want to satisfy:
(x1^x2^x3)&(x1^x2)&(x1^x3).
And now we want to satisfy this formula. We start with writing it as system of boolean equations modulo 2:
 |1 1 1|   |x_1|   |1|
 |0 1 1| * |x_2| = |1| mod 2
 |1 0 1|   |x_3|   |1|
Now solve it with Gaussian elimination.
First, add the first and the second rows to the third:
1 1 1   1      1 1 1   1
0 1 1   1   -> 0 1 1   1
1 0 1   1      0 0 1   1   // for RHS, 1+1+1 = 1 mod 2
Second, back-substitute: x1 = 0, x2 = 0, x3 = 1, which is obviously the answer.
So, the main complexity here is to program Gaussian elimination process.
1
solved Can one turn all the light bulbs with m switches or not