1) I have the following problem that I would like to first solve optimally but I have not been able to express it in a way that can be accepted by Matlab optimization functions.
$$\begin{array}{ll} \text{maximize} & {\bf c}^T{\bf x} \\ \text{subject to: } & {\bf A}_1{\bf x} = {\bf q} \\ & [{\bf B}_1 {\bf x}] \circ [{\bf B}_2 {\bf x}] = {\bf 0}\\ & [{\bf C}_1 {\bf x}] \circ [{\bf C}_2 {\bf x}] = {\bf 0}\\ & [{\bf D}_1 {\bf x}] \circ [{\bf D}_2 {\bf x}] = {\bf 0}\\ & x_i = \{0, 1\} \end{array}$$
The symbol $\circ$ represents the Hadamard product. The elements of ${\bf c}$ are all non-negative. The elements of matrix ${\bf A}$ are all non-negative. The matrices ${\bf B}_1, {\bf B}_2, {\bf C}_1$ and ${\bf C}_2$ are sparse matrices with either $0 $ or $1$ as elements. Any idea about how to transform this system to be a valid input to any optimizer?
2) Then, I have decided to relax the constraints involving a Hadamard product by multiplying ${\bf 1}^T$ at each side of the equality (unweighted sum of constraints). Thus, I obtain
$$\begin{array}{ll} \text{maximize} & {\bf c}^T{\bf x} \\ \text{subject to: } & {\bf A}_1{\bf x} = {\bf q} \\ & {\bf x}^T {\bf A}_2 {\bf x} = 0 \\ & {\bf x}^T {\bf A}_3 {\bf x} = 0 \\ & {\bf x}^T {\bf A}_4 {\bf x} = 0 \\ & x_i = \{0, 1\} \end{array}$$
where ${\bf A}_2, {\bf A}_3$ and ${\bf A}_4$ are also sparse matrices having $0$ or $1$ as their elements. Since these matrices are neither positive semidefinite or negative semidefinite, the relaxed problem is still non-convex (to the best of my knowledge). Can you please recommend a way to work this problem out a bit more? Is there any framework that allows me to further manipulate the problem and perhaps reach a simplified solution by hand?
Quadratic Forms
I believe these equations can be linearized by exploiting that $x$ are binary variables. Let's do the general equation $x^T A x= b$. This can be written as
$$ \sum_{i,j} x_i x_j a_{i,j} = b $$
Let's introduce a binary variable $y_{i,j} = x_i x_j$. This product can be linearized as:
$$\begin{align} &y_{i,j} \le x_i\\ &y_{i,j} \le x_j\\ &y_{i,j} \ge x_i+x_j-1 \end{align}$$
So now we can write:
$$ \sum_{i,j} y_{i,j} a_{i,j} = b $$
This will now give you a linear MIP model, which you can solve with
intlinprog.Now for some bonus points. If we want to be frugal we can create and exploit some symmetry as follows:
$$\begin{align} &y_{i,j} \le x_i&&\forall i<j\\ &y_{i,j} \le x_j&&\forall i<j\\ &y_{i,j} \ge x_i+x_j-1&&\forall i<j\\ &\sum_i x_i a_{i,i} + \sum_{i<j} y_{i,j} (a_{i,j}+a_{j,i}) = b \end{align}$$
This saves us a bunch of variables and equations.
Hadamard Products
(Updated) Make them quadratic forms and see above.
Or if you adventurous, there is another formulation possible. I am a bit less sure about this, but I think we have $(Ax) \circ (Bx) = 0$ or:
$$\begin{align} &u_i = \sum_j a_{i,j} x_j\\ &v_i = \sum_j b_{i,j} x_j\\ &u_i v_i = 0 \end{align}$$
Only the last equation is nonlinear. We can write this as $u_i=0 \text{ or } v_i=0$. Not so nice, but this can be linearized by:
$$\begin{align} &-M\delta_i \le u_i \le M \delta_i\\ &-M(1-\delta_i) \le v_i \le M (1-\delta_i)\\ &\delta_i \in \{0,1\} \end{align}$$
where $M$ is a bound on $u,v$. We may be able to find good bounds by carefully looking at the matrices $A,B$ (or even by running a series of small optimization problems).