This is a case of Power Group Enumeration with the group permuting
the slots being the eight symmetries $G_N$ of the $N\times N$ square
and the group acting on the $Q$ colors being the symmetric group
$S_Q$.
The cycle indices for $G_N$ were carefully documented and computed at
the following MSE link I.
The cycle index of the symmetric group can be computed from the
classical recurrence by Lovasz.
It then remains to apply the Power Group Enumeration formula /
algorithm as documented at the following
MSE link II.
We get for the case of coloring a square with at most two
interchangeable colors
$$1, 4, 51, 4324, 2105872, 4295327872, 35184441295872,
\\ 1152921514807410688,\ldots$$
which is [OEIS A182044](https://oeis.org/A182044)
where a closed formula can be found.
For at most three colors we get the sequence
$$1, 6, 490, 901012, 17653123147, 3126972103187548,
\\ 4985402694248850150928, 71535079589434063488162675274,
\\ 9238051838396620455005158025879826301,\ldots $$
Finally for at most four colors we get the sequence
$$1, 7, 1555, 22396971, 5864091026091, 24595658827938966187,
\\ 1650586719048786316922366635, 1772303994379887884341412962742479531,
\\ 30447950777727144129269702965544605972525918891,\ldots$$
The sequence of colorings using any number of available interchangeable colors is also quite interesting and starts
$$1, 7, 2966, 1310397193,
\\ 579823814813639193, 477464341236883456112705749206,
\\ 1340767144321669800049265230788088027597024910,
\\ 21516767919669856366796245458194341840929552797762722429430679631,
\ldots$$
An implementation of this algorithm is included below.
with(combinat);
pet_cycleind_symm :=
proc(n)
option remember;
local l;
if n=0 then return 1; fi;
expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;
pet_cycleind_grid :=
proc(N)
option remember;
local cind;
if type(N, even) then
cind :=
1/8*(a[1]^(N^2)+3*a[2]^(N^2/2)+
2*a[1]^N*a[2]^((N^2-N)/2) + 2*a[4]^(N^2/4));
else
cind :=
1/8*(a[1]^(N^2)+4*a[1]^N*a[2]^((N^2-N)/2)+
a[1]*a[2]^((N^2-1)/2) + 2*a[1]*a[4]^((N^2-1)/4));
fi;
cind;
end;
gridcols :=
proc(N, Q)
option remember;
local idx_slots, idx_cols, res, term_a, term_b,
v_a, v_b, inst_a, inst_b, len_a, len_b, p, q;
if N > 1 then
idx_slots := pet_cycleind_grid(N);
else
idx_slots := [a[1]];
fi;
if Q > 1 then
idx_cols := pet_cycleind_symm(Q);
else
idx_cols := [a[1]];
fi;
res := 0;
for term_a in idx_slots do
for term_b in idx_cols do
p := 1;
for v_a in indets(term_a) do
len_a := op(1, v_a);
inst_a := degree(term_a, v_a);
q := 0;
for v_b in indets(term_b) do
len_b := op(1, v_b);
inst_b := degree(term_b, v_b);
if len_a mod len_b = 0 then
q := q + len_b*inst_b;
fi;
od;
p := p*q^inst_a;
od;
res := res +
lcoeff(term_a)*lcoeff(term_b)*p;
od;
od;
res;
end;
Let $V$ be the vector space over the field with two elements, $\Bbb F_2$, generated by all matrices of shape $N\times N$, $N=18$, which have exactly one row, or exactly one column with one entries, all other entries being zero. We consider the following map from $V$ to $\Bbb F_2$:
$X$ in $V$ is mapped to
$$
f(X)=
(x_{11}+x_{22}+\dots+x_{N-1,N-1}+x_{N,N})
+
(x_{12}+x_{23}+\dots+x_{N-1,N}+x_{N,1})\ .
$$
(Sum on two consecutive fixed diagonals of $X$, where the indices are considered modulo $N$.)
Then each generator of $V$ is mapped to zero, since the "bits" taken in $f(X)$ are exactly two in each row and/or column. So $f$ vanishes on $V$. But it does not vanish on any quadratic submatrix $A=A(I)$ with ones exactly on the positions $I\times I$, $I$ being a proper interval of $\{1,2,\dots,N\}$. (For instance $I=\{1,\dots,16\}$ for our $N=18$.)
So in our case the answer to the problem is negative.
(The "change of the color on a line/columns" corresponds to adding a generator of $V$ to a matrix. We are starting with a zero matrix. The matrices that can be obtained are exactly those in the vector space $V$ of dimension $N+N-1$.)
Note: I can provide some simple computer check in sage, if needed.
Later edit: It is simpler to test the equations
$$
x_{is}-x_{js}-x_{it}+x_{jt}=0\ ,\qquad 1\le i<j\le N\ ,\ 1\le s<t\le N\ ,
$$
that can be extracted from all $2\times 2$ minors of a matrix $X$ in order to test/decide if a matrix $X$ can be realized. (The minus is also a plus, but makes it simpler to see that e.g. each of the above "one row matrix of ones" satisfy the equation, explicitly $1-1-0+0=0$, and correspondingly for the column matrices in the base of $V$...)
Best Answer
Your analysis for the $8 \times 8$ square amounts to saying that each operation leaves the parity of the total number of black squares unchanged. This assertion is accurate, and therefore, the analysis is valid. The assertion critically depends on the idea that $(8-1)$ has the same parity as $(1)$.
For the $3 \times 3$ square, this specific approach is invalid, because (for example) $(3-1)$ does not have the same parity as $(1)$.
At this point, it is unclear whether the same conclusion is accurate in the $3 \times 3$ square. That is, there may exist a similar but different valid argument that involves $\pmod{n}$, where $n$ is some positive integer other than $(2)$.
For what it's worth, an alternative (convoluted and therefore inferior) approach for the $8 \times 8$ square is to recognize that the computation of $(W - B)$ (i.e. # of white squares - number of black squares) is unchanged $\pmod{4}$ because (for example) $[(+6) - (-6)] \equiv 0\pmod{4}$.
This yields a similar contradiction, as your analysis of the $8 \times 8$ square, because $(63 - 1) \equiv 2\pmod{4}.$