Generate Bernoulli vector with given covariance matrix

Question

I am from different background, so please forgive me if the answer is so well known.

Let $C=(c_{ij})$ be a given $n\times n$ matrix. Do we have a way to generate samples of random Bernoulli vectors with covariance matrix equal to $C$?

More specifically, are there functions that take in Bernoulli random vectors with i.i.d. entries, and output Bernoulli random vectors with covariance matrix equaling $C$.

Answer 1

I'm going to provide two algorithms here:

1. A super-exponential-time solution that works in all cases.
2. 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$:

1. Produce a sample $(a_1,...,a_{n+1})$ from a Gaussian random variable with covariance matrix $H$.
2. Threshold the vector by setting $b_i=sign(a_i)$.
3. Set $c_i=b_i b_{n+1}$ for $i=1,...,n$.
4. 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.

Answer 2

This post will essentially be a list of literature pointers, and as such does not have much mathematical content as bibliographical. This problem has drawn some interest in the past; the best method I am aware of for this task is the "arcsine law", which has a very long story, starting from the 1966 article:

Van Vleck, J. H., & Middleton, D. (1966). The spectrum of clipped noise. Proceedings of the IEEE, 54(1), 2-19.

and, e.g., reused in the synthesis of binary textures nearly 20 years ago:

Jacovitti, G., Neri, A., & Scarano, G. (1998). Texture synthesis-by-analysis with hard-limited Gaussian processes. Image Processing, IEEE Transactions on, 7(11), 1615-1621.

I can provide you at least two more recent references from an author I know well:

Caprara, A., Furini, F., Lodi, A., Mangia, M., Rovatti, R., & Setti, G. (2014). Generation of antipodal random vectors with prescribed non-stationary 2-nd order statistics. Signal Processing, IEEE Transactions on, 62(6), 1603-1612.

Rovatti, R., Mazzini, G., Setti, G., & Vitali, S. (2008, May). Linear probability feedback processes. In Circuits and Systems, 2008. ISCAS 2008. IEEE International Symposium on (pp. 548-551). IEEE.

From my experience this sequence synthesis method is not capable of synthesising all possible correlation matrices; in fact, quite obviously all Bernoulli random vectors are forced to have unit non-central second moment, i.e., any random vector ${\bf x}$ defined over the configuration space $\Omega = \{-1,1\}^n$ is so that $\forall i \in [n], \mathbb{E}[x^2_i] = 1$ (straightforward, by applying the definition of expectation and using conditioning over $x_j : j \neq i$). So for sure you will not be able to design such a sequence with arbitrary correlation.

Answer 3

Given $C$, the task is to either generate a Bernoulli random vector with covariance matrix $C$ or report that it's impossible to do so. What follows is a 2-step solution:

Step 1: Find consistent marginal and pairwise probabilities.

Following Robert Israel, we seek probabilities $\{p_i\}$ and $\{p_{ij}\}$ which are consistent with $C$, that is:

(1) $p_i(1-p_i)=C_{ii}$ for every $i$

(2) $\min\{p_i,p_j\}\geq C_{ij}+p_ip_j\geq\max\{0,p_i+p_j-1\}$ for every $i,j$ with $i\neq j$

(3) $p_{ij}=C_{ij}+p_ip_j$ for every $i,j$ with $i\neq j$

If no consistent probabilities exist, we report "impossible." To find the probabilities, Robert suggests you consider each of the $2^n$ possibilities for $\{p_i\}$ that satisfy (1) until you find one that satisfies (2), and then select $\{p_{ij}\}$ according to (3). In practice, it may be better to attempt locally minimizing

$$\sum_{i=1}^n\big(p_i(1-p_i)-C_{ii}\big)^2 \quad\mbox{subject to}\quad (2)$$

Presumably, it would be good to somehow initialize $p_i$ in terms of $C_{ii}$. Assuming consistency, you know the value of this program is 0, so you can re-initialize if you find a suboptimal local minimum. I haven't attempted to run this, so I can't speak to it's performance.

Step 2: Simulate random variables from marginal and pairwise probabilities.

Here, the big idea is to produce characteristic functions $\chi_i\colon \Omega\rightarrow\{0,1\}$ such that $\|\chi_i\|^2=p_i$ for each $i$ and $\langle \chi_i,\chi_j\rangle =p_{ij}$ for each $i,j$ with $i\neq j$. Then if we draw $\omega$ uniformly from $\Omega$, $\{\chi_i(\omega)\}$ will be the desired random vector. One may attempt to find such characteristic functions by seeking a non-negative matrix factorization of the matrix $V$ defined by $V_{ii}=p_i$ and $V_{ij}=p_{ij}$.

Answer 4

I'm assuming "Bernoulli vector" here means a vector with all entries in $\{0,1\}$. Thus we have the configuration space $\Omega = \{0,1\}^n$.

For example, let's try the case $n=2$. The covariance matrix is $$ C = \pmatrix{ p_1 (1-p_1) & p_{12} - p_1 p_2\cr p_{12} - p_1 p_2 & p_2 (1-p_2)\cr}$$ where $p_1 = P(X_1 = 1)$, $p_2 = P(X_2 = 1)$, $p_{12} = P(X_1 = 1, X_2 = 1)$. We have $$\min(p_1, p_2) \ge p_{12} \ge \max(0, p_1 + p_2 - 1)$$ Given $C_{11}$ with $0 \le C_{11} \le 1/4$ there are two possible values of $p_1$ with $C_11 = p_1 (1-p_1)$, and similarly for $C_{22}$ and $p_2$. For each of the four resulting pairs $(p_1, p_2)$ we can then check whether $$ \min(p_1, p_2) \ge C_{12} + p_1 p_2\ge \max(0, p_1 + p_2 - 1) $$ (but by symmetry, only two of the four need to be checked).

Answer 5

The Bernoulli distribution is characterized by a single parameter, $p$. If the covariance matrix $C$ is valid, it will completely characterize the distributions of $n$ random variables and thus it is possible to generate Bernoulli random vectors.

The following are the conditions for $C$ to be a valid covariance matrix:

• It should be a diagonal matirx.
• All the diagonal elements must be of the form $p*(1-p) \, \forall \, p \, \epsilon [0,1]$.

Your function does not need an input of sample Bernoulli random vectors.