title: How can I further optimize this solver of a variant of "Lights Out"?
body: <p>I'm still solving <a href="https://math.stackexchange.com/questions/1052609/minimal-number-of-moves-needed-to-solve-a-lights-out-variant">this problem</a>, taken from the current "Google Foobar" challenge. It's a variation of the "Lights Out" game, in which pressing a light will flip the state of every light on the same row and the same column.</p>
  
  <p>I previously tried <a href="https://stackoverflow.com/questions/27337694/avoid-a-deepcopy-when-doing-a-bfs">using a BFS</a>, which turned out to be too slow for n > 6, while I need to handle 2 &lt; n &lt; 16. I currently have a program that can handle all even n and all odd numbers except 13 and 15. Here's what it does:</p>
  
  <ol>
  <li><p>I use the <a href="https://math.stackexchange.com/a/441697/4471">strategy</a> outlined by @Aryabhata to find a special solution x' of some system Ax = b that can be associated with an instance of this problem (see <a href="https://math.stackexchange.com/a/441588/4471">here</a> for details).</p></li>
  <li><p>Having found a <a href="https://math.stackexchange.com/questions/1056944/characterize-the-nullspace-of-a-given-matrix-in-mathbbf-2n-times-n?lq=1">base</a> of the null space of A, I compute all sums of x' plus a linear combination of the vectors of the base.</p></li>
  <li><p>The set of those sums is the set of all solutions of the original problem, therefore I can find by brute-force the solution that achieves the minimum.</p></li>
  <li><p>It should be noted that, for n even, the null space is empty (A is invertible), therefore x' achieves the minimum because it's the only solution. If n is odd the number of vectors in a base of the null space is 2n - 2, therefore the search space has size 2^(2n - 2), which is 2^28 in the worst case (n = 15).</p></li>
  </ol>
  
  <p>Here's my program:</p>
  
  <pre><code>from itertools import product
  
  
  MEMO = {}
  
  
  def bits(iterable):
      bit = 1
      res = 0
      for elem in iterable:
          if elem:
              res |= bit
          bit &lt;&lt;= 1
      return res
  
  
  def mask(current, n):
      if (current, n) in MEMO:
          return MEMO[(current, n)]
  
      result = 0
      if current &lt; n:
          for j in xrange(n):
              result += (2 ** ((current - 1)*n + j) + 2 ** (current*n + j))
      else:
          for i in xrange(n):
              result += (2 ** (i*n + current - n) + 2 ** (i*n + current - n + 1))
  
      MEMO[(current, n)] = result
  
      return result
  
  
  # See: https://math.stackexchange.com/a/441697/4471
  def check(matrix, n):
      parities = [sum(row) % 2 for row in matrix]
      for i in xrange(n):
          parities.append(sum([row[i] for row in matrix]) % 2)
  
      return len(set(parities)) == 1
  
  
  def minimize(matrix, current, n):
      if current == 0:
          # See: https://stackoverflow.com/a/9831671/374865
          return bin(matrix).count("1")
      else:
          return min(minimize(matrix ^ mask(current, n), current - 1, n),
                     minimize(matrix, current - 1, n))
  
  
  def solve(matrix, n):
      result = [0 for i in xrange(n) for j in xrange(n)]
  
      for i, j in product(xrange(n), repeat=2):
          if matrix[i][j]:
              for k in xrange(n):
                  result[i*n + k] ^= 1
                  result[k*n + j] ^= 1
              result[i*n + j] ^= 1
  
      if n % 2 == 0:
          return sum(result)
      else:
          return minimize(bits(result), 2*n - 2, n)
  
  
  def answer(matrix):
      n = len(matrix)
  
      if n % 2 == 0:
          return solve(matrix, n)
      else:
          if check(matrix, n):
              return solve(matrix, n)
          else:
              return -1
  </code></pre>
  
  <p>I've already tried optimizing it: for instance, matrices are encoded as binary numbers by the function <code>bits</code>, while the function <code>mask</code> creates binary masks that are used to add a single element of the base to x'. Those masks are also memoized because they are frequently used, so that they are calculated only once.</p>
  
  <p>The number of ones is then counted using the idiom <code>bin(n).count('1')</code>, which <a href="https://stackoverflow.com/questions/9829578/fast-way-of-counting-bits-in-python/9831671#9831671">should be</a> the fastest implementation (I checked it against the classical one by Kernighan).</p>
  
  <p>So, what else can I do to squeeze more performance out of my program? Here are a few test cases:</p>
  
  <pre><code>print answer([
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
  ]), 1
  
  print answer([
      [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
      [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0]
  ]), 14
  
  print answer([
      [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
      [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
      [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
      [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
      [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
      [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
      [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
      [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
      [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]
  ]), 15
  
  print answer([
      [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
      [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
      [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
      [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
      [1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
      [1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
      [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
      [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
      [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
      [1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
      [1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
      [1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
      [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
      [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
  ]), 14
  
  print answer([
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
      [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
      [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
      [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
      [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
      [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
      [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
      [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
      [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
      [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
  ]), 15
  </code></pre>
  
  <p><strong>EDIT:</strong> I passed this round. This implementation solves correctly 4 out of 5 test cases, then I brute-forced the fifth. I'm still interested in further optimizations or different algorithms!</p>
  
  <p><strong>EDIT 2:</strong> <a href="https://math.stackexchange.com/a/1093527/4471">This answer</a>, and in particular <a href="http://people.math.umass.edu/~hajir/m499c/vardy-complexity.pdf" rel="nofollow noreferrer">this paper</a> give a proof that this particular problem is NP-hard (Section 3), which hints that we shouldn't be looking for a polynomial algorithm. So the question becomes: "What is the best exponent we can get?".</p>