I'm going to provide two algorithms here:

- A super-exponential-time solution that works in all cases.
- A polynomial-time solution that applies if the mean values are known and if a certain matrix is positive semi-definite.

The polynomial-time solution essentially recycles the end of the proof of Goemans-Williamson's result for approximating MAX-CUT. I include the super-exponential-time case just because (Dustin G. Mixon's clear sketch notwithstanding) the other solutions do not seem to specify their algorithms completely.

# General setup

Let $x=(x_1,...,x_n)$ be the random variable with covariance matrix $C$.

Note that by a "random Bernoulli vector" we might mean $x_i\in\{0,1\}$ or $x_i\in\{-1,1\}$. We can convert the covariance matrix from the former version to the latter by multiplying by 4, so we can adopt either convention without issue. I'll consider $\{-1,1\}$.

Let $\overline{x_i}=\mathbb{E}[x_i]$. We are given the covariances:
$$c_{i,j}=\mathbb{E}[(x_i-\overline{x_i})(x_j-\overline{x_j})]=\mathbb{E}[x_ix_j]-\overline{x_i}\,\overline{x_j}$$

Note that $c_{i,i}=\mathbb{E}[x_i^2]-\overline{x_i}^2=1-\overline{x_i}^2$, which implies that
$$\overline{x_i}=\pm\sqrt{1-c_{i,i}}$$

# Super-exponential-time algorithm

With an outer loop of size $2^n$, we can try all possible sign combinations for the $\overline{x_i}$.

For the inner loop, then, we can assume that we know $\overline{x_i}$. Suppose we are in the inner loop.

Let $\mathcal{B}$ be the domain of $x$ (so, $\mathcal{B}$ has $2^n$ elements). Consider a linear program with variables $p_B$ for each $B\in \mathcal{B}$. This represents the probability that $x=B$. All the constraints can be expressed as linear combinations of these variables.

To see how, let $B=(b_1,...,b_n)$. Then the constraint that we have a probability distribution is:
$$0\leq p_B \leq 1$$
$$\sum_{B\in\mathcal{B}} p_B=1$$
The value of $\overline{x_i}$ (given to us by the outer loop):
$$\left(\sum_{B\in\mathcal{B}, b_i=1}p_B\right)-\left(\sum_{B\in\mathcal{B}, b_i=-1}p_B\right)=\overline{x_i}$$
The covariances (for $i\neq j$):
$$\left(\sum_{B\in\mathcal{B}, b_i=b_j}p_B\right)-\left(\sum_{B\in\mathcal{B}, b_i\neq b_j}p_B\right)=c_{i,j}+\overline{x_i}\overline{x_j}$$

If there exists a distribution of Bernoulli vectors consistent with $C$, then the linear program will be feasible (which we can determine in time polynomial in $|\mathcal{B}|=2^n$), and we exit, returning the distribution. On the other hand,
if all the linear programs are infeasible, then $C$ is not consistent with any random variable over Bernoulli vectors

# Polynomial-time algorithm

We provide an algorithm that works under the assumptions that (1) the $\overline{x_i}$ are known, and (2) a certain matrix (specified below) is positive semi-definite.

Let $f(x)=\sin(\pi x/2)$, and note that $f:[-1,1]\rightarrow [-1,1]$. We will abuse notation and apply $f$ element-wise to matrices, i.e. for any matrix $X$ with elements in $[-1,1]$, we write $f(X)=(f(x_{i,j}))$.

Consider, the second moments:
$$d_{i,j}=\mathbb{E}[x_ix_j]$$
Let $D=(d_{i,j})$ be the matrix of second moments.

Note that if we are given first moments $\overline{x_i}$, we can translate between $C$ and $D$:
$$d_{i,j}-\overline{x_i}\,\overline{x_j}=c_{i,j}$$

Note that $d_{i,i}=1$, so that in some sense $C$ contains more information than $D$. To address that imbalance, consider the $(n+1)\times (n+1)$ matrix $G=(g_{i,j})$. For $1\leq i,j \leq n$, we set $g_{i,j}=d_{i,j}$. For $1\leq i \leq n$, we set $g_{i,n+1}=g_{n+1,i}=\overline{x_i}$. Finally, $g_{n+1,n+1}=1$.

