Given a matrix of integers, $ \bf C $, find row vector $ \bf x $ of binaries and column vector $ \bf y $ of binaries that maximize the product $ \bf xCy $.
Is there an efficient algorithm to either find a solution or approximate solution? I've been researching (Integer) Linear Programming and can't seem to find an algorithm for this particular variant.
Here's an example. You're given matrix $$ \textbf{C} = \begin{bmatrix} -2 & 4 & -4\\ 8 & 4 & -7\\ -5&3&1\\ 2&2&2 \end{bmatrix} $$.
The solution of binary row/column vectors is:
$$ \textbf{x} = \begin{bmatrix} 1 & 1 & 0 & 1 \end{bmatrix}, \textbf{y} = \begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix} $$
where $ \textbf{xCy} = 18$
No, there is (likely) no efficient algorithm. The problem is NP-hard. You can prove it by reduction from MAX-CUT.
Here is the reduction. Given a graph $G$, let $A$ be the adjacency matrix, and let $C=A- c \cdot \text{Id}$, where $c$ is a large constant (larger than $|V|^2$). Consider the vectors $x,y$ that maximize $x^\top Cy$. It is easy to see that there is no index $i$ with $x_i=y_i=1$ (since that would incur a penalty of $-c$ to the objective function). It is also easy to see that the maximum is attained at $x,y$ such that $x+y=(1,1,\dots,1)$ (if $x_i+y_i<1$ for any $i$, then you can bump up either $x_i$ or $y_i$ to 1, without decreasing the value of the objective function). Therefore, each such $x,y$ corresponds to a cut, and vice versa. Now the value of the objective function is exactly the value of the cut, which proves that any algorithm for solving your problem could also be used to solve MAX-CUT.
You could use any standard method for dealing with NP-hardness: hope that your $C$ falls into a special case where MAX-CUT is easy; using an approximation algorithm; or use algorithms whose running time is exponential in the worst case.
I suspect that the most pragmatic solution will likely be to formulate this as an instance of some standard problem and use an off-the-shelf solver. One option is to recognize it as an instance of MIQP and use an off-the-shelf MIQP solver, as V.S.e.H. suggests. Alternatively, you could formulate it as an instance of ILP, and use an off-the-shelf solver. Specifically, introduce binary variables $t_{ij}$, with the intended meaning that $t_{ij}=x_i \land y_j$. Then the objective function is $\sum_{i,j} C_{ij} t_{ij}$. Also you can enforce validity of $t_{ij}$ by the method in https://cs.stackexchange.com/q/12102/755 for encoding logical-AND. If all entries in $C$ are non-negative, it suffices to constrain $t_{ij} \le x_i$ and $t_{ij} \le x_j$, which is simpler.