[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:
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$.