Solving system of polynomial matrix equations over $\Bbb Z_2$

Question

Let $A, B, C, D$ be $4 \times 4$ matrices over $\mathbb{Z}_2$.

Suppose they satisfy \begin{equation} \begin{cases} A^2+ BC+ BCA+ ABC+ A= I_4\\ AB+ BCB+ ABD=0_{4 \times 4} \\ CA+ DCA+ CBC= 0_{4 \times 4} \\ DCB+ CBD= I_4 \end{cases}\,. \end{equation}

Does anyone have an idea on how to find matrices $A, B, C, D$ either theoretically or numerically? I'm thinking that we can find a solution since there are four equations and four unknowns.

Hints would suffice. Thank you so much.

Edit: I found out that the 148 solutions is not yet the total number of complete solutions. To possibly decrease the number of solutions I added another constraint:

$$\sum_{k=1}^4\sum_{l=1}^4 A_{1,k}A_{k,l}A_{l,1}+ \sum_{k=1}^4\sum_{l=1}^4 B_{1,k}C_{k,l}A_{l,1}+ \sum_{k=1}^4\sum_{l=1}^4 A_{1,k}B_{k,l}C_{l,1}+ \sum_{k=1}^4\sum_{l=1}^4 B_{1,k}D_{k,l}C_{l,1}=2Y_{5,1,1}.$$

I added this constraint to the last part of the code:

