Skip to main content
Extended numerical experiments and accelerated rejection algorithm
Source Link

[Edit: After improving the rejection sampling algorithm and running it on a more powerful computer, I was able to extend my earlier numerical experiments. Improvements are described below in square brackets.]

This post does not answer the original question, but discusses a sub-problem: what is the probability $p$ that an $n$ by $n$ Bernoulli matrix is singular? (This bears on the original problem because it gives a better sense of how effective rejection sampling is.)

This does not provide any theoretical speed-up, since the operations are all $O(n^3)$ (ignoring fast matrix multiplication tricks), but in practice it's very helpful because the constants are so much better for matrix multiplication than either of the decompositions. Unfortunately, if $k$ is too large, this operation can become numerically unstable. The instability does not effect correctness, but nullifies the speed benefit. I found that instability occurred when $k$ was around 7 on my machine; to be conservative, I set $k=4$ and got about a 4x increase in speed.

[Edit: I was able to improve this trick as follows: Note that the product of the 0/1 matrices is an integer matrix. Therefore, we can compute this product modulo $P$ for any prime $P$. If the product is invertible mod $P$, then it is also invertible over the reals. For large $P$, there is approximately a $1/P$ probability of a false positive. Because there are no numerical stability concerns operating mod $P$, we can take $k$ to be arbitrarily large; our only practical constraint is the $1/P$ probability of a false positive per 0/1 matrix. I took $P=2^16-15=65521$ and $k=100$; this provides a speed-up on the order of 15x over the naive method of computing the rank of each matrix.]

As Ofer Zeitouni mentioned in the comments above, it's straightforward to construct $O(n^2 2^{-n})$ lower bounds on $p$ by considering the occurrence of co-linear columns. More precisely, the probability that at least one column is all zero is $1-(1-2^{-n})^n$. The probability that all columns are non-zero and precisely two columns are identical is ($n$ choose $2$)$\times (1-2^{-n})^2\prod_{i=3}^{n}(1-(i-1)2^{-n})$. Adding together the probability of these two disjoint events gives a lower bound for $p$ for all $n\geq 2$.

[Edit: We can provide an asymptotic lower bound by considering the event that a single row or column is all zero, or a pair of rows or columns is equal, and treating all four events as independent. (The events are not independent, but become increasingly independent as $n\rightarrow \infty$.) The bound is $n(n+1)2^{-n}$.]

For $n=6,...,38$$n=6,...,40$, I performed a series of Monte Carlo experiments to estimate $p$. Here are the results: Probability that an n by n Bernoulli matrix is singular, for n=1 to 38Probability that an n by n Bernoulli matrix is singular, for n=1 to 40

The green line crosses the sample means. (I computed error bars by assuming a beta distribution on $p$, with an initial uniform prior. The error bars then refer to two standard deviations above or below.) The last few values are fairlysomewhat uncertain; for example, for $n=38$$n=40$ I examined 474,061,83214.7 billion matrices and found that only 6 were24 singular ones. The blue line "Lower ("Lower"asymp)" is the asymptotic lower bound discussed above; the red line "Upper ("Upper"asymp)" is $2^{-n/2}$, i.e. I set the $o(1)$ term to zero. Ignoring the $o(1)$ term means that the "upper" bound is not actually an upper bound for $n<24$.

This post does not answer the original question, but discusses a sub-problem: what is the probability $p$ that an $n$ by $n$ Bernoulli matrix is singular? (This bears on the original problem because it gives a better sense of how effective rejection sampling is.)

This does not provide any theoretical speed-up, since the operations are all $O(n^3)$ (ignoring fast matrix multiplication tricks), but in practice it's very helpful because the constants are so much better for matrix multiplication than either of the decompositions. Unfortunately, if $k$ is too large, this operation can become numerically unstable. The instability does not effect correctness, but nullifies the speed benefit. I found that instability occurred when $k$ was around 7 on my machine; to be conservative, I set $k=4$ and got about a 4x increase in speed.

