# Probability two products are equal

I am interested in the following simple looking problem on which I am stuck. Let $M$ be a fixed $m$ by $n$ matrix with $\pm1$ elements. Let $x$ and $y$ be two independently sampled random $n$-dimensional vectors whose elements are chosen i.u.d. from $\{-1,1\}$.

Assuming that $m<n$ and both $m$ and $n$ are large, can we give a good estimate or bounds for $P(Mx = My)$?

Clearly $P(Mx= My) \geq P(x=y) = 2^{-n}$.

Both $Mx$ and $My$ have elements which are distributed as a simple symmetric random walk and have covariance matrix $MM^T$.

It feels like $\det(MM^T)$ should be part of any solution but this is just a guess and it doesn't work directly. I also attempted to solve the problem by relating the distributions of $Mx$ and $My$ to multivariate Gaussians but this didn't work out.

To narrow the question a little, let us assume that rank$(M) = m$.

Is the question more tractable if all the rows of $M$ are orthogonal?

• You are asking for the probability that $x-y$ is in the null-space of $M$. Since $x-y$ has entries in $\{-2,-1,0,1,2\}$, I first will try to understand how many vectors of this kind can be in the null-space of $M$, for example using dimensions and linear independence reasoning... – user40023 Jan 4 '16 at 11:14
• @Fry Thank you! I think $x-y$ has entries in $\{-2,0,2\}$ however. – Lembik Jan 4 '16 at 11:17
• One thing that jumps out at me immediately is that the rank of $M$ or one of its submatrices, viewed in isolation, might not be particularly helpful especially if it is much less than $n$. This is because for two matrices $M_1$ and $M_2$ with the same rank, depending on how the entries are shared by the subsets of rows, the answer could be dramatically different. – Yi Liu Jan 13 '16 at 11:37
• @YiLiu This is very interesting. Would you say the same thing about orthogonality of the rows? Is there some other property that might be more useful? For example, if we ensure that all the columns are distinct as well even up to sign flips? – Lembik Jan 13 '16 at 11:43
• After thinking about this spread across weeks, I am now convinced that the exact version of this problem is #P-hard. I am a novice and have not been able to find a reduction, but I think it deserves the attention of experts. – Yi Liu Jan 19 '16 at 2:27

I came across this problem a few days ago, and must say it was a lot of fun to think about. Thank you!

The following is, I believe, a quite complete solution of the problem for a "typical", general $M$, making no assumptions about the relative size of $n$ and $m$ (even $m>n$ is allowed), or about the rank of $M$.

We want to find the probability $P(Ms=0)$, where the $n$ components of the vector $s=\frac{x-y}2$ are i.i.d. random variables which have values $1$ or $-1$ with probability $\frac1 4$ and value $0$ with probability $\frac1 2$. $P(Ms=0)$ can then be expressed as the expectation value of a product of Kronecker deltas for the $m$ components of $Ms$: $$P(Ms=0) = \left\langle \prod_{i=1}^m \delta_{0,\sum_{j=1}^nM_{ij}s_j} \right\rangle_s \ , \hspace{2cm} (1)$$ where $\left\langle \cdots \right\rangle_s$ indicates the average over all $s_j$.

For any integer $t$ with $|t|\leq n$ we can substitute $\delta_{0,t} = \frac1{n+1} \sum_{r=0}^n \omega^{rt}$ with the primitive $(n+1)^{\rm th}$ root of unity $\omega = {\rm e}^\frac{2\pi{\rm i}}{n+1}$ ("discrete Fourier transform"), so that $(1)$ becomes $$P(Ms=0) = \left\langle \prod_{i=1}^m \left( \frac1{n+1} \sum_{r_i=0}^n \omega^{r_i\sum_{j=1}^n M_{ij}s_j} \right) \right\rangle_s = \left\langle \prod_{i=1}^m \prod_{j=1}^n \omega^{M_{ij} r_i s_j} \right\rangle_{r,s} \ , \hspace{2cm} (2)$$ where in the second step the $r_i$ have been re-interpreted as i.u.d. random variables with values in $\{0,...,n\}$, so that $\left\langle \cdots \right\rangle_{r_i} = \frac1{n+1} \sum_{r_i=0}^n (\cdots)$, and $\left\langle \cdots \right\rangle_{r,s}$ indicates the simultaneous average over all $r_i$ and $s_j$.