proc optmodel;
   num n = 4;
   set ROWS = 1..n;
   set COLS = ROWS;
   set CELLS = ROWS cross COLS;

   var A {CELLS} binary;
   var B {CELLS} binary;
   var C {CELLS} binary;
   var D {CELLS} binary;

   var Y {1..5, CELLS} >= 0 integer;
   
    /* M[i,j,k] = A[i,k]*A[k,j] */
   var M {1..n, 1..n, 1..n} binary;
   con Mcon1 {i in 1..n, j in 1..n, k in 1..n}:
      M[i,j,k] <= A[i,k];
   con Mcon2 {i in 1..n, j in 1..n, k in 1..n}:
      M[i,j,k] <= A[k,j];
   con Mcon3 {i in 1..n, j in 1..n, k in 1..n}:
      M[i,j,k] >= A[i,k] + A[k,j] - 1;
      
    /* X[i,j,k] = B[i,k]*C[k,j] */
   var X {1..n, 1..n, 1..n} binary;
   con Xcon1 {i in 1..n, j in 1..n, k in 1..n}:
      X[i,j,k] <= B[i,k];
   con Xcon2 {i in 1..n, j in 1..n, k in 1..n}:
      X[i,j,k] <= C[k,j];
   con Xcon3 {i in 1..n, j in 1..n, k in 1..n}:
      X[i,j,k] >= B[i,k] + C[k,j] - 1;
      
     /* O[i,j,k,l] = B[i,k]*C[k,l]*A[l,j] */
   var O {1..n, 1..n, 1..n, 1..n} binary;
   con Ocon1 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      O[i,j,k,l] <= B[i,k];
   con Ocon2 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      O[i,j,k,l] <= C[k,l];
   con Ocon3 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      O[i,j,k,l] <= A[l,j];   
   con Ocon4 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      O[i,j,k,l] >= B[i,k]+C[k,l]+A[l,j] - 2;  
      
      /* P[i,j,k,l] = A[i,k]*B[k,l]*C[l,j] */
   var P {1..n, 1..n, 1..n, 1..n} binary;
   con Pcon1 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      P[i,j,k,l] <= A[i,k];
   con Pcon2 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      P[i,j,k,l] <= B[k,l];
   con Pcon3 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      P[i,j,k,l] <= C[l,j];   
   con Pcon4 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      P[i,j,k,l] >= A[i,k]+B[k,l]+C[l,j] - 2;

   /* A^2 + BC + BCA + ABC + A = I_4 */
   con Con1 {<i,j> in CELLS}:
      sum {k in 1..n} M[i,j,k] 
    + sum {k in 1..n} X[i,j,k] 
    + sum {k in 1..n, l in 1..n} O[i,j,k,l] 
    + sum {k in 1..n, l in 1..n} P[i,j,k,l]
    + A[i,j] 
   = 2*Y[1,i,j] + (i=j);
   
   /* E[i,j,k] = A[i,k]*B[k,j] */
   var E {1..n, 1..n, 1..n} binary;
   con Econ1 {i in 1..n, j in 1..n, k in 1..n}:
      E[i,j,k] <= A[i,k];
   con Econ2 {i in 1..n, j in 1..n, k in 1..n}:
      E[i,j,k] <= B[k,j];
   con Econ3 {i in 1..n, j in 1..n, k in 1..n}:
      E[i,j,k] >= A[i,k] + B[k,j] - 1;
      
    /* F[i,j,k,l] = B[i,k]*C[k,l]*B[l,j] */
   var F {1..n, 1..n, 1..n, 1..n} binary;
   con Fcon1 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      F[i,j,k,l] <= B[i,k];
   con Fcon2 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      F[i,j,k,l] <= C[k,l];
   con Fcon3 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      F[i,j,k,l] <= B[l,j];   
   con Fcon4 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      F[i,j,k,l] >= B[i,k]+C[k,l]+B[l,j] - 2;
      
      /* G[i,j,k,l] = A[i,k]*B[k,l]*D[l,j] */
   var G {1..n, 1..n, 1..n, 1..n} binary;
   con Gcon1 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      G[i,j,k,l] <= A[i,k];
   con Gcon2 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      G[i,j,k,l] <= B[k,l];
   con Gcon3 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      G[i,j,k,l] <= D[l,j];   
   con Gcon4 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      G[i,j,k,l] >= A[i,k]+B[k,l]+D[l,j] - 2; 

   /* AB + BCB + ABD = 0 */
   con Con2 {<i,j> in CELLS}:
      sum {k in 1..n} E[i,j,k] 
    + sum {k in 1..n, l in 1..n} F[i,j,k,l] 
    + sum {k in 1..n, l in 1..n} G[i,j,k,l]
   = 2*Y[2,i,j];
   
   /* H[i,j,k] = C[i,k]*A[k,j] */
   var H {1..n, 1..n, 1..n} binary;
   con Hcon1 {i in 1..n, j in 1..n, k in 1..n}:
      H[i,j,k] <= C[i,k];
   con Hcon2 {i in 1..n, j in 1..n, k in 1..n}:
      H[i,j,k] <= A[k,j];
   con Hcon3 {i in 1..n, j in 1..n, k in 1..n}:
      H[i,j,k] >= C[i,k] + A[k,j] - 1;
      
    /* Q[i,j,k,l] = D[i,k]*C[k,l]*A[l,j] */
   var Q {1..n, 1..n, 1..n, 1..n} binary;
   con Qcon1 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      Q[i,j,k,l] <= D[i,k];
   con Qcon2 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      Q[i,j,k,l] <= C[k,l];
   con Qcon3 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      Q[i,j,k,l] <= A[l,j];   
   con Qcon4 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      Q[i,j,k,l] >= D[i,k]+C[k,l]+A[l,j] - 2; 
      
    /* R[i,j,k,l] = C[i,k]*B[k,l]*C[l,j] */
   var R {1..n, 1..n, 1..n, 1..n} binary;
   con Rcon1 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      R[i,j,k,l] <= C[i,k];
   con Rcon2 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      R[i,j,k,l] <= B[k,l];
   con Rcon3 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      R[i,j,k,l] <= C[l,j];   
   con Rcon4 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      R[i,j,k,l] >= C[i,k]+B[k,l]+C[l,j] - 2; 

   /* CA + DCA + CBC = 0 */
   con Con3 {<i,j> in CELLS}:
      sum {k in 1..n} H[i,j,k]
    + sum {k in 1..n, l in 1..n} Q[i,j,k,l] 
    + sum {k in 1..n, l in 1..n} R[i,j,k,l]
   = 2*Y[3,i,j];
   
   /* W[i,j,k,l] = D[i,k]*C[k,l]*B[l,j] */
   var W {1..n, 1..n, 1..n, 1..n} binary;
   con Wcon1 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      W[i,j,k,l] <= D[i,k];
   con Wcon2 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      W[i,j,k,l] <= C[k,l];
   con Wcon3 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      W[i,j,k,l] <= B[l,j];   
   con Wcon4 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      W[i,j,k,l] >= D[i,k]+C[k,l]+B[l,j] - 2; 
      
      /* T[i,j,k,l] = C[i,k]*B[k,l]*D[l,j] */
   var T {1..n, 1..n, 1..n, 1..n} binary;
   con Tcon1 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      T[i,j,k,l] <= C[i,k];
   con Tcon2 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      T[i,j,k,l] <= B[k,l];
   con Tcon3 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      T[i,j,k,l] <= D[l,j];   
   con Tcon4 {i in 1..n, j in 1..n, k in 1..n, l in 1..n}:
      T[i,j,k,l] >= C[i,k]+B[k,l]+D[l,j] - 2;

   /* DCB + CBD = I_4 */
   con Con4 {<i,j> in CELLS}:
      sum {k in 1..n, l in 1..n} W[i,j,k,l] 
    + sum {k in 1..n, l in 1..n} T[i,j,k,l]
   = 2*Y[4,i,j] + (i=j);
   
   /* U[1,1,k,l] = A[1,k]*A[k,l]*A[l,1] */
   var U {1..n, 1..n} binary;
   con Ucon1 {k in 1..n, l in 1..n}:
      U[k,l] <= A[1,k];
   con Ucon2 {k in 1..n, l in 1..n}:
      U[k,l] <= A[k,l];
   con Ucon3 {k in 1..n, l in 1..n}:
      U[k,l] <= A[l,1];   
   con Ucon4 {k in 1..n, l in 1..n}:
      U[k,l] >= A[1,k]+A[k,l]+A[l,1] - 2;
      
    /* V[1,1,k,l] = B[1,k]*D[k,l]*C[l,1] */
   var V {1..n, 1..n} binary;
   con Vcon1 {k in 1..n, l in 1..n}:
      V[k,l] <= B[1,k];
   con Vcon2 {k in 1..n, l in 1..n}:
      V[k,l] <= D[k,l];
   con Vcon3 {k in 1..n, l in 1..n}:
      V[k,l] <= C[l,1];   
   con Vcon4 {k in 1..n, l in 1..n}:
      V[k,l] >= B[1,k]+D[k,l]+C[l,1] - 2;  
   
   con 
   sum {k in 1..n, l in 1..n} U[k,l] 
    + sum {k in 1..n, l in 1..n} V[k,l]+ sum {k in 1..n, l in 1..n} O[1,1,k,l] + sum {k in 1..n, l in 1..n} P[1,1,k,l] 
   = 2*Y[5,1,1] ;
   

   solve noobj with milp / maxpoolsols=100 soltype=best;
   create data out(drop=s) from [s]=(1.._NSOL_) {i in 1..n, j in 1..n} <col('A_'||i||''||j)=A[i,j].sol[s]> {i in 1..n, j in 1..n} <col('B'||i||''||j)=B[i,j].sol[s]> {i in 1..n, j in 1..n} <col('C'||i||''||j)=C[i,j].sol[s]> {i in 1..n, j in 1..n} <col('D'||i||'_'||j)=D[i,j].sol[s]>;  
   
   print A;
   print B;
   print C;
   print D;
