Q:

How to sample uniformly from singular matrices

I would like to uniformly sample from all singular $n$ by $n$ Bernoulli matrices (that is each entry is $1$ or $0$ with probability $1/2$). I could of course just sample from all $n$ by $n$ Bernoulli matrices and reject those that are non-singular over $\mathbb{R}$ but for any moderate n that is extremely inefficient.
Is there an efficient way to do this?
Previously asked at https://stackoverflow.com/questions/21191831/how-to-uniformly-sample-from-singular-matrices where an answer over finite fields was given. 

A:

[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.)
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.
[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 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,...,40$, I performed a series of Monte Carlo experiments to estimate $p$.  Here are the results:

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$.
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).

