Is there an algorithm in $\mathbf{P}$ that solves this overlapping-XOR-grid problem? Suppose you have an $m \times n$ array which can hold boolean values. It is initially empty, full of uniform $0$s.
The only way to change the array is by supplying two pairs of coordinates, defining a rectangular area within the array; every bit in that area is then flipped. A first such move would give you something like this:

(I'm not sure of the ideal notation for that, but it's (2,4)-(8,15) by my reckoning, defining opposite inclusive corner coordinates.)
However, observe what happens when you play a second move, (4,4)-(6,17):

Being an exclusive-or, the middle bits flip back to $0$ while the overhang on the right flips to $1$. Naively, one could look at this second picture and conclude that three moves were necessary to reach that configuration if starting from all zeros, one move for each rectangle. As we've seen, this was done with only two.
This is a straightforward example, but when you combine even a handful of overlapping moves together, it appears to become extremely difficult to disentangle. If I'm not mistaken, this is not the case with an analogous one-dimensional array, where it's simple to derive the optimal series of moves (that is, XORing line segments) to recreate a configuration in $O(n)$ time. This is because every unequal adjacent pair implies an XOR segment with an endpoint between that pair.
The goal is to calculate the minimum number of moves needed to recreate any given final configuration, starting with an empty array. I'll consider my question answered if someone can tell me whether an efficient algorithm is known for this or something similar, can sketch one themselves, or can provide an informed opinion on why such an algorithm is or is not likely to exist. By efficient algorithm, I mean one that runs in polynomial time in the general case.
(I realize this is arguably a CS question, but I'm trying here first.)

Confirmed RobPratt's nice solution.
For future reference, this Mathematica code will run a test with random parameters and give you output as in the screenshot below.
a[i_, j_, i1_, j1_, i2_, j2_] := Boole[i1 <= i <= i2 && j1 <= j <= j2]
x[coords_] := Boole@MemberQ[rects, coords]
m = Table[0, 9, 13];
randrect := Module[{x1, y1, x2, y2},
   {x1, x2} = Sort[RandomInteger[{1, Length@First@m}, 2]];
   {y1, y2} = Sort[RandomInteger[{1, Length@m}, 2]];
   {y1, x1, y2, x2}
   ];
rects = Table[randrect, RandomInteger[{1, 13}]]
Do[
  m[[rect[[1]] ;; rect[[3]], rect[[2]] ;; rect[[4]]]] = 
    1 - m[[rect[[1]] ;; rect[[3]], rect[[2]] ;; rect[[4]]]];
  , {rect, rects}];
ArrayPlot[m, Mesh -> True]
t = Minimize[{Sum[
     x[{i1, j1, i2, j2}], {i2, Length@m}, {j2, Length@First@m}, {i1, 
      i2}, {j1, j2}], 
    And @@ Flatten@
      Table[y[i, j] \[Element] Integers && y[i, j] >= 0 && 
        Sum[a[i, j, i1, j1, i2, j2] x[{i1, j1, i2, j2}], {i2, 
           Length@m}, {j2, Length@First@m}, {i1, i2}, {j1, j2}] == 
         2 y[i, j] + m[[i, j]], {i, Length@m}, {j, Length@First@m}]}, 
   Flatten[Table[y[i, j], {i, Length@m}, {j, Length@First@m}]]];
TextGrid[{{Length@rects, "rectangles drawn"}, {First[t], 
   "rectangles calculated"}}]
ArrayPlot[Partition[Last /@ Last@t, Length@First@m], Mesh -> True]

Example output:

 A: As in the Lights Out game, you can solve the problem via Gaussian elimination mod 2, which is $O(k^3)$ for a $k \times k$ matrix.  See How can I solve a vector equation in Z2?

Suppose the given $m\times n$ matrix has entries $B_{i,j}\in\{0,1\}$.  Let binary decision variable $x_{i_1,j_1,i_2,j_2}$ indicate whether the rectangular submatrix with top-left corner $(i_1,j_1)$ and bottom-right corner $(i_2,j_2)$ is selected.  Let $a_{i,j,i_1,j_1,i_2,j_2}\in\{0,1\}$ indicate whether entry $(i,j)$ is in the rectangular matrix defined by $(i_1,j_1)$ and $(i_2,j_2)$.  With Iverson bracket notation,
$$a_{i,j,i_1,j_1,i_2,j_2}=[i_1 \le i \le i_2][j_1 \le j \le j_2].$$
The problem is to minimize $$\sum_{i_1,j_1,i_2,j_2} x_{i_1,j_1,i_2,j_2}$$ subject to the $mn$ equations
$$\sum_{i_1,j_1,i_2,j_2} a_{i,j,i_1,j_1,i_2,j_2} x_{i_1,j_1,i_2,j_2} \equiv B_{i,j} \pmod2 \quad \text{for all $i,j$} \tag1$$
You can reformulate $(1)$ with linear equations by introducing nonnegative integer decision variables $y_{i,j}$ and replacing
$$\equiv B_{i,j} \pmod2$$
with
$$= 2y_{i,j} + B_{i,j}$$
Then call an integer linear programming solver.