quit;

The software found 100 solutions, but when I checked the output data, it is the same solution replicated 100 times. Is there something wrong with the syntax of the constraint that I added to the code?

Answer 1

You can formulate this as an integer linear programming problem, and here is a feasible solution: \begin{align} A&= \begin{pmatrix} 0 &0 &0 &1 \\ 1 &0 &1 &0 \\ 1 &1 &1 &1 \\ 1 &0 &0 &1 \\ \end{pmatrix} \\ B&= \begin{pmatrix} 1 &1 &0 &0 \\ 0 &1 &1 &0 \\ 1 &0 &1 &1 \\ 0 &0 &1 &0 \\ \end{pmatrix} \\ C&= \begin{pmatrix} 1 &1 &0 &1 \\ 1 &1 &0 &0 \\ 1 &0 &0 &0 \\ 1 &0 &1 &0 \\ \end{pmatrix} \\ D&= \begin{pmatrix} 0 &0 &1 &1 \\ 0 &0 &0 &1 \\ 0 &0 &0 &0 \\ 0 &0 &0 &0 \\ \end{pmatrix} \end{align} There are 64 constraints, 64 binary decision variables $(A_{i,j},B_{i,j},C_{i,j},D_{i,j})$, and 64 nonnegative integer variables ($Y_{t,i,j}$ for $t\in\{1,2,3,4\}$ and all $(i,j)$). Explicitly, the constraints are \begin{align} \sum_{k=1}^n \left(A_{i,k} A_{k,j} + B_{i,k} C_{k,j}\right) + \sum_{k=1}^n \sum_{\ell=1}^n \left(B_{i,k} C_{k,\ell} A_{\ell,j} + A_{i,k} B_{k,\ell} C_{\ell,j}\right) + A_{i,j} &= 2 Y_{1,i,j} + [i=j] \\ \sum_{k=1}^n A_{i,k} B_{k,j} + \sum_{k=1}^n \sum_{\ell=1}^n \left(B_{i,k} C_{k,\ell} B_{\ell,j} + A_{i,k} B_{k,\ell} D_{\ell,j}\right) &= 2 Y_{2,i,j} \\ \sum_{k=1}^n C_{i,k} A_{k,j} + \sum_{k=1}^n \sum_{\ell=1}^n \left(D_{i,k} C_{k,\ell} A_{\ell,j} + C_{i,k} B_{k,\ell} C_{\ell,j}\right) &= 2 Y_{3,i,j} \\ \sum_{k=1}^n \sum_{\ell=1}^n \left(D_{i,k} C_{k,\ell} B_{\ell,j} + C_{i,k} B_{k,\ell} D_{\ell,j}\right) &= 2 Y_{4,i,j} + [i=j] \end{align} These can be linearized by replacing each product with a new binary variable, together with additional linear constraints, as shown in https://or.stackexchange.com/questions/37/how-to-linearize-the-product-of-two-binary-variables

