We get the same value for $p'=1-\delta-p$ instead of $p$, i.e. the graph is symmetric with respect to $p=\frac{1-\delta}{2}$. What would be sufficient to prove is that there are no more local minimizers. I guess we have convexity.

There do exist, for example $(i,j)=(7,1)$ and $(i,j)=(6,2)$

I agree with you and tried to simplify the matrix as much as possible. (However without permutations.)

I parametrized the matrix as $U*V'$ with $U,V\in \mathbb{R}^{n\times r}$ and solved the corresponding system of nonlinear equations. When I had a numerical solution I tried to find integer solutions with the same 0-1 patterns. For this I filled $r$ columns randomly with integers $(-3,..,3)$ and tested if it can be completed to a solution. After a lot of repetions the example above appeared.