As Ofer Zeitouni mentioned in the comments above, it's straightforward to construct $O(n^2 2^{-n})$ lower bounds on $p$ by considering the occurrence of co-linear columns. More precisely, the probability that at least one column is all zero is $1-(1-2^{-n})^n$. The probability that all columns are non-zero and precisely two columns are identical is ($n$ choose $2$)$\times (1-2^{-n})^2\prod_{i=3}^{n}(1-(i-1)2^{-n})$. Adding together the probability of these two disjoint events gives a lower bound for $p$ for all $n\geq 2$.

For $n=6,...,38$, I performed a series of Monte Carlo experiments to estimate $p$. Here are the results: Probability that an n by n Bernoulli matrix is singular, for n=1 to 38

The green line crosses the sample means. (I computed error bars by assuming a beta distribution on $p$, with an initial uniform prior. The error bars then refer to two standard deviations above or below.) The last few values are fairly uncertain; for example, for $n=38$ I examined 474,061,832 matrices and found that only 6 were singular. The blue line ("Lower") is the lower bound discussed above; the red line ("Upper") is $2^{-n/2}$, i.e. I set the $o(1)$ term to zero. Ignoring the $o(1)$ term means that the "upper" bound is not actually an upper bound for $n<24$.

[Edit: After improving the rejection sampling algorithm and running it on a more powerful computer, I was able to extend my earlier numerical experiments. Improvements are described below in square brackets.]

This post does not answer the original question, but discusses a sub-problem: what is the probability $p$ that an $n$ by $n$ Bernoulli matrix is singular? (This bears on the original problem because it gives a better sense of how effective rejection sampling is.)

This does not provide any theoretical speed-up, since the operations are all $O(n^3)$ (ignoring fast matrix multiplication tricks), but in practice it's very helpful because the constants are so much better for matrix multiplication than either of the decompositions. Unfortunately, if $k$ is too large, this operation can become numerically unstable. The instability does not effect correctness, but nullifies the speed benefit. I found that instability occurred when $k$ was around 7 on my machine; to be conservative, I set $k=4$ and got about a 4x increase in speed.

[Edit: I was able to improve this trick as follows: Note that the product of the 0/1 matrices is an integer matrix. Therefore, we can compute this product modulo $P$ for any prime $P$. If the product is invertible mod $P$, then it is also invertible over the reals. For large $P$, there is approximately a $1/P$ probability of a false positive. Because there are no numerical stability concerns operating mod $P$, we can take $k$ to be arbitrarily large; our only practical constraint is the $1/P$ probability of a false positive per 0/1 matrix. I took $P=2^16-15=65521$ and $k=100$; this provides a speed-up on the order of 15x over the naive method of computing the rank of each matrix.]

As Ofer Zeitouni mentioned in the comments above, it's straightforward to construct $O(n^2 2^{-n})$ lower bounds on $p$ by considering the occurrence of co-linear columns. More precisely, the probability that at least one column is all zero is $1-(1-2^{-n})^n$. The probability that all columns are non-zero and precisely two columns are identical is ($n$ choose $2$)$\times (1-2^{-n})^2\prod_{i=3}^{n}(1-(i-1)2^{-n})$. Adding together the probability of these two disjoint events gives a lower bound for $p$ for all $n\geq 2$.

[Edit: We can provide an asymptotic lower bound by considering the event that a single row or column is all zero, or a pair of rows or columns is equal, and treating all four events as independent. (The events are not independent, but become increasingly independent as $n\rightarrow \infty$.) The bound is $n(n+1)2^{-n}$.]

For $n=6,...,40$, I performed a series of Monte Carlo experiments to estimate $p$. Here are the results: Probability that an n by n Bernoulli matrix is singular, for n=1 to 40