---

As requested, here is the SAS code. The NOOBJ option specifies that there is no objective, and the LINEARIZE option automatically performs the linearization.

proc optmodel;
   num n = 4;
   set ROWS = 1..n;
   set COLS = ROWS;
   set CELLS = ROWS cross COLS;

   var A {CELLS} binary;
   var B {CELLS} binary;
   var C {CELLS} binary;
   var D {CELLS} binary;

   var Y {1..4, CELLS} >= 0 integer;

   /* A^2 + BC + BCA + ABC + A = I_4 */
   con Con1 {<i,j> in CELLS}:
      sum {k in 1..n} A[i,k]*A[k,j] 
    + sum {k in 1..n} B[i,k]*C[k,j] 
    + sum {k in 1..n, l in 1..n} B[i,k]*C[k,l]*A[l,j] 
    + sum {k in 1..n, l in 1..n} A[i,k]*B[k,l]*C[l,j]
    + A[i,j] 
   = 2*Y[1,i,j] + (i=j);

   /* AB + BCB + ABD = 0 */
   con Con2 {<i,j> in CELLS}:
      sum {k in 1..n} A[i,k]*B[k,j] 
    + sum {k in 1..n, l in 1..n} B[i,k]*C[k,l]*B[l,j] 
    + sum {k in 1..n, l in 1..n} A[i,k]*B[k,l]*D[l,j]
   = 2*Y[2,i,j];

   /* CA + DCA + CBC = 0 */
   con Con3 {<i,j> in CELLS}:
      sum {k in 1..n} C[i,k]*A[k,j] 
    + sum {k in 1..n, l in 1..n} D[i,k]*C[k,l]*A[l,j] 
    + sum {k in 1..n, l in 1..n} C[i,k]*B[k,l]*C[l,j]
   = 2*Y[3,i,j];

   /* DCB + CBD = I_4 */
   con Con4 {<i,j> in CELLS}:
      sum {k in 1..n, l in 1..n} D[i,k]*C[k,l]*B[l,j] 
    + sum {k in 1..n, l in 1..n} C[i,k]*B[k,l]*D[l,j]
   = 2*Y[4,i,j] + (i=j);

   solve noobj linearize;

   print A;
   print B;
   print C;
   print D;
