How to sample uniformly from singular matrices

Question

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.

Answer 1

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

Answer 2

This is an MCMC algorithm for uniform sampling over singular $n$ by $n$ Bernoulli matrices.

Let $H$ (for "hypercube") be the set of all 0/1 vectors of length $n$.

One step of the MCMC algorithm is as follows:

1. Generate an $(n-1)$ by $n$ matrix $A$, filled with 0/1 iid Bernoulli samples. This will be the first $n-1$ rows of our proposal matrix.

2. Find $K$, the kernel of $A$. Let $r$ be the rank of $K$ (and note that $r>0$).

3. Consider extending $A$ to an $n$ by $n$ matrix by adding a final 0/1 row. Let $L$ be the number of possible completions.

   1. If $r>1$, then $L=2^n$. Generate a proposal matrix $M$ with the same first $n-1$ rows as $A$, and the last row selected uniformly at random from $H$ (transposed).

   2. If $r=1$, let $L=|K\cap H|$. Generate a proposal matrix $M$ with the same first $n-1$ rows as $A$, and the last row selected uniformly at random from $K\cap H$. (We will discuss how to enumerate $K\cap H$ below.)

4. Let $L'$ be the corresponding value for $L$ for the previous accepted proposal. Let $p=\min(1, L/L')$. Accept the new proposal matrix with probability $p$.

(To generate the initial proposal, we can just follow steps 1-3.)

This appears to provide a straightforward MCMC algorithm for uniform sampling from singular Bernoulli matrices. All of the steps are efficient and straightforward except for one: enumerating $K \cap H$. How can we compute that in a (somewhat) practical way?

Observe that the probability that $r>1$ is less than the probability that a uniform Bernoulli matrix is singular. My earlier post on this question addresses suggests that this probability is on the order of $2^{-n}$ for large $n$ (e.g. $n=100$). Therefore, we will essentially always be in the case $r=1$.

We can find an element in $K\cap H$ by using a mixed (binary) integer linear program. Specifically, we try to find a binary integer vector that is in the kernel $K$ (which is a linear constraint). We only need a feasible solution (i.e. "phase 1" of the MILP); the objective function doesn't matter. This produces a single element of $K\cap A$ (or shows that no such element exists).

Note that, for any $z\in H$, we can construct a linear function $f_z$ on $R^n$ that coincides with the Hamming distance from any on vector $H$ to $z$, to wit: $f_z(x)=\sum_i^{n}(z_i+(-1)^z_i x_i)$. That is, for any $x,z\in H$, $f_z(x)$ is equal to the Hamming distance between the vectors.

This provides us with a method to enumerate all elements of $K\cap A$. First, we use the MILP to enumerate one solution, $z\in H$. Then we add a linear constraint to the MILP forcing the Hamming distance from $z$ to be $\geq 1$. This removes $z$ from the feasible set but keeps all other integer solutions. We keep repeating this until the MILP finds no feasible binary solutions. This procedure will enumerate $K\cap H$.

Note that the all-zero vector provides us with an initial feasible solution (so $|K\cap H|\geq 1$). In the typical case that the first $n-1$ rows are distinct and non-zero, there are at least $n$ solutions: the all-zeroes vector, and each of the other $n-1$ rows. Therefore, we can usually skip the first $n$ steps and just add those constraints to the MILP directly.

For a particularly bad $A$, this solution may take exponentially much time. For example, if $A$ is the first $n-1$ rows of an identity matrix, then $r=1$ but $|K\cap H|=2^{n-1}$, so the outer loop may take exponentially many steps before terminating. Moreover, the MILP might take exponentially long to identify a single feasible solution.

In practice (i.e. for a random $A$), I imagine that the $n$ solutions mentioned above are almost certainly the only solutions, so we probably only need to solve a single MILP (which will prove that there are no more feasible solutions). The practical question then becomes how long this MILP takes. I think we would just need to try it to find out. In my experience, CPLEX tends to be much faster than GLPK at solving MILPs, but I don't have a license handy to try it out.

There are a few optimizations that might be worth mentioning.

1. We could modify the MCMC algorithm so that with probability 1/2 it proceeds as above, but with probability 1/2 it reuses the same $A$ and just resamples the last row from $K\cap H$. Since $K\cap H$ has already been enumerated, and since these states are all equally likely, this second possibility takes essentially no compute time and MCMC always accepts. The net effect is much faster "local" mixing but longer scale dependence on $A$. Whether that is a good idea or not depends on what the samples are being used for.
2. Similarly, with a certain probability we can permute the rows, columns, or transpose the matrix (or some combination of all three).

Finally, it might be worth mentioning an alternate approach. We might call the previous suggestion the "kernel/MILP" approach, and the following suggestion the "cokernel/dynamic program" approach. (I'm not sure if I should break this off as a separate post, but the thought of providing 3 separate "answers" to this question made me feel sheepish.)

Note that if the matrix $M$ is singular, then there exists some row that is a linear combination of the other rows; therefore, there is some row vector $v\in R^n$ such that $vA=0$. I believe that in the preponderance of cases, $v$ is extremely sparse and has small entries. For example, if a row is all zeros, then $v$ has weight 1 and coefficient 1; if two rows are co-linear, then $v$ has weight 2 and coefficients $1$ and $-1$. (By "weight", I mean the number of non-zero entries in $V$.)

Suppose that we make a proposal for $v\in Z^n$, where we (strongly) bias our distribution to favor sparse vectors with small coefficients. For all the zero entries in $v$, we choose the entries of the corresponding rows of $M$ as uniformly Bernoulli samples. For the remaining rows, we can compute all valid sets of 0/1 entries such that the dot product with $v$ is zero. This involves solving the subset sum problem. This is NP-complete, but it is pseudo-polynomial-time solvable (with a dynamic program). We have "stacked the deck" in our favor by sampling from vectors with small entries, so this step will be efficient (that is, we can achieve polynomial expected work).

If $v$ has weight $k$, the solution to the dynamic program tells us the set $S$ of valid $k$-tuples for each column; the number of consistent solutions is then $S^n$. We can then perform an MCMC random walk with acceptance probability proportional to $|S|$ (eliding details). I imagine that the dynamic program would be vastly faster than the MILP in practice (in addition to being more theoretically tractable), so an approach along these lines would probably be much faster.

By the way, it is possible to bound the maximum possible value for an entry of $v$. This is essentially equivalent to Hadamard's maximum determinant problem, but unfortunately the bound is $(n+1)^{(n+1)/2}2^{-n}$ (so the corresponding dynamic program would be exponentially large.)

That said, I have gotten stuck with some technical issues when trying to construct the MCMC, so I presented the "kernel/MILP" approach above.