The green line crosses the sample means. (I computed error bars by assuming a beta distribution on $p$, with an initial uniform prior. The error bars then refer to two standard deviations above or below.) The last few values are somewhat uncertain; for example, for $n=40$ I examined 14.7 billion matrices and found 24 singular ones. The blue line "Lower (asymp)" is the asymptotic lower bound discussed above; the red line "Upper (asymp)" is $2^{-n/2}$, i.e. I set the $o(1)$ term to zero. Ignoring the $o(1)$ term means that the "upper" bound is not actually an upper bound for $n<24$.

Source Link

This post does not answer the original question, but discusses a sub-problem: what is the probability $p$ that an $n$ by $n$ Bernoulli matrix is singular? (This bears on the original problem because it gives a better sense of how effective rejection sampling is.)

Let's begin by mentioning a useful trick for accelerating Monte Carlo (or rejection) sampling of singular matrices. The straightforward technique would be: generate a random $n$ by $n$ Bernoulli matrix $M_1$, check if the rank is $<n$, repeat. Checking the rank involves either an QR decomposition or an SVD, both of which are kind of expensive.

Instead, we can generate $k$ candidate matrices $M_1, M_2,...,M_k$, compute their product, and check if the product is singular. If the product is singular, the check each of the $M_i$. If it's non-singular, then all the $M_i$ are guaranteed to be non-singular, so you can skip all those SVDs or QR decompositions. Because the $M_i$ are almost always non-singular for large $n$, this is a large acceleration.

This does not provide any theoretical speed-up, since the operations are all $O(n^3)$ (ignoring fast matrix multiplication tricks), but in practice it's very helpful because the constants are so much better for matrix multiplication than either of the decompositions. Unfortunately, if $k$ is too large, this operation can become numerically unstable. The instability does not effect correctness, but nullifies the speed benefit. I found that instability occurred when $k$ was around 7 on my machine; to be conservative, I set $k=4$ and got about a 4x increase in speed.

As Ofer Zeitouni mentioned in the comments above, it's straightforward to construct $O(n^2 2^{-n})$ lower bounds on $p$ by considering the occurrence of co-linear columns. More precisely, the probability that at least one column is all zero is $1-(1-2^{-n})^n$. The probability that all columns are non-zero and precisely two columns are identical is ($n$ choose $2$)$\times (1-2^{-n})^2\prod_{i=3}^{n}(1-(i-1)2^{-n})$. Adding together the probability of these two disjoint events gives a lower bound for $p$ for all $n\geq 2$.

For an asymptotic upper bound, a paper by Bourgain, Vu and Wood show that $p<(1/\sqrt{2}+o(1))^n)$ (see Corollary 3.3 here).

Note that on a log scale, the lower bound $O(2^{-n})$ and the upper bound $O(2^{-n/2})$ are not tight. Ofer Zeitouni mentioned that it was widely believed that $p$ asymptotically approaches the lower bound. Here are a few numerical experiments that support this belief.

For $n=1,...,5$, it is easy to compute $p$ exactly. The results are:

n, p=singular/total
1, 1/2              =0.500000
2, 10/16            =0.625000
3, 338/512          =0.660156
4, 42976/65536      =0.655762
5, 21040112/33554432=0.627044

For $n=6,...,38$, I performed a series of Monte Carlo experiments to estimate $p$. Here are the results: Probability that an n by n Bernoulli matrix is singular, for n=1 to 38

The green line crosses the sample means. (I computed error bars by assuming a beta distribution on $p$, with an initial uniform prior. The error bars then refer to two standard deviations above or below.) The last few values are fairly uncertain; for example, for $n=38$ I examined 474,061,832 matrices and found that only 6 were singular. The blue line ("Lower") is the lower bound discussed above; the red line ("Upper") is $2^{-n/2}$, i.e. I set the $o(1)$ term to zero. Ignoring the $o(1)$ term means that the "upper" bound is not actually an upper bound for $n<24$.

It's a bit strange, but the maximum value of $p$ appears to occur at $n=3$ (that is, $p$ is not monotonically decreasing, as you might expect).