quit;

Another approach to this feasibility problem is to introduce nonnegative surplus and slack variables and minimize their use. The original problem is feasible if and only if all surplus and slack variables are 0.

   var Surplus {1..4, CELLS} >= 0;
   var Slack   {1..4, CELLS} >= 0;

   min Infeasibility = sum {t in 1..4, <i,j> in CELLS} (Surplus[t,i,j] + Slack[t,i,j]);

   /* A^2 + BC + BCA + ABC + A = I_4 */
   con Con1 {<i,j> in CELLS}:
      sum {k in 1..n} A[i,k]*A[k,j] 
    + sum {k in 1..n} B[i,k]*C[k,j] 
    + sum {k in 1..n, l in 1..n} B[i,k]*C[k,l]*A[l,j] 
    + sum {k in 1..n, l in 1..n} A[i,k]*B[k,l]*C[l,j]
    + A[i,j] 
    - Surplus[1,i,j] + Slack[1,i,j]
   = 2*Y[1,i,j] + (i=j);

   /* AB + BCB + ABD = 0 */
   con Con2 {<i,j> in CELLS}:
      sum {k in 1..n} A[i,k]*B[k,j] 
    + sum {k in 1..n, l in 1..n} B[i,k]*C[k,l]*B[l,j] 
    + sum {k in 1..n, l in 1..n} A[i,k]*B[k,l]*D[l,j]
    - Surplus[2,i,j] + Slack[2,i,j]
   = 2*Y[2,i,j];

   /* CA + DCA + CBC = 0 */
   con Con3 {<i,j> in CELLS}:
      sum {k in 1..n} C[i,k]*A[k,j] 
    + sum {k in 1..n, l in 1..n} D[i,k]*C[k,l]*A[l,j] 
    + sum {k in 1..n, l in 1..n} C[i,k]*B[k,l]*C[l,j]
    - Surplus[3,i,j] + Slack[3,i,j]
   = 2*Y[3,i,j];

   /* DCB + CBD = I_4 */
   con Con4 {<i,j> in CELLS}:
      sum {k in 1..n, l in 1..n} D[i,k]*C[k,l]*B[l,j] 
    + sum {k in 1..n, l in 1..n} C[i,k]*B[k,l]*D[l,j]
    - Surplus[4,i,j] + Slack[4,i,j]
   = 2*Y[4,i,j] + (i=j);

   solve linearize;

To find more than one solution, you can use the MAXPOOLSOLS= option.

---

Here is an example of how you can manually linearize the first product:

   /* M[i,j,k] = A[i,k]*A[k,j] */
   var M {1..n, 1..n, 1..n} binary;
   con Mcon1 {i in 1..n, j in 1..n, k in 1..n}:
      M[i,j,k] <= A[i,k];
   con Mcon2 {i in 1..n, j in 1..n, k in 1..n}:
      M[i,j,k] <= A[k,j];
   con Mcon3 {i in 1..n, j in 1..n, k in 1..n}:
      M[i,j,k] >= A[i,k] + A[k,j] - 1;

You would then replace A[i,k]*A[k,j] with M[i,j,k] throughout the original nonlinear constraints.

You can optionally declare the new variables as >= 0 rather than binary.

---

To exclude a given set of solutions, first read them in:

set SOLUTIONS;
num Asol {CELLS, SOLUTIONS};
num Bsol {CELLS, SOLUTIONS};
num Csol {CELLS, SOLUTIONS};
num Dsol {CELLS, SOLUTIONS};
read data out into SOLUTIONS=[_N_] 
  {i in 1..n, j in 1..n} <Asol[i,j,_N_]=col('A_'||i||''||j)> 
  {i in 1..n, j in 1..n} <Bsol[i,j,_N_]=col('B_'||i||''||j)>
  {i in 1..n, j in 1..n} <Csol[i,j,_N_]=col('C_'||i||''||j)>
  {i in 1..n, j in 1..n} <Dsol[i,j,_N_]=col('D_'||i||''||j)>;

Then impose a no-good constraint:

con NoGood {sol in SOLUTIONS}:
  sum {<i,j> in CELLS} (
     (if Asol[i,j,sol] = 0 then A[i,j] else 1 - A[i,j])
   + (if Bsol[i,j,sol] = 0 then B[i,j] else 1 - B[i,j])
   + (if Csol[i,j,sol] = 0 then C[i,j] else 1 - C[i,j])
   + (if Dsol[i,j,sol] = 0 then D[i,j] else 1 - D[i,j])
  ) >= 1;

Then call the solver:

solve noobj;

If you do this for the 148 solutions you found and the resulting solve comes back infeasible, then you have found all solutions.

Answer 2

I have made an algebraic attempt that might help you somewhere. Lets start with the second and third equation: $$AB(I+D) = -BCB \rightarrow AB = -BCB(I+D)^{-1}$$ $$(I+D)CA = -CBC \rightarrow CA = -(I+D)^{-1}CBC$$ Now plug these results in the first equation: $$A^2 + BC - B(I+D)^{-1}CBC - BCB(I+D)^{-1}C+A = I$$ From this, two things can be observed: $$A = -B(I+D)^{-1}C, ~~ \text{and} ~~ BCA = ABC$$ Next, lets postmultiply the second equation with $C$ and the premultiply the third with $B$: $$ABC + BCBC + ABDC = 0, ~~ BCA + BDCA + BCBC = 0$$ $$BCA + BDCA + BCBC = ABC + BCBC + ABDC \rightarrow BDCA = ABDC$$ Now if we substitute the definition for $A$ in this result the following can also be deduced: $$BDCA = ABDC \rightarrow -BDCB(I+D)^{-1}C = -B(I+D)^{-1}CBDC$$ $$DCB(I+D)^{-1}=(I+D)^{-1}CBD$$ All these commutations are starting me to believe that $A$ is something like $\alpha I$, and therefore $-(I+D)^{-1} = \alpha (CB)^{-1}$. Hope this does help you get to an answer.

Answer 3

Here are two simplifications:

1. Up to similarity, there are only six choices of $D$. Note that the fourth equation is a Sylvester equation of the form $DX+XD=I$. By considering the Jordan form of $D$ over the algebraic closure of $\mathbb Z_2=GF(2)$, we see that $DX+XD=I$ is solvable only if all Jordan blocks of $D$ have even sizes. It follows that the characteristic polynomial $p(x)$ of $D$ can always be factorised as $$ p(x)=(x+\lambda)^2(x+\mu)^2 =(x^2+\lambda^2)(x^2+\mu^2) =x^4+(\lambda^2+\mu^2)x^2+\lambda^2\mu^2.\tag{$\ast$} $$ Since $p$ has coefficients in $GF(2)$, there are only four possible choices of $p$, namely $x^4,\,x^2(x+1)^2,\,(x+1)^4$ or $x^4+x^2+1$. In other words, up to similarity, we can always assume that $D$ is equal to one of the five Jordan forms $$ J_4(0),\,\pmatrix{J_2(0)\\ &J_2(0)},\,\pmatrix{J_2(0)\\ &J_2(1)}\,\pmatrix{J_2(1)\\ &J_2(1)}\, J_4(1) $$ (where $J_m(\lambda)$ denotes the $m\times m$ Jordan block with eigenvalue $\lambda$) or the companion matrix $$ \pmatrix{0&0&0&1\\ 1&0&0&0\\ 0&1&0&1\\ 0&0&1&0}. $$
2. You may replace $D$ by $D+I$. With the new $D$, the first and the forth equations then remain unchanged, but the second and the third equations become \begin{cases} BCB+ABD=0,\\ DCA+CBC=0,\\ \end{cases} which involve fewer summands.

There are also three observations of which I don't know the usefulness. First, (with the old $D$ or the new $D$) the second and the third equations together give the Sylvester equation $D(CAB)+(CAB)D=0$. Second, the second or the third equation implies that at least one of $B$ or $C$ is necessarily singular when one of $A$ or $D$ is singular. Finally, the LHS of the first equation can be rewritten as $$ (I+A)A+A(I+BC)+(I+BC)(I+A)=0, $$ where the LHS is a cyclic sum.