Note that all the entries of $G$ lie within $[-1,1]$. Set
$$H=f(G)$$

If $H$ is a positive semidefinite matrix, then the following algorithm will produce a Bernoulli vector with covariance matrix $C$:

- Produce a sample $(a_1,...,a_{n+1})$ from a Gaussian random variable with covariance matrix $H$.
- Threshold the vector by setting $b_i=sign(a_i)$.
- Set $c_i=b_i b_{n+1}$ for $i=1,...,n$.
- Return $(c_1,...c_n)$.

Why does this work? Using a Cholesky decomposition, if $H$ is rank $r$, we can write $H=J^TJ$, where $J$ is an $r\times n$ matrix. Since the main diagonal of $H$ is all ones, each column of $J$ is a vector that lies on the unit sphere in $R^r$.

Consider two columns, $s,t$, lying on the unit sphere in $R^r$. Select a random $v$ (as an i.i.d. normal random vector). Let $\hat{s}=sign(v \cdot s)$, and $\hat{t}=sign(v \cdot t)$. In other words, if we put a hyperplane perpendicular to $v$ through the origin and split the sphere into two halves, $\hat{s}=1$ if it lies in the same half as $v$, and -1 otherwise.

Consider the two-dimensional plane spanned by $s$ and $t$. The vectors lie on the unit circle. The hyperplane divides the circle into two halves. If the angle between $s$ and $t$ is $\theta=\arccos(s \cdot t)$, then the covariance between $\hat{s}$ and $\hat{t}$ is (noting that $\mathbb{E}[\hat{s}]=\mathbb{E}[\hat{t}]=0)$:
$$\mathbb{E}[\hat{s}\hat{t}]$$
which is the probability that they lie on the same side of the circle $(1-\theta/\pi)$ minus the probability that they lie on opposite sides $(\theta/\pi)$, which is:
$$=1-\frac{2}{\pi}\theta$$
$$=\frac{2}{\pi}\arcsin(s \cdot t)=f^{-1}(s\cdot t)$$

Some steps of the algorithm may make more sense now:

- The $f()$ function tells us what Gaussian covariance is necessary to map to the target covariance on the Bernoulli random variables.
- With a little squinting, you can recognize that the business with $v$ and the hyperplane is the same as selecting a Gaussian random vector and thresholding it.

The algorithm above will produce $(b_1,...,b_{n+1})$ with the second moment matrix given by $G$. Unfortunately, $\mathbb{E}[b_i]=0$ for all $i$.
To recover our target mean values $\overline{x_i}$,
we use $b_{n+1}$. Letting $c_i=b_ib_{n+1}$, note that the covariance of $(c_i,c_j)$ is the same as $(b_i,b_j)$. However, the mean of $c_i$ has shifted to $\overline{x_i}$.

# Other notes

Because of paywalls, I haven't been able to follow V.C.'s references yet, but the presence of the arcsine above makes me suspicious that I'm reinventing those wheels. The definition of Van Vleck's arcsine law given here seems related but a little different from what we need. When I manage to track down more of V.C.'s references I'll update this post.

If we are given the second moment matrix $D$ directly, we do not need to make any assumptions about knowing the $\overline{x_i}$; we only use those values to convert from $C$ to $D$.

If the main diagonal of the covariance matrix is identically one, that implies that $\overline{x_i}=1$ or all $i$, so in that case we can deduce all the $\overline{x_i}$.

The extra row and column in $G$ is a bit of a hack; it seems like it should be possible to incorporate the mean more directly and perhaps apply to a wider class of covariance matrices.

Rather than guessing all $2^n$ possibilities for the signs of the $\overline{x}_i$, we could also express this problem as a semidefinite program with a rank 1 constraint, and remove the rank 1 constraint to make the problem tractable. However, I don't know any guarantees for how well that method would work.

I suspect that a polynomial time solution for all $C$ would probably solve MAX-CUT (and hence imply that $P=NP$). So such an algorithm is unlikely to exist.