Now, obviously $P(Ms=0)$ has not the same value for all $(-1,1)$-matrices $M$ of a given size $m \times n$. As an example, if $M$ has rank $1$ (i.e., each of its column vectors is equal to either $v$ or $-v$, with some fixed $(-1,1)$-vector $v$) then the problem is equivalent to the case $m=1$, and the result is $$P(Ms=0) = 2^{-n} \sum_{l=0}^{\lfloor n/2 \rfloor} 2^{-2l} \begin{pmatrix} {n} \\ {2l} \end{pmatrix} \begin{pmatrix} {2l} \\ {l} \end{pmatrix} \ , \hspace{2cm} (3)$$ whereas usually it will be much smaller (for instance, if $M$ has rank $n$ then trivially $P(Ms=0) = 2^{-n}$). Thus $(2)$ cannot be evaluated without further assumptions on $M$. Alternatively, we can try to calculate a "typical" value of $(2)$, assuming that $M$ itself is a random sample out of a distribution of matrices with i.i.d. matrix elements $M_{ij}$, and averaging $P(Ms=0)$ over this distribution. This approach leads to $$\Big\langle P(Ms=0) \Big\rangle_M = \left\langle \prod_{j=1}^n \prod_{i=1}^m \left\langle \omega^{M_{ij} r_i s_j} \right\rangle_M \right\rangle_{r,s} = \left\langle \prod_{j=1}^n \prod_{i=1}^m \frac{ \omega^{r_i s_j} + \omega^{-r_i s_j} }2 \right\rangle_{r,s} \ . \hspace{1cm} (4)$$ Note that we have also changed the order of the product operations over $i$ and $j$ here. Now the average over the $s_j$ can be performed easily: for each $j$, $$\left\langle \prod_{i=1}^m \frac{ \omega^{r_i s_j} + \omega^{-r_i s_j} }2 \right\rangle_{s_j} = \frac1 2 \left( 1 + \prod_{i=1}^m \frac{\omega^{r_i} + \omega^{-r_i}}2 \right) \ ,$$ which is independent of $j$, so that $$\Big\langle P(Ms=0) \Big\rangle_M = 2^{-n} \left\langle \left( 1 + \prod_{i=1}^m \frac{\omega^{r_i} + \omega^{-r_i}}2 \right)^n \right\rangle_{r} \ . \hspace{2cm} (5)$$ To evaluate the right hand side of $(5)$ we expand the $n^{\rm th}$ power, $$\Big\langle P(Ms=0) \Big\rangle_M = 2^{-n} \sum_{k=0}^n \begin{pmatrix} {n} \\ {k} \end{pmatrix} 2^{-km} \prod_{i=1}^m \left\langle \Big( \omega^{r_i} + \omega^{-r_i} \Big)^k \right\rangle_{r} \ , \hspace{2cm} (6)$$ and likewise the $k^{\rm th}$ power in $(6)$: $$\left\langle \Big( \omega^{r_i} + \omega^{-r_i} \Big)^k \right\rangle_{r} = \sum_{l=0}^k \begin{pmatrix} {k} \\ {l} \end{pmatrix} \left\langle \omega^{(k-2l)r_i} \right\rangle_{r} = \begin{cases} \begin{pmatrix} {k} \\ {k/2} \end{pmatrix} , \ k\ {\rm even} \\ \quad \ 0 \ , \quad k\ {\rm odd} \end{cases} \ , \hspace{2cm} (7)$$ where $\left\langle \omega^{(k-2l)r_i} \right\rangle_{r} = \delta_{k,2l}$ has been used. Inserting $(7)$ into $(6)$ then leads to $$\Big\langle P(Ms=0) \Big\rangle_M = 2^{-n} \sum_{l=0}^{\lfloor n/2 \rfloor} 2^{-2lm} \begin{pmatrix} {n} \\ {2l} \end{pmatrix} \begin{pmatrix} {2l} \\ {l} \end{pmatrix} ^m \ . \hspace{2cm} (8)$$ Note the close similarity to Eq. $(3)$. The $l=0$ term in $(8)$ yields $\big\langle P(Ms=0) \big\rangle_M = 2^{-n}$, which is the lower bound mentioned in the question. For large $n$, $m$, and $l$, Stirling's approximation gives to leading order $$2^{-2lm} \begin{pmatrix} {n} \\ {2l} \end{pmatrix} \begin{pmatrix} {2l} \\ {l} \end{pmatrix} ^m \sim \sqrt{\frac{n}{4\pi l(n-2l)}} \ {\rm e}^{n I(\frac{2l}n)} (\pi l)^{-\frac m2} \ , \hspace{2cm} (9)$$ where $I(p)$ denotes the "information entropy" of a 2-state system with probabilities $p$ and $1-p$: $$I(p) = - p \ln(p) - (1-p) \ln(1-p) \ . \hspace{2cm} (10)$$ Note that the divergence of the right hand side of $(9)$ for $l \to 0$ is an artefact due to Stirling's approximation; the original expression on the left hand side has a $l \to 0$ limit equal to $1$. To find the dominant term(s) in the sum of $(8)$, we look for the maximum of the exponential part of $(9)$ as a function of $l$, disregarding for the moment the square root. That is, setting $p = \frac{2l}n$ and $q = \frac{m}n$, we look for the maximum of $${\rm e}^{n \left( I(p) - \frac q 2 \ln\frac{\pi p n}2 \right)} \hspace{2cm} (11)$$ as a function of $p \in ]0,1]$. The sum in $(8)$ is either dominated (for large enough $q$) by the $l=0$ term $1$, or otherwise by a term with $l \approx \frac{pn}2$, where $p$ is a solution of the equation $$\frac q{2p} = I'(p) = \ln\frac{1-p}p \ . \hspace{2cm} (12)$$ If $q$ is larger than about $0.56$, $(11)$ has no solution $p \in ]0,1]$. For $q$ smaller than about $0.56$ there are two solutions, the smaller of the two corresponding to a local minimum of ${\rm e}^{n \left( I(p) - \frac q 2 \ln\frac{\pi p n}2 \right)}$ and the larger one to a local maximum. For large $n$, this local maximum dominates the sum in $(8)$ if and only if $I(p) - \frac q 2 \ln\frac{\pi p n}2 > 0$ for this $p$ value. For large $n$, the maximal $q$ that meets this criterion decreases logarithmically with increasing $n$, so in that case we may assume $\frac q{2p} \ll 1$ and thus obtain $p \approx \frac1 2$ from $(12)$. Inserting this back into the exponential expression $(11)$ yields an estimate of its maximum for large $n$, $m$, and $l$ $$\max_{p \gg \frac1n} \left\{ {\rm e}^{n \left( I(p) - \frac q 2 \ln\frac{\pi p n}2 \right)} \right\} \approx {\rm e}^{n \left( \ln2 - \frac q 2 \ln\frac{\pi n}4 \right)} \ . \hspace{2cm} (13)$$ Finally, we can combine this with $(8)$, $(9)$ and obtain, to leading order for large $n$, $$2^{n} \Big\langle P(Ms=0) \Big\rangle_M \approx 1 + \sqrt{\frac{2}{\pi n}} \ {\rm e}^{n \left( \ln2 - \frac q 2 \ln\frac{\pi n}4 \right)} = 1 + 2^{n-\frac1 2} \left( \frac{\pi n}4 \right)^{-\frac{qn+1} 2} \ . \hspace{2cm} (14)$$ As long as the exponent in $(13)$ is negative, i.e. $q > \frac{\ln4}{\ln\frac{\pi n}4}$ (or $m > n \frac{\ln4}{\ln\frac{\pi n}4}$), the first term $=1$ in $(14)$ dominates, which means that only the trivial solution $s=0$ of $Ms=0$ has nonnegligible probability. Otherwise, the second term dominates, quantifying how the number of solutions increases when $m$ becomes small enough compared to $n$.

• This is a wonderful answer. Thank you! Take the case where $q > \ln{4}/\ln(\pi n/4)$. I find it very interesting if there is no easily or efficiently identifiable property of $M$ which will determine if the probability is approximately of the form $2^{-n}$ (or indeed $c^{-n}$ for any constant $c>1$). Do you think that is the case? – Lembik Jan 13 '16 at 18:50
• Thanks a lot for your appreciation! I am not sure whether I fully understand your comment, but I believe the rank of $M$ might be such a property. Assume $m<n$ and let $m'$ be the rank of $M$. We can select $m'$ linearly independent rows of $M$ and form a new $m'\times n$-matrix $M'$. Then $M's=0$ is equivalent to $Ms=0$, i.e. in particular the probabilities are equal. My result above applied to $M'$ tells us that, for large $n$, among the matrices with $m'/n > \ln{4}/\ln(\pi n/4)$ only a negligibe fraction has a probability $P(M's=0)$ larger than $2^{-n}$. – Dierk Bormann Jan 13 '16 at 20:46
• On the other hand, the majority of matrices with $m'/n < \ln{4}/\ln(\pi n/4)$ have a probability of the approximate form $P(M's=0) \approx \left(\pi n/4\right)^{- m'/2} = c^{-n}$, with $c = \left(\pi n/4\right)^{m'/2n}$ so that $1<c<2$. Of course this $c$ is not a constant, but at least it varies more slowly than exponentially with $n$. – Dierk Bormann Jan 13 '16 at 20:55
• nice answer. it seems one can't go far without relaxing the fixed matrix assumption. – kodlu Jan 13 '16 at 21:04
• My follow up question was slightly against the spirit of your answer. Let me put it another way. Say you have a particular large matrix $M$ with approx. $10n/\ln{n}$ rows and $n$ columns with rank $10n/\ln{n}$, say. Is there some property of the matrix that one can test from which one could work out if the probability will be approximately $2^{-n}$ or from which one could give some upper bound? I was originally hoping that orthogonality might be sufficient to get small probability but I think kodlu argued that that wouldn't be enough. – Lembik Jan 13 '16 at 22:24

We note that the brute force way to determine the exact answer involves $O(m^2 2^n)$ operations, where you exhaustively evaluate and tabulate the results of $Mv$ for all $v \in \{-1, 1\}^n$, and you sum the squares of the counts. This answer improves upon the brute force answer by providing an intuitive elementary lower bound that can be quickly computed.

As a preliminary observation, for fixed integer $n$ and vector $v$, we count the number of ways in which of $v^t x = v^t y$ where $v, x, y \in \{-1,1\}^n$. From elementary combinatorics, This count is $$\sum_{i = 0}^{n} {n \choose i}^2$$, where $i$ ranges over the total counts of entries in which $x$ (and hence $y$) differs from fixed $v$.

Now we tackle the matrix case. Let $M$ be given. Let the column indices of M be partitioned into $J$ non-empty equivalence classes $K_1, ... K_J$ whereby two columns indices are in the same equivalence class iff they are +/-1 multiples of each other. In our notation, the $K_j$ are disjoint subsets of $\{1, ... n\}$, and $\displaystyle \sum_{j=1}^J |K_j| = n$. Denote each $M_j$ as the submatrix of $M$ restricted to the columns $K_j$, and denote the subvector of $x$ and $y$ restricted to columns in $K_j$ as $x_j$ and $y_j$. To result in $Mx = My$, it suffices that $M_j x_j = M_j y_j$ for all $j \in 1 \ldots J$. Because the criterion is not necessary, the subsequent bound we provide will not be tight. We proceed by manually counting, for each $j$, the number of times that $M_j x_j = M_j y_j$ holds for $x_j, y_j \in \{-1, 1\}^{|K_j|}$. Using the preliminary observation, this is $$\sum_{i=0}^{|K_j|} {|K_j| \choose i}^2$$

The above provides the number of ways by which $M_j x_j = M_j y_j$ for a single $j$. Across all columns, the total count of ways to specify $((x_1, \ldots x_J), (y_1, \ldots, y_J))$ is $$\prod_{j=1}^J \sum_{i=0}^{|K_j|} {|K_j| \choose i}^2$$. Expressed as a probability, there are $2^{2n}$ ways for $x$ and $y$ to take values. Thus, we arrive at the lower bound

$$P(Mx=My)\ge 2^{-2n} \prod_{j=1}^J \sum_{i=0}^{|K_j|} {|K_j| \choose i}^2$$

Given the equivalence classes, the running time of the above equation is $O(m+n)$. The time needed to identify the equivalence classes for the columns of $M$ is $O(mn \log(mn))$.

Examples:

Note that the bound is exact in some notable cases. When $M$ is square and has rank $n$, we have $J = n$, and the bound evaluates to: $$2^{-2n} \prod_{j=1}^J \sum_{i=0}^{|K_j|} {|K_j| \choose i}^2 = 2^{-2n }\prod_{j=1}^{n} (1+1) = 2^{-n}$$.

When $M$ has rank 1, we have $J=1$ and the computation reduces to the preliminary observation $$2^{-2n} \prod_{j=1}^J \sum_{i=0}^{|K_j|} {|K_j| \choose i}^2 = 2^{-2n }\sum_{i=0}^n {n \choose i}^2$$.

• Thank you. Is it possible to get an upper bound using this or some other approach do you think? – Lembik Jan 13 '16 at 9:57
• I can't think of a way to generate high quality upper bounds at the moment using this particular technique. – Yi Liu Jan 13 '16 at 12:36

Not a complete answer, too long as a comment: Firstly, the difference of $x-y$ has entries in $\{-2,0,2\}$ but we can consider vectors over $\{-c,0,c\}$ for any $c>0$ by scaling, since all we care about is whether they belong to the nullspace.

What is the rank distribution of $M$? If $m\ll n$ it is likely to be $m$ or very nearly $M$. If we impose the uniform probability on the space of $m\times n$ $\{\pm 1\}$ matrices, and let $A$ have rows $a_1,\ldots,a_m$ ($M$ can be taken as a realization of $A$) we can argue recursively and find $$Pr[rank([a_1|\cdots|a_m])^{T}=k],\quad 1\leq k\leq m.$$

Essentially think of it as a discrete Markov chain where $X_v=rank([a_1|\cdots|a_v])$ and given $X_v=k,$ we have $X_{v+1}=k+1$ with probability $2^{n-k}/2^n=2^{-k}$ and $X_{v+1}=k$ otherwise. The relative size of the orthogonal complement of the $k-$dimensional subspace spanned by $a_1,\ldots,a_v$ determine this probability distribution for $v=m.$

From here, the question becomes, for a fixed rank $k$ what is the probability that a $\{-1,0,1\}$ vector lies in the subspace determined by that rank. If you could stay in $\{-1,1\}$ (equivalently over $F_2$) the number of subspaces of a given rank is also known. Here, your given $\{-1,0,1\}$ vector can actually be taken to be on the unit sphere to answer this question since if the Hamming weight of the vector is $w,$ we can normalize to a $\{-1/\sqrt{w},0,1/\sqrt{w}\}$-vector. Such a unit vector, if it was drawn at random from the set of all unit vectors would have an essentially gaussian distribution of Hamming weight, with standard deviation $\sqrt{n}.$

I believe that concentration of measure would then imply that for $m\approx \log n$ with a fixed nonzero probability all unit vectors would be in a subspace of dimension $\log n$ and you'd only need to scale by the total number of subspaces of that dimension.

There are experts here who can answer this properly, I am sure.

• The probability could be much larger than $2^{-n}$ couldn't it, depending on the entries of $M$ and how large $m$ is relative to $n$? $M$ is not a random matrix in my question. Also, as rank($M$) < n, we are guaranteed that there exist non-zero vectors in its nullspace unless I have misunderstood what you meant. – Lembik Jan 4 '16 at 13:43
• @Lembik you're right! I shouldn't post at 1am... I'll modify the answer. – kodlu Jan 4 '16 at 22:58
• I made a couple of additions to the question which I hope might make it more tractable. I am particularly interested when $m$ is not too small compared to $n$. Let us say $m > n^c$ for some $0< c < 1$ – Lembik Jan 8 '16 at 11:34

I will consider $n$ even for simplicity, and try to argue for an approximate answer. The difference $z=x-y$ has independent identically distributed entries $z_i$ in $\{-2,0,2\}$ with distribution $Pr[z=-2]=Pr[z=+2]=1/4$ and $Pr[z=0]=1/2.$ Consider an arbitrary vector $a$ in $\{\pm 1\}.$ Then the sum $$\langle a, z\rangle=\sum_{i=1}^n a_i z_i$$ has terms $u_i=a_i z_i$ with the same distribution. Thus $\langle a, z\rangle$ has a multinomial distribution. Let $Z$ be the random variable denoting the number of $u_i=0.$ Then $$Pr[\langle a,z\rangle=0]=\sum_{0\leq Z\leq n:Z\equiv n~(mod~2)} 2^{-Z} \binom{n}{Z} Pr[number(u_i=+1)=(n-Z)/2],$$ since the +1 and -1 must balance. From here we get $$Pr[\langle a,z\rangle=0]=\sum_{0\leq Z\leq n:Z\equiv n~(mod~2)} 2^{-n} \binom{n}{Z} \frac{\binom{n-Z}{(n-Z)/2}}{2^{n-Z}},$$ which simplifies to $$Pr[\langle a,z\rangle=0]=\sum_{0\leq Z\leq n:Z\equiv n~(mod~2)} \frac{\binom{n}{Z} \binom{n-Z}{(n-Z)/2}}{2^{2n-Z}}=\frac{\Gamma(\frac{1}{2}+n)}{\sqrt{\pi}\Gamma(n+1)}\stackrel{\triangle}{=} q_n$$
and this is asymptotic to $1/\sqrt{\pi n}$ as $n$ gets large.i

The answer would be $$Pr[\langle a_i,z\rangle=0],\quad i=1,\ldots,m$$which would be $q_n^m$ if independence could be assumed.

Edit: Under the assumption that the rows $a_i$ are orthogonal, consider two such rows $a,a'$ and let $W$ be the number of coordinates where $a$ and $a'$ have matching entries (both -1 or both +1). Since $a,a'$ are orthogonal, $W=n/2$ and $n$ must be even under the orthogonality assumption. Let $I_W\subset\{1,\ldots,n\}$ be the support of the coordinates where equality is achieved. Since $a$ and $a'$ are orthogonal we can (by multiplying coordinates of $a$ and $a'$ by $-1$ if necessary) assume $a$ to be the all 1 vector and $a'$ to have (say) its first $n/2$ coordinates +1 and the rest -1.

Given $a$ is orthogonal to $z$, it is clear that if the sum of the coordinates of $z$ on the first $n/2$ components is $S$, it is $-S$ on the second $n/2$ components. Multiplying componentwise by entries of $a'$ this would give $2S$. Thus, for both $a$ and $a'$ to be orthogonal to $z$ the coordinates of $z$, $S=0$ is needed. This is a similar sum as before, but on $n/2$ coordinates, so $$Pr[\langle a,z\rangle=0]\times Pr[\langle a',z\rangle=0~|~\langle a,z\rangle=0]=q_n q_{n/2}$$

Proceeding by induction with the same kind of coordinate partitioning, and assuming $2^{m-1}$ divides $n$, one should get $$Pr[\langle a_1,z\rangle=0]\times Pr[\langle a_2,z\rangle=0~|~\langle a,z\rangle=0] \times \cdots \times \quad \quad \quad$$ $$Pr[\langle a_m,z\rangle=0~|~\langle a_1,z\rangle=0, \ldots, \langle a_{m-1},z\rangle=0]=q_n q_{n/2} q_{n/4} \cdots q_{n/2^{m-1}}$$

• @DouglasZare: do you have any idea how to proceed. Maybe a worst case distribution for the $a_i$'s? – kodlu Jan 5 '16 at 4:03
• In the first half of your answer before the Edit,I think there is a typo. $Pr[\langle a,z\rangle=0] \sim 1/\sqrt{2\pi\sigma^2 n}$ where $\sigma^2 = 1/2$ in this case. You are missing the $1/\sqrt{n}$ term. – Lembik Jan 9 '16 at 6:56
• Thank you for the edit. It seems your answer requires $m \approx \log{n}$. Is there anything we can do if $n >m > n^c$ for some $0<c<1$ ? That is if $2^{m-1} \gg n$. – Lembik Jan 9 '16 at 7:07
• I think you are now missing the $1/\sqrt{\pi}$ term. – Lembik Jan 9 '16 at 12:20
• @Lembik:since you find my non-answer useful, an upvote would be appreciated :-) – kodlu Jan 13 '16 at 22:46