text:
\begin{document}

\maketitle

\begin{abstract}
Recently Ermon et al. (2013) pioneered an ingenuous way to practically compute approximations to large scale counting or discrete integration problems by using random hashes. The hashes are used to reduce the counting problems into many separate discrete optimization problems. The optimization problems can be solved by an NP-oracle, and if they allow some amenable  structure then commercial SAT solvers or linear programming (LP) solvers can be used in lieu of the NP-oracle. In particular, Ermon et al. has shown that if the domain of integration is  $\{0,1\}^n$ then it is possible to obtain a $16$-approximation by this technique.

In many crucial counting tasks, such as computation of partition function of ferromagnetic Potts model, the domain of integration is naturally $\{0,1,\dots, q-1\}^n, q>2$. A straightforward extension of Ermon et al.'s work would allow a $q^2$-approximation for this problem, assuming the existence of an optimization oracle. In this paper, we show that it is possible to obtain a $(2+\frac2{q-1})^2$-approximation to the discrete integration problem, when $q$ is a power of an odd prime (a similar expression follows for even $q$). We are able to achieve this by using an idea of optimization over multiple bins of the hash functions, that can be easily implemented by inequality constraints, or even in unconstrained way. Also the burden on the NP-oracle is not increased by our method (an LP solver can still be used). Furthermore, we provide a close-to-4-approximation for the permanent of a matrix by extending our technique. Note that, the domain of integration here is the symmetric group. Finally, we provide memory optimal hash functions that uses minimal number of random bits for the above purpose. We are able to obtain these structured hashes without sacrificing the amenability of the NP-oracle. We provide experimental simulation results to support our algorithms.

\end{abstract}

\section{Introduction}
Large scale counting problems, such as computing the permanent of a matrix or computing the partition function of a graphical probabilistic generative model, come up often in variety of inference tasks. These problems can be, without loss of any generality, written as a discrete integration:  the summation of evaluations of a nonnegative function $w:\Omega \to \reals$ over all elements of $\Omega$:
\begin{align}\label{eq:main}
S_\Omega(w) \equiv \sum_{\sigma \in \Omega} w(\sigma).
\end{align}

These problems are computationally intractable because of the exponential (and sometime super-exponential) size of $\Omega$. A special case is of course the set of problems \#P, counting problems associated with the decision problems in NP. For example, one might ask how many variable assignments a given CNF (conjunctive normal form) formula satisfies. The complexity class \#P was defined by Valiant \cite{valiant1979complexity}, in the context of computing the permanent of a matrix. Indeed, the permanent of a matrix $A$ can be represented as,
\begin{align}\label{eq:perm}
{\rm Perm}(A) \equiv \sum_{\sigma \in S_n} \prod_{i=1}^n A_{i,\sigma(i)},
\end{align}
where $S_n$ is the symmetric group of $n$ elements and $A_{i,j}$ is the $(i,j)$-th element of $A$. Clearly, here $S_n$ is playing the role of $\Omega$, and   $w(\sigma) = \prod_{i=1}^n A_{i,\sigma(i)}$. Therefore computing permanent of a matrix is a canonical example of a problem defined by eq.~\eqref{eq:main}.

Similar counting problems arise when one wants to compute the partition functions of the well-known probabilistic generative models of statistical physics, such as the Ising model, or more generally the Ferromagnetic Potts Model \cite{potts1952some}. Given a graph $G(V,E)$, and a label-space $Q \equiv \{0,1,2,\dots, q-1\}$, the partition function  $Z(G)$ of the Potts model is given by,
\begin{align}\label{eq:part}
Z(G)  \equiv \sum_{\sigma \in Q^{|V|}}& \exp\Big(- \zeta\Big(J\sum_{(u,v) \in E} \delta(\sigma(u),\sigma(v)) + H\sum_{u \in V} \delta(\sigma(u),0)  \Big) \Big),
\end{align}
where $\zeta$, $J$ and $H$ are system-constants (representing the temperature, spin-coupling and external force),   $\delta(x,y)$ is the delta-function that is $1$ if and only if $x=y$ and otherwise $0$, and $\sigma$ represents a label-vector, where $\sigma(u)$ is the label of vertex $u$.
It has been shown that, under the availability of an NP-oracle, every problem in \#P can be approximated within a factor of $(1\pm \epsilon), \epsilon >0$, with high probability via a randomized algorithm \cite{stockmeyer1985approximation}. This result says \#P can be approximated by ${\rm BPP}^{\rm NP}$ and the power of an NP-oracle and randomization is sufficient.  However, depending on the weight function $w(\cdot)$,  eq.~\eqref{eq:main} may or may not be in \#P.

Recently, Ermon et al. proposed an alternative approach to solve these counting problems \cite{ermon2013taming,ermon2013optimization} by breaking them into multiple optimization problems. Namely, they  use families of hash functions $h: \Omega \to \tilde{\Omega}, |\tilde{\Omega}| < |\Omega|,$ and use a (possibly NP) oracle that can return the correct solution of the  optimization problem:
$
\max_{\sigma: h(\sigma) = a} w(\sigma).
$
We call this oracle a MAX-oracle. In particular, when $\Omega$ can be represented as $\{0,1\}^n$, and $h(\cdot)$ is a random hash function, assuming the availability of a MAX-oracle, Ermon et al. \cite{ermon2013taming} propose a randomized algorithm that approximates the discrete sum within a factor of  sixteen (a 16-approximation) with high probability.
Ermon et al. use simple linear sketches over $\ff_2$, i.e., the hash function $h_{A,b}:\ff_2^n \to \ff_2^m, A \in \ff_2^{m \times n}, b \in \ff_2^m$ is defined to be
\begin{align}\label{eq:hash}
h_{A,b}(x) = Ax + b,
\end{align}
where the arithmetic operations are over $\ff_2$. The matrix $A$ and the vector $b$ are randomly and uniformly chosen from the respective sample spaces. The MAX-oracle in this case simply provides solutions to the optimization problem:
$
\max_{\sigma \in \ff_2^n: A\sigma = b} w(\sigma).
$

The constraint space $\{\sigma \in \ff_2^n: A\sigma = b\}$ is nice since it is a coset of the nullspace of $A$, and experimental results showed them to be manageable by  optimization softwares/SAT solvers. In particular it was observed that being Integer Programming constraints, real-world instances are often solved quickly. Since the implementation
of the hash function heavily affects the runtime, it makes sense to keep constraints of the MAX-oracle as an affine space as above. These constraints are also called parity constraints. The idea of using such constraints to show reduction among class of problems appeared in several papers before, including \cite{sipser1983complexity,valiant1986np,gomes2006model,thurley2011approximation,gomes2007counting} among others.

The key property that the hash functions $\{h_{A,b}\}$ satisfy is that they are pairwise independent. This property can be relaxed somewhat - and in a subsequent paper Ermon et al. show that a hash family would work even if the matrix $A$ is sparse and random, thus effectively reducing the randomness as well as making the problem more tractable empirically \cite{ermon2014low}. Subsequently, Achlioptas and Jiang \cite{achlioptas2015stochastic} have shown another way of achieving similar guarantees. Instead of arriving at the set $\{\sigma \in \ff_2^n: A\sigma = b\}$ as a solution of a system of linear equations (over $\ff_2$), they view the set as the image of a lower-dimensional space. This is akin to the generator matrix view of a linear error-correcting code as opposed to the parity-check matrix view. This viewpoint allows their MAX-oracle  to solve just an unconstrained optimization problem.

\noindent{\bf Our contributions.} Note that, some crucial counting problems, such as computing the partition function of the Ferromagnetic Potts model of Eq.~\eqref{eq:part}, naturally have $\Omega= \ff_q^n$ for $q>2$. To use the binary-domain algorithm of \cite{ermon2013taming} for any $\Omega = \ff_q^n$, we need to use a look-up table to map $q$-ary numbers to binary\footnote{Also it is worth noting that while there exists polynomial time approximation (FPRAS) for the Ising model ($q=2$),  FPRAS for general Potts model ($q >2$) is significantly more challenging (and likely impossible \cite{goldberg2012approximating})}. In this process the number of variables (and also the number of constraints) increases by a factor of $\log q$. This makes the MAX-oracle significantly slower, especially when $q$ is large. Also, for the permanent problem, where $|\Omega|=\exp(n\log n)$,  this can create a big problem. It would be useful to extend the method of   \cite{ermon2013taming} for  $\Omega= \ff_q^n$ without increasing the number of variables. An extension of  the method  of \cite{ermon2013taming} for $q$-ary is possible: it  provides a $q^2$-approximation at best which is  particularly bad if $q$ is growing with $n$.
For the binary setting, it has been noted in \cite[section 5.3]{ermon2013taming} that the approximation ratio can be improved to any  $\alpha>1$ by increasing the number of variables, which extends to this $q$-ary setting. However this also results in an increase in number of variables by a factor of $\log_\alpha (q^2)$ which is undesirable.

Our first contribution in this paper is to provide a new algorithm by modifying and generalizing the algorithm of \cite{ermon2013taming}.

For any $\Omega = \ff_q^n, q$ is a power of prime, our algorithm provides a $4(1+\frac{1}{q-1})^2$-approximation, when $q$ is odd, and $4(1+\frac{2}{q-2})^2$-approximation, when $q>2$ is even, to the optimization problem of \eqref{eq:main} assuming availability of the MAX-oracle.
Our algorithm utilizes an idea of using optimization over multiple bins of the hash function that can be easily implemented via inequality constraints.  The constraint space of the MAX-oracle remains an affine space and can be represented as an integer linear program.

In general, for arbitrary $\Omega$, if represented as as $\{0,1\}^n$, the approximation factor is at best $16$ by the hashing technique. But by having it represented as $\{0,1,\dots, q-1\}^n$  the approximation factor can be improved by our technique.
Our multi-bin technique can also be use to extend the generator-matrix based algorithm of Achlioptas and Jiang \cite{achlioptas2015stochastic}. As a result, we need the MAX-oracle to only perform unconstrained maximization, as opposed to constrained. This lead to significant speed-up in the system, while resulting in the same approximation guarantees.

Secondly, we show that by using our technique and some modifications to the MAX-oracle, it is possible to obtain close-to-$4$-approximation to the problem of computing permanent of nonnegative matrices. The NP-oracle still is amenable to be implemented in an LP solver. It is to be noted that our idea of optimization over multiple bins is crucial here, since the straightforward generalization of Ermon et al.'s result would have given an approximation factor of $\Omega(n^2)$. While there already exists a polynomial time randomized approximation scheme ($1\pm \epsilon$-approximation) of permanent of a nonnegative matrix \cite{jerrum2004polynomial}, the runtime there is $\tilde{O}(n^{10})$. Since we are delegating the hard task to a professional optimization solver, our method can still be of interest here.

We note that it is possible to derandomize the hash families based on parity-constraints to the optimal extent while maintaining the essential properties necessary for their performance. Namely, it can be ensured that the hash family can still be represented as $\{x \mapsto Ax+b\}$ while using information theoretically optimal memory to generate them.

Finally, we show the performance of our algorithms to compute the partition function of the ferromagnetic Potts model by running experiments on both synthetic datasets and real-worlds datasets. We also use our algorithm to compute the Total Variation distance between two joint probability distributions over a large number of variables. The algorithm to compute permanent is also validated experimentally. All the experiments exhibit  good performance guarantees.

\noindent{\bf Organization.}
The paper is organized as follows. In Section~\ref{sec:back} we describe the technique by \cite{ermon2013taming}, and then elaborate  our new ideas and  main results. In Section~\ref{sec:wish}, we provide an improvement of the \texttt{WISH} algorithm by \cite{ermon2013taming} that lead to an improved approximation. We delegate an algorithm with unconstrained optimization oracle (similar to  \cite{achlioptas2015stochastic}) and its analysis to the appendix.  Section~\ref{sec:perm} is devoted  to computation of permanent of a matrix. In Section~\ref{sec:der}, we show how to optimally derandomize the hash function used in our algorithm.  The experimental results on computation of partition functions, total variation distance, and the permanent is provided in Section~\ref{sec:exp}.

\section{Background and our techniques}\label{sec:back}
In this section we describe the main ideas developed by \cite{ermon2013taming} and provide an overview of the techniques that we use to arrive at our new results.

First of all, notice that from \eqref{eq:main} we obtain:
$
S_\Omega(w) = \sum_{u\ge 0}u |\{\sigma \in \Omega: w(\sigma) =u \}|
 = \sum_{u\ge 0} |\{\sigma \in \Omega: w(\sigma) \ge u \}| = \sum_{u\ge 0} T(u),
$
where $T(u) \equiv |\{\sigma \in \Omega: w(\sigma) \ge u \}| $ is the tail distribution of weights and a nonincreasing function of $u$. Note that, $0 \le T(u) \le |\Omega|$. We can split the range of $T(u)$ into geometrically growing values $1, q, q^2, \dots, q^{n'}$ such that $q^{n'}\ge |\Omega|$. Let $\beta_i =u: T(u) = q^i, i=0, 1, \dots, n'$. Clearly $\beta_0 > \beta_1 > \dots > \beta_{n'}$. As we have not made any assumption on the values of the weight function, $\beta_i$ and $\beta_{i+1}$ can be far from each other and they are hard to bound despite the fact that $T(u)$  is monotonic in nature. On the other hand we can try to bound the area under the curve $T(u)$ by bounding the area of the slice between $\beta_i$ and $\beta_{i+1}$. This area is at least $q^i(\beta_i-\beta_{i+1})$ and  at most $q^{i+1}(\beta_i-\beta_{i+1})$. Therefore:
$
\sum_{i=0}^{n'-1} q^{i}(\beta_{i}-\beta_{i+1})+ q^{n'}\beta_{n'}  \le S_\Omega(w) \le \sum_{i=0}^{n'-1} q^{i+1}(\beta_{i}-\beta_{i+1})+ q^{n^{'}}\beta_{n'}$ which implies
\vspace{-0.1in}
\begin{align} \beta_0 +(q-1) \sum_{i=1}^{n'}q^{i-1}\beta_{i}  &\le S_\Omega(w) \le  \beta_0 +(q-1) \sum_{i=1}^{n'}q^{i}\beta_{i}.\label{eq:bound}
\end{align}
Hence $\beta_0 +(q-1) \sum_{i=1}^{n'}q^{i-1}\beta_{i}$ is a $q$-factor approximation of $S_\Omega(w)$ and if we are able to find a $k$ approximation of each value of $\beta_{i}$ we will be able to obtain a $kq$-factor approximation of $S_\Omega(w)$. In \cite{ermon2013taming}, subsequently the main idea is to estimate the coefficients $\{\beta_i, 0\le i \le n'\}$.

Now note that,
$
q^i = |\{\sigma \in \Omega: w(\sigma) \ge \beta_i \}|.
$
Suppose, using a random hash function $h:\Omega \to \{0,1, \dots, q^i-1\}$ we compute hashes of all elements in $\Omega$. The pre-image of an entry in $\{0,1, \dots, q^i-1\}$ is called the {\em bin} corresponding to that value, i.e., $\{\sigma: h(\sigma) = x\}$ is the bin corresponding to the value $x\in \{0,1, \dots, q^i-1\}$. In every  bin for the hash function, there is on average one element $\sigma$ such that $w(\sigma) \ge \beta_i$.  So for a randomly and  arbitrarily chosen bin if $\sigma^\ast = \max_{\sigma: h(\sigma) = x} w(\sigma)$, then $\sigma^\ast \ge \beta_i$ on expectation. However suppose one performs this random hashing $\ell = O(\log n')$ times and then take the aggregate (in this case the median) value of $\sigma^\ast$s. That is say, $\hat{\sigma^\ast} = {\rm median}(\sigma^\ast_1, \dots, \sigma^\ast_\ell)$.
Then by using the independence of the hash functions, it can be shown that the aggregate is an upper bound on $\beta_i$ with high probability.
Indeed, without loss of generality,  if we assume that the configurations within $\Omega$ are ordered according to the value of the function $w$, i.e., $w(\sigma_1) \ge w(\sigma_2) \ge \dots \ge w(\sigma_{|\Omega|})$ then we can take $\beta_i = w(\sigma_{q^i})$. If the hash family is pairwise independent,  then by using the Chebyshev inequality it can be shown that
$\hat{\sigma^\ast} \in [\beta_{i+c}, \beta_{i-c}]$ with high probability, $c \ge 2$. This lead to a $q^{2c}$-approximation for $S_\Omega(w)$. For $c=2$ this leads to the $16$-approximation, because \cite{ermon2013taming} identified $\Omega$ with $\ff_2^n$ and took $q=2$.
The \texttt{WISH} algorithm proposed by  \cite{ermon2013taming} makes use of the above analysis and provides a $16$-approximation of $S_w(\Omega)$. If we naively extend this algorithm by identifying $S_w(\Omega)$ with $\ff_q^n$ then the approximation factor we achieve is $q^{2c}, q>2, c\ge 1$. Note that, for $q=2$, it was not possible to take $c=1$, but as we will see later that it is possible to take $c=1$ when $q>2$, and  for $q=3$, this observation immediately gives a $9$-approximation to $S_w(\Omega)$.

Instead of using a straightforward  analysis for the $q$-ary case, in this paper we use a MAX-oracle that can optimize over multiple bins of the hash function. Using this oracle we proposed a modified \texttt{WISH} algorithm and call it \texttt{MB-WISH} (Multi-Bin \texttt{WISH}). Just as in the case of  \cite{ermon2013taming,ermon2014low}, the MAX-oracle constraints can be integer linear programming constraints and commercial softwares such as CPLEX can be used.

The main idea of using an optimization over multiple bins is that it boosts the probability that the $\sigma^\ast$ we are getting above is close to $b_i$. However if we restrict ourselves to the binary alphabet then (as will be clear later) there is no immediate way to represent such multiple bins in a compact way in the MAX-oracle. For the non-binary case, it is possible to represent multiple bins of the hash function as simple inequality constraints.
{\em This idea lead to an improvement in the approximation factor of $S_w(\Omega)$ to $4+\epsilon$, where $\epsilon$ decays to $0$ proportional to $q^{-1}$. Note that we need to choose $q$ to be a power of prime so that $\ff_q$ is a field.}

In \cite{achlioptas2015stochastic}, the bins (as described above) are produced as images of some function, and not as pre-images of hashes. Since we want the number of bins to  be $q^i$, this can be achieved by looking at images of $g: \ff_q^{n-i} \to \Omega$ where $|\{g(\sigma): \sigma \in \ff_q^{n-i} \}| = q^{n-i}.$ The rest of the analysis of \cite{achlioptas2015stochastic} is almost same as above. The benefit of this approach is that the MAX-oracle just has to solve an unconstrained optimization  here.
Implementing our multi-bin idea for this perspective of \cite{achlioptas2015stochastic} is not straight-forward as we can no longer use inequality constraints for this. However, as we show later, we found a way to combine bins here in a succinct way generalizing the design of $g$. As a result, we get the same approximation guarantee as in \texttt{MB-WISH}, with the oracle load heavily reduced (the details of this part can be found in the Appendix, we call this the \texttt{Unconstrained MB-WISH}).

Coming back to the discussion on \texttt{MB-WISH}, for computing  the {\em permanent}, the domain of integration  is the symmetric group $S_n$. However $S_n$ can be embedded in  $\ff_q^n$ for a $q>n$. Therefore we can try to use \texttt{MB-WISH} algorithm and same set of hashes on elements of $S_n$ treating them as $q$-ary vectors, $q>n$. We need to be careful though since it is essential that the MAX-oracle returns a permutation  and not an arbitrary vector. The modified MAX-oracle for permanents therefore must have some additional constraints. However those being affine constraints, it turns out MAX-oracle is still implementable in common optimization softwares quite easily.

For the analysis of  \cite{ermon2013taming,ermon2014low} to go through, we needed a family of hash functions that are pairwise independent\footnote{It is sufficient to have the hash family satisfy some weaker constraints, such as being pairwise negatively correlated.}. A hash family  $\Hc = \{h:\Omega \to \tilde{\Omega}\}$ is called uniform and pairwise independent if the following two criteria are met for a randomly and uniformly chosen $h$ from $\Hc$:

1) for every $x \in \Omega,$ $h(x)$ is uniformly distributed in $\tilde{\Omega}$ and 2) for any two $x,y \in \Omega$ and $u,v \in \tilde{\Omega}$, $\Pr(h(x) =u, h(y) =v) = \Pr(h(x) =u)\Pr(h(y) =v).$

By identifying $\Omega$ with $\ff_2^n$ (and $\tilde{\Omega}$ with $\ff_2^m$) and by using a family of hashes $\{x\mapsto h_{A,b}(x) = Ax + b: A \in \ff_2^{m \times n}, b \in \ff_2^m\}$ defined in \eqref{eq:hash}, Ermon et al. \cite{ermon2013taming} show the family to be pairwise independent and thereby achieve their objective.

Finally, the size of the hash family $\Hc$ determines how many random bits are required for the randomized algorithm to work. Although by defining the hash family in the above way Ermon et al. reduce the number of random bits from potentially $m2^n$ to $mn +m = m(n+1)$ bits (see, p.~3 of \cite{ermon2013taming}). {\em We show that it is possible to come up with hash functions of the same form (i.e., $Ax+b$) and same pairwise independence properties using only $m+n$ random bits.}
If one uses a random sparse Toeplitz matrix, the construction of hash family takes only $O(n)$ random bits and the hash operation can be faster because of the structure of the matrix.

\section{The \texttt{MB-WISH} algorithm and its analysis} \label{sec:wish}

Let us assume $\Omega = \ff_q^n$ where $q$ is a prime-power. Let us also fix an ordering among the elements of $\ff_q\equiv  \{\alpha_0, \alpha_1,\dots , \alpha_{q-1}\}$ and write $\alpha_0< \alpha_1 < \dots <\alpha_{q-1}$. In this section, the symbol `$<$' just signifies a fixed  ordering and has no real meaning over the finite field. Extending this notation, for any two vectors $x, y \in \ff_q^m$, we will say $x <y$ if and only if the $i$th coordinates of $x$ and $y$, satisfy $x_i<y_i$ for all $i=1, \dots, m$. Below $\mathbf{1}$ denotes an all-one vector of a dimension that would be clear from context. Also, for any event $\Ec$ let $\mathds{1}[\Ec]$ denote the indicator random variable for the event $\Ec$.

The MAX-oracle for \texttt{MB-WISH} performs the following optimization, given $A \in \ff_2^{m \times n}; b,s \in \ff_q^m$:
\vspace{-0.04in}
\begin{align}
\max_{\sigma \in \ff_q^n : A\sigma +b < s} w(\sigma).
\end{align}

The modified \texttt{WISH} algorithm is presented as Algorithm \ref{algo:mbwish}. The main result of this section is below.
\begin{algorithm}
\caption{\texttt{MB-WISH} algorithm for $\Omega = \ff_q^n$\label{algo:mbwish}}

\begin{algorithmic}
      \REQUIRE $r$, $\gamma= \frac{q}{3r}(\frac{1}{2}-\frac{r}{q})^{2}$, $\ell = \lceil \frac{1}{\gamma}\ln \frac{2n}{\delta} \rceil$, $n'=\lceil n \log_{q/r} q\rceil$
      \STATE $M_0 \equiv \max_{\sigma\in \ff_q^n} w(\sigma)$
       \FOR{$i\in\{1,2,\dots,n'\}$}
         \FOR{$k \in \{1,\dots,\ell\}$}
           \STATE Sample hash functions $h_i \equiv h_{A^{i},b^{i}}$ uniformly at random from $\Hc_{i,n}$ as defined in \eqref{eq:fam}
           \STATE $w_{i}^{(k)}=\max_{\sigma: A^{i}\sigma+b^{i} < \alpha_r \cdot \mathbf{1}} w(\sigma) $
        \ENDFOR
         \STATE $M_{i}={\rm Median}(w_{i}^{(1)},w_{i}^{(2)},\dots,w_{i}^{(\ell)})$
       \ENDFOR
       \STATE Return $M_{0}+(\frac{q}{r}-1)\sum_{i=0}^{n'-1}M_{i+1}\big(\frac{q}{r}\big)^{i} $
\end{algorithmic}
\end{algorithm}

\begin{theorem}\label{thm:mb}
Let $q$ be a power of prime and  $r=\lfloor{\frac{q-1}{2}}\rfloor$. Let $\Omega = \ff_q^n$. For any  $\delta> 0$ and a positive constant $\gamma \equiv \frac{q}{3r}(\frac{1}{2}-\frac{r}{q})^{2}$, Algorithm \ref{algo:mbwish} makes $\Theta(n \log \frac{n}{\delta})$ calls to the MAX-oracle and, with probability at least $1-\delta$ outputs a $(\frac{q}{r})^{2}$-approximation of $S_w(\Omega)$.
\end{theorem}

The theorem will be proved by a series of lemmas.
The key trick that we are using is to ask the MAX-oracle to solve an optimization problem over not a single bin, but {\em multiple bins} of the hash function. The hash family is defined in the following way. We have $h_{A,b}:\ff^n \to \ff^m: x \mapsto Ax+ b$, the operations are over $\ff_q$.
Let
\begin{align}\Hc_{m,n} = \{h_{A,b}: A \in \ff_q^{m \times n}, b \in \ff_q^m\}.\label{eq:fam}
\end{align}
The coding theoretic intuition behind our technique is following. The set of configurations $\{\sigma \in \ff_q^n : A\sigma ={\bf 0}\}$ forms a linear code of dimension $n-m$. The bins of the hash function define the cosets of this linear code. We would like to chose $q^r$ cosets of a random linear code and the find the optimum value of $w$ over the configurations of these cosets as the MAX-oracle.
To choose a hash function uniformly and randomly from $\Hc$, we can just choose the entries of $A$ and $b$ uniformly at random from $\ff_q$ independently.

Note that, the hash family $\Hc_{m,n}$ as defined in \eqref{eq:fam} is uniform and pairwise independent.

It  follows from the following more general result.

\begin{lemma}\label{lem:pair}
Let us define  $Z_{\sigma}$ to be the indicator random variable denoting $A\sigma+b < \alpha_r\cdot \mathbf{1}$ for some $r\in \{0, \dots, q-1\}$ and $A,b$ randomly and uniformly sampled from $\Hc_{m,n}$. Then $\Pr(Z_{\sigma}=1)=\big(\frac{r}{q}\big)^{m}$ and for any two configurations $\sigma_{1},\sigma_{2}\in \ff_q^n$ the random variables $Z_{\sigma_{1}}$ and $Z_{\sigma_{2}}$ are independent.
\end{lemma}

The proof is delegated to Appendix~\ref{sec:proof}.

Fix an ordering of the configurations $(\sigma_i,1\leq i \leq q^n)$ such that $1\leq j \leq q^n, w(\sigma_j)\geq w(\sigma_{j+1})$. We can also interpolate the space of configuration to make it continuous by the following technique. For any positive real number $x= z+f$, where $z =\lfloor x\rfloor $ is the integer part and $f = x-z$ is the fractional part, define $w(\sigma_x) = w(\sigma_z)$.  For $i \in \{0,1,2,\ldots, n'\equiv \lceil n \log_{q/r} q\rceil\}$, define $\beta_i=w(\sigma_{t^i}) = w(\sigma_{\lfloor t^i\rfloor})$, where $t =\frac{q}{r}$. We take $w(\sigma_k) =0$ for $k > q^n$.
To prove Theorem \ref{thm:mb} we need the following crucial lemma as well.
\begin{lemma}\label{lem:cru}
Let $M_{i}={\rm Median}(w_{i}^{(1)},\dots,w_{i}^{(\ell)})$ be defined as in the Algorithm \ref{algo:mbwish}. Then for $\gamma = \frac{q}{3r}(\frac{1}{2}-\frac{r}{q})^{2}$, we have,
$\Pr \bigg(M_i \in [\beta_{\min(i+1,n')},\beta_{\max(i-1,0)}] \bigg) \ge 1-2\exp(-\gamma \ell).$
\end{lemma}
\vspace{-0.1in}
The proof of this lemma is again delegated to Appendix~\ref{sec:proof} along with the full proof of Theorem \ref{thm:mb}.
Here we provide a sketch of the remaining steps.
From Lemma \ref{lem:cru}, the output of the algorithm  lies in the range $[L',U']$ with probability at least $1-\delta$ where
$
L'=\beta_0+(t-1)\sum_{i=0}^{n'-1}\beta_{\min\{i+2,n'\}}t^i$ and $U'=\beta_0+(t-1)\sum_{i=0}^{n'-1}\beta_it^i$. $L'$ and $U'$ are a factor of $t^2$ apart.
Now, following an argument similar to \eqref{eq:bound}, we can show $L' \le S_w(\Omega) \le U'.$

Therefore Algorithm \ref{algo:mbwish} provides a $t^2$-approximation to $S_\Omega(w)$. The total number of calls to the MAX-oracle is
$n'\ell +1 = O(n\log(n/\delta))$.

As an example of this result, suppose $q=3$. In this case the algorithm provides a $9$-approximation. Later, in the experimental section, we have used a ferromagnetic Potts model with $q=5$. The approximation guarantee that \texttt{MB-WISH} provides in that case is a factor of $\frac{25}{4}=6.25$. Note that, for a $5$-ary Potts model, it is only natural to use our algorithm instead of converting it to binary in conjunction with the original algorithm of Ermon et al.

Instead of pairwise independent hash families, if we employ $k$-wise independent families, it leads to a better decay probability of error. However it does not improve the approximation factor.

\paragraph{\texttt{MB-WISH} with unconstrained optimization oracle.} We can modify and generalize the results of Achlioptas and Jiang \cite{achlioptas2015stochastic} to formulate a version of \texttt{MB-WISH} that can use unconstrained optimizers as the MAX-oracle. The MAX-oracle for this algorithm performs an unconstrained optimization of the  form:  $
\max_{\sigma \in B}w(A\sigma +b)$, given $A\in \ff_q^{m\times n}, b \in \ff_q^n$ and a set $B \subseteq \ff_q^m$.

The aim is to carefully design $B$ so that all the desirable statistical properties are satisfied.
This part is quite different from the hashing-based analysis  and not an immediate extension of \cite{achlioptas2015stochastic}. We provide the algorithm (\texttt{Unconstrained MB-WISH}) and its analysis in Appendix~\ref{sec:mba}.

\section{\texttt{MB-WISH} for computing permanent}\label{sec:perm}
Recall the permanent of a matrix as defined in Eq.~\eqref{eq:perm}: ${\rm Perm}(A) \equiv \sum_{\sigma \in S_n} \prod_{i=1}^n A_{i,\sigma(i)}$. We will show that it is possible to  approximate the permanent with a modification of the \texttt{MB-WISH} algorithm and our idea of using multiple bins for optimization in the calls to MAX-oracle.
Also, recall from Section~\ref{sec:wish} that we set $\ff_q\equiv  \{\alpha_0, \alpha_1,\dots , \alpha_{q-1}\}$ where there exists a fixed ordering among the elements. We set $q \ge n$ and consider any $\sigma \in S_n$ as an $n$-length vector  over $\ff_q$ (that is by identifying $1, 2, \dots, n$ as $\alpha_0, \alpha_1, \dots, \alpha_{n-1}$ respectively). Then we define a modified hash family
$\Hc_{m,n} = \{h_{A,b}:  A \in \ff_q^{m \times n}, b \in \ff_q^m \}$ with $h_{A,b}: S_n \to \ff_q^m: \sigma \mapsto A\sigma + b,$ the operations are over $\ff_q$.

However, when calling the MAX-oracle, we need to make sure that we are getting a permutation as the output. Hence the modified MAX-oracle for computing permanent will be:
\begin{align}
\max_{\sigma \in \ff_q^n} w(\sigma),
\text{ subject to, }
 A\sigma +b < \alpha_r\cdot \mathbf{1} ;
 \sigma < \alpha_{n-1}\cdot \mathbf{1} ;
\sigma(i) \ne \sigma(j) \forall i \ne j.
\end{align}
These constraints, which are all linear, ensures that the MAX-oracle returns a permutation over $n$ elements.
With this change we propose Algorithm~\ref{algo:permanent} to compute permanent of a matrix and
call it \texttt{PERM-WISH}.  The full algorithm is provided in Appendix~\ref{sec:perm_}.

The main result of this section is the following.
\begin{theorem}\label{thm:perm}
Let $D$ be any $n\times n$ matrix. Let $q>n$ be a power of prime and  $r=\lfloor{\frac{q-1}{2}}\rfloor$.  For any  $\delta> 0$ and a positive constant $\gamma = \frac{q}{3r}(\frac{1}{2}-\frac{r}{q})^{2}$, Algorithm \ref{algo:permanent} makes $\Theta(n \log \frac{n}{\delta})$ calls to the MAX-oracle and, with probability at least $1-\delta$ outputs a $(\frac{q}{r})^{2}=(4+O(1/n))$-approximation of $\mathrm{Perm}(D)$.
\end{theorem}

The proof of Theorem~\ref{thm:perm} follows the same trajectory as in Theorem~\ref{thm:mb}. The constraints in MAX-oracle ensures that a permutation is always returned. So in the proof of Theorem~\ref{thm:mb}, the $w_{i}^{(k)}$s can be though of as permutations instead in this setting.
It should be noted that, we must take $q>n$ for \texttt{PERM-WISH} to work. That is the reason we get a $(4+O(1/n))$-approximation for the permanent.

It also has to be noted that, since $q$ is large, the straightforward extension of \texttt{WISH} algorithm would have provided only a $q^2 = n^2$-approximation of the permanent. Therefore the idea of using optimizations with multiple bins  are crucial here as it lead to a close to $4$-approximation.

\begin{remark}[Constraints]
It turns out that the constraints in the MAX-oracle in Algorithm \ref{algo:permanent} are linear/affine. Therefore they are still easy to implement in different CSP softwares.
\end{remark}

\section{Derandomization with structured hash families}\label{sec:der}
In this section, we show that it is possible to construct pairwise independent hash family $\{\ff_2^n\to \ff_2^m\}$ using only $O(n)$ random bits such that any hash function from the family still has the structure $h(x) = Ax+b$. Both of our constructions can be easily extended to $q$-ary alphabets. The proofs of the  propositions below can all be found in Appendix~\ref{app:der} in the supplementary material.

\noindent{\bf Construction 1:}
Let $f(x) \in \ff_2[x]$ be an irreducible polynomial of degree $n$. We construct the finite field $\ff_{2^n}$ with the $\zeta$, root of $f(x)$ as a generator of $\ff_{2^n}^\ast$. Now, any $x \in \ff_2^n$ can be written as a power of $\zeta$ via a natural map $\phi: \ff_2^n \to \ff_{2^n}$. Indeed, for any element $\zeta^{k} \in \ff_{2^n}^\ast$ consider the polynomial $\zeta^k \mod f(\zeta)$ of degree $n-1$. The coefficients of this polynomial from an element of $\ff_2^n$. $\phi$ is just the inverse of this map. Also, assume that the all-zero vector is mapped to $0$ under $\phi$.

Let $x\in \ff_2^n$ be the configuration to be hashed. Suppose the hash function is $h_{\nu,b}$, indexed by $\nu\in \ff_2^n$ and $b \in \ff_2^m$. The hash function is defined as follows:

Let  $\nu \in \ff_{2}^n$.
 Compute $z=\phi^{-1}(\phi(x) \cdot \phi(\nu) \mod f(\zeta)) \in \ff_2^n.$
 Let $y\in \ff_2^m$ be the first $m$ bits of  $z$.
  Finally, output $y+b$, where $b \in \ff_2^m$.

\begin{proposition}\label{prop:aff}
The hash function $h_{\nu,b}$ can be written as an affine transform ($x \mapsto Ax+b$) over $\ff_2^n$.
\end{proposition}

Note that, to chose a random and uniform  hash function from $\{h_{\nu,b}, \nu \in \ff_2^n, b \in \ff_2^m\}$, one needs $m+n$ random bits.

It follows that the  hash family is  pairwise independent.

\begin{proposition}\label{prop:pair}
The hash family $\{h_{\nu,b}, \nu \in \ff_2^n, b \in \ff_2^m\}$ is uniform and pairwise independent.
\end{proposition}
Moreover the randomness used to construct this hash function is also optimal. It can be shown that,

the size of a pairwise independent hash family $\{h:\{0,1\}^n \to \{0,1\}^m\}$ is at least $2^{m+n}-2^{n}+1$  (see,  \cite{stinson1996connections}).

This implies that $m+n$ random bits were essential for the construction.

\noindent{\bf Construction 2: Toeplitz matrix.}
In \cite{ermon2013optimization}, a Toeplitz matrix was used as the hash function. In a Toeplitz matrix, each descending diagonal from left to right is fixed, i.e., if $A_{i,j}$ is the $(i,j)$th entry of a Toeplitz matrix, then $A_{i,j} = A_{i-1,j-1}$. So to specify an $m \times n$ Toeplitz matrix one needs to provide only $m+n-1$ entries (entries of the first row and first column).
Consider the  random $m \times n$ Toeplitz matrix $A_T$ where each of the entries of the first row and first column are chosen with equal probability from $\{0,1\}$, i.e., each entry in the first row and column is a Bernoulli($0.5$) random variable. The random hash function $h_{A_T,b}: x\mapsto A_T x+ b$, is constructed by chosing a uniformly random $b\in \ff_2^m.$

\begin{proposition}\label{prop:top}
The hash family $\{h_{A_T, b}\}$  is uniform and pairwise independent \cite{ermon2013optimization}.
\end{proposition}
Note that, since the number of random bits required from this construction is $2m+n-1$, in terms of memory requirement this hash family is slightly suboptimal to the previous construction. However, Toeplitz matrix allow for much faster computation of the hash function (matrix-vector multiplication with Toeplitz matrix takes only $O(n \log n)$ time compared to $\Omega(mn)$ for unstructured matrices).

 We remark that {\em sparse Toeplitz Matrices} also can be used as our hash family, further reducing the randomness.

In particular, we could  construct a Toeplitz matrix with Bernoulli($p$) entries for $p <0.5$.  While the pairwise independence of the hash family is lost, it is still possible to analyze the \texttt{MB-WISH} algorithm with this family of hashes since they form a {\em strongly universal} family \cite{stinson1996connections}.  The number of random bits used in this hash family is $(m+n-1)h(p)+m$.  This construction allows us to have sparse rows in the Toeplitz matrix for small values of $p$, which can lead to further speed-up.

Both the constructions of this section extend to $q$-ary alphabet straightforwardly.

\section{Experimental results} \label{sec:exp}
All the experiments were performed in a shared parallel computing environment that is equipped with 50 compute nodes with 28 cores Xeon E5-2680 v4 2.40GHz processors with 128GB RAM.

\subsection{Experiments on simulated Potts model}

We implemented our algorithm to estimate the partition function of  Potts Model (a generalization of Ising Model). Recall that the partition function of the  Potts model on a graph $G=(V,E)$   is given in Eq.~\eqref{eq:part}.

For our simulation, we have randomly generated the graph $G$ with number of nodes $n \equiv |V|$ varying in $4,5,6,7,8,9,10,11,12$ and corresponding regular degree $d=2,2,4,4,4,4,6,6,8$ using a python library \texttt{networkx}. We took the number of states of the Potts model $q=5$, the external force $H$ and the spin-coupling $J$   to be 0.1 and then varied the values of $\zeta$.
The partition functions for different cases are calculated using both brute force technique and our algorithm. We have used a python module \texttt{constraint} (for Constraint Satisfaction Problems (CSPs)) to handle the constrained optimization for MAX-oracle. The obtained approximation factors for different values of $\zeta$ are listed in Table~\ref{tab:part}.

For $n=10, 11, d=6$ the approximation factor for \texttt{MB-WISH} is exactly 1 (up to the precision of the number system used). However the time taken by \texttt{MB-WISH} is much higher.
For $n=12, d=8$, \texttt{MB-WISH} gives an approximation factor of $2.5$ after running for eight hours in a parallel computing environment that is equipped with 50 compute nodes with 28 cores Xeon E5-2680 v4 2.40GHz processors with 128GB RAM.

\begin{table*}[t]
\begin{center}
	\begin{tabular}{| c| c |c |c | c |c|c|  }
		\hline
$\zeta $ & $n=4, d=2$ & $n=5,d=2$ & $n=6,d=4$ & $n=7,d=4$ & $n=8,d=4$ & $n=9,d=4$ \\
\hline
 0 & 0.976 & 1.220 & 0.610 & 1.907 &  0.953 & 1.192\\
 \hline
 5& 0.580 & 0.708 & 1.639  &  0.755 &  0.630 & 0.599\\
\hline
 10 & 0.7470 & 1.191 & 3.271 & 0.989 & 1.875 &  1.25\\
\hline
 15 &   1.430 & 1.036 & 1.013 & 1.224 & 1.399 & 1.692\\
 \hline
 20 &  1.032 & 1.590 & 1.141 & 1.173 & 1.365 & 1.491\\
 \hline
 25 & 0.839 & 1.118 & 1.339 & 1.035 & 1.429 & 1.326\\
 \hline
 30 & 0.510 & 4.0562 & 2.226 & 1.060 & 0.690 & 2.122 \\
 \hline
 35 & 1.073 & 5.442 & 0.489 & 2.871 & 1.639 & 1.263\\
 \hline
 40 &  1.210 &  2.434 &  0.980 &  0.582 &  0.666 & 0.969 \\
 \hline
 45 & 1.127 &  4.640 &  2.348 &  1.336&  0.3673 & 1.341\\
 \hline
 50 &  1.152 &  1.025 &  2.511 &  3.4307 &  1.1522 &2.636 \\
\hline
	\end{tabular}
\caption{The ratio of the partition function calculated by \texttt{MB-WISH} and the actual value calculated by brute force: $\frac{\hat{Z}}{Z}$. \label{tab:part}}
\end{center}
\end{table*}

For graphs with larger number of vertices, it is difficult to compute the partition function of Potts Model by brute force. Therefore, we compare the partition function computed by  \texttt{MB-WISH} ($\hat{Z}$) with the one ($Z$) computed by a belief propagation algorithm in the in the PGMPY library in \texttt{python} \cite{PGMPY}.

Again, for our simulation, we have randomly generated the graph $G$ with number of nodes $n \equiv |V|$ varying in $10,15,20,25,30,40,50,$ and with regular degree $d=4$ using a python library \texttt{networkx}.
 We took the number of states of the Potts model $q = 5$, the external force $H$ and the spin-coupling   to be   $ 0.1$ and then varied the values of  $\zeta $. In our experiments each optimization  instances are run with a timeout of $3,5,10,15 ,20,20,25 $ minutes for $n=10,15,20,25,30,40,50$ respectively. The results are summarized in Table \ref{tab:part1}.
\begin{table*}[h!]
\begin{center}
	\begin{tabular}{|c|c|c|c|}
		\hline
			&$\zeta=1$ & $\zeta=2$ &   $\zeta=5$ 	\\
			\hline
		$|V|=10$ & $0.7027854$   & $0.537426$     &$0.607568$      \\
		\hline
		$|V|=15$ & $0.843724$   & $0.893989$    &$0.669195$     \\
		\hline
	$|V|=20$ &  $1.0117$  &  $1.021$   & $0.819284$      \\
		\hline
		$|V|=25$ & $0.60450$   & $0.967849$     &$0.602233$      \\
		\hline
		$|V|=30$ & $0.72632$   & $0.65447$    &$1.12595$     \\
		\hline
		$|V|=40$ & $0.83407$   & $0.69569$    &$1.1989$     \\
		\hline
		$|V|=50$ & $0.760586$   & $0.80906$    &$2.1912$     \\
		\hline
	\end{tabular}
\end{center}
\caption{The ratio of the partition function calculated by \texttt{MB-WISH} and the approximate value calculated by belief propagation: $\frac{\hat{Z}}{Z}$. \label{tab:part1}}
\end{table*}

\subsection{Real-world CSPs}
 Many  instances of real-world graphical models are available in

\url{http://www.cs.huji.ac.il/project/PASCAL/showExample.php}.
Notably, some of them (e.g., image alignment, protein folding) are defined on non-Boolean domains, which justify the use of \texttt{MB-WISH}. We have computed the partition functions for several datasets from this repository.

The  dataset \texttt{Network.uai}
is a Markov network with $120$ nodes each having a binary value. A configuration here is a binary sequence of length $120$. To calculate the partition function, we need to find the sum of weights for $2^{120}$ different configurations. In order to use \texttt{Unconstrained MB-WISH}, we view each configuration as a $16$-ary string of length $30$.  Our results for the log-partition came out to be $156.00$ with one hour time out for each call to the MAX-oracle. The benchmark for the log-partition function is provided to be $163.204$.

The  \texttt{Object detection} dataset comprised of $60$ nodes each having a $11$-ary value and by \texttt{Unconstrained MB-WISH} we found the log-partition function to be $-38.9334$.  The \texttt{CSP} dataset is a Markov network with $30$ node having a ternary value: we found the log partition function to be $-39.9933$.  For these datasets there were no baselines.

\subsection{Permanent}
We demonstrate the \texttt{PERM-WISH} algorithm by calculating the permanent of randomly generated matrix of size $n \times n$ where $n=5,6,7,8,9,10$ and  compute the ratio of the value achieved by our method ($\hat{Z}$) and actual permanent ($Z$) by brute force. For the purpose of constrained optimization, we use the python modules described above. For the experiment we take $q=11$ and $r=5$. We find $\hat{Z}/Z = 2.137, 1.103, 1.245, 2.407, 0.730, 1.527$ for $n =5,6,7,8,9,10$ respectively.

\subsection{Experiments to estimate Total Variation distance}
 We use \texttt{Unconstrained MB-WISH} to compute total variation distances between two high dimensional (up to dimension $100$) probability distributions, generated via an iid model and two different Markov chains (domain size up to $5^{100}$). There is no baseline to compare with for the Markov chain case, but for the product distributions the method shows good accuracy. The detailed results are provided below.

 The total variation (TV) distance between any two discrete distributions $P$ and $Q$ with common sample space $\Pc$ is defined to be
$$
\|P-Q\|_{TV} = \sup_{A\subseteq \Pc}|P(A) -Q(A)| = \frac12 \sum_{\sigma \in \Pc}|P(\sigma) -Q(\sigma)|.
$$
It is difficult to estimate the TV distance between two distribution over product spaces. Even if the random variables are independent, there is no easy way to analytically compute the total variation distance. Simply consider finding TV distance between joint distributions of $n$ random variables that can take value in $\{0, 1, \dots, q-1\}$. In that case, we seek to find,
$$
\frac12 \sum_{\sigma \in \{0,1, \dots, q-1\}^n}|P^{n}(\sigma) -Q^{n}(\sigma)|,
$$
which is in the exact form of Eq.~\eqref{eq:main}.  Therefore we can use \texttt{MB-WISH} algorithm to estimate the total variation distance. To see how good the performance is, we compare the computed value with the known upper and lower bound on TV distance based on Hellinger distance: $h(P,Q)^2 \equiv \sum_{\sigma \in \Pc}(\sqrt{P(\sigma)}-\sqrt{Q(\sigma)})^2$. It is known that\footnote{See \url{https://stanford.edu/class/stats311/Lectures/full_notes.pdf}.},
$$
\frac12 h(P,Q)^2\leq \|P-Q\|_{TV}\leq h(P,Q)\sqrt{1- h(P,Q)^2/4}.
$$

For the experiments, we  choose two distribution over $5$ points ($q=5$): $P=[0.2,0.2,0.2,0.2,0.2]$   and $Q=[0.2,0.2+\epsilon,0.2-\epsilon,0.2+\epsilon,0.2-\epsilon]$, where $\epsilon=10^{-2}$. The distributions are chosen very close to each other. Here the distribution $P^n$ and $Q^n$ are supported on $\{0,1,2,3,4\}^n$ where $n$ can be any natural number.
It is known that,
\begin{align*}
 h(P^n,Q^n)^2 &=2-2\prod_{i=1}^{n}(1-\frac{1}{2}h(P_i,Q_i)^2).
\end{align*}

To compare our method with the upper and lower bound with Hellinger distance we choose the random variables to be independent, although that our method is capable of approximation irrespective of any such conditions.
We do our experiments in a time constrained manner (10 minute for each calls to MAX-oracle).
The results are summarized in Table~\ref{tab:hel}.

We also ran the experiment when the joint distributions are not of independent random variables and therefore no standard upper or lower bounds are available. We assume the random variables $X_1,X_2,\ldots,X_n $ form a Markov chain over $5$ states. We consider two different first order Markov chains with initial probability distributions $P=[0.2,0.2,.02,0.2,0.2 ] $ and $ Q=[0.2,0.21,0.19,0.19,0.21]$ respectively and the state-transition matrices given by:
\[
T=
\begin{bmatrix}
0.2204 & 0.2813 & 0.2285  & 0.02819 & 0.24155 \\
0.11137 & 0.1277 & 0.04995 & 0.3949 & 0.3159 \\
 0.01068 & 0.3220 & 0.2.387 & 0.11459 & 0.3139 \\
0.3714 & 0.3397 & 0.1599 & 0.02134 & 0.1074 \\
0.02691 & 0.38807 & 0.2.869 & 0.05981 & 0.2382
\end{bmatrix}
\]
and $T+D$ where every row of $D$ is given by the vector  $[0,10^{-4},-10^{-4},-10^{-4},10^{-4}].$
We estimate the total variation distance between $P(X_1,X_2,\ldots X_n)$  and $Q(X_1,X_2,\ldots X_n)$ in Table~\ref{tab:Markov} below.

\begin{table*}
\begin{center}

	\begin{tabular}{|c|c|c|c|}
		\hline
			& Lower Bound & \texttt{MB-WISH} & Upper Bound 	\\
			\hline
		$n=4$ & $0.001$   & $0.023$     &$0.0316$      \\
		\hline
		$n=6$ & $0.0015$   & $0.025$    &$0.0387$     \\
		\hline
	$n=8$ &  $0.002$  &  $0.0255$   & $0.0447$      \\
		\hline
	$n=10$&    $0.0025$& $0.045$    & $0.05$      \\
	\hline
	$n=12$ & $0.003$   &$0.0186$     &$0.0547$     \\
	\hline
	$n=30$ & $0.007492$   & $0.1097$     &$0.086557$      \\
		\hline
		$n=40$ & $0.00998$   & $0.1147$    &$0.099917$     \\
		\hline
	$n=50$ &  $0.012471$  &  $0.07265$   & $0.1116$      \\
		\hline
	$n=100$&    $0.02486$& $0.1767$    & $0.15768$      \\
	\hline
	\end{tabular}
\end{center}
\caption{Computation of TV distance of product distribution by  \texttt{MB-WISH} compared with upper bound (UB) and lower bound (LB) based on Hellinger distance. Here $\epsilon = 10^{-2}$.\label{tab:hel}}
\end{table*}

 \begin{table}
\begin{center}
	\begin{tabular}{|c|c|c|c|}
		\hline
			& \texttt{Unconstrained MB-WISH} 	\\
			\hline
		$n=30$ &   $0.031307$    \\
		\hline
		$n=40$ & $0.03622$    \\
		\hline
	$n=50$ &     $0.045950$   \\

	\hline
	\end{tabular}
\end{center}
\caption{Computation of TV distance for Markov case  by  \texttt{Unconstrained MB-WISH}. Each call to MAX-oracle is time limited to 40 minutes. \label{tab:Markov}}
\end{table}

\bibliographystyle{plain}

\bibliography{references}

\appendix

\section{Omitted proofs from Section~\ref{sec:wish}}\label{sec:proof}

\begin{proof}[Proof of Lemma~\ref{lem:pair}]
Let $A_i$ denote the $i$th row of $i$ and $b_i$ denote the $i$th entry of $b$. Then
$\mathds{1}[A\sigma+b < \alpha_r\cdot \mathbf{1}]=\bigwedge_{i=1}^{m}\mathds{1}[A_{i}\sigma+b_i <\alpha_r]$.

For all configurations $\sigma \in \Omega, \forall i$, we must have
\begin{align*}
\Pr(A_i\sigma +b_i < \alpha_r) =\sum_{j=0}^{r-1} \Pr(A_i\sigma +b_i =\alpha_j)=\frac{r}{q}.
\end{align*}

As  $A_i,b_i$  are independent $1\le i \le m$, we must have that $\Pr(A\sigma +b < \alpha_r\cdot \mathbf{1})=(\frac{r}{q})^m$. Now for any two configurations $\sigma_1 ,\sigma_2\in \ff_q^n$,
\begin{align*}
& \Pr(A_i\sigma_1+b_i < \alpha_r \wedge A_i\sigma_2+b_i < \alpha_r ) \\
&=\sum_{k=0}^{r-1}\sum_{j=0}^{r-1}\Pr(A_i\sigma_1+b_i =\alpha_k \wedge A_i\sigma_2+b_i = \alpha_j )\\
&=\sum_{k=0}^{r-1}\sum_{j=0}^{r-1}\Pr(A_i\sigma_1+b_i =\alpha_k)\\
& \qquad \cdot\Pr(A_i\sigma_2+b_i = \alpha_j|A_i\sigma_1+b_i =\alpha_k)\\
&=\sum_{k=0}^{r-1}\sum_{j=0}^{r-1}\Pr(A_i\sigma_1+b_i =\alpha_k)\Pr(A_i(\sigma_2-\sigma_1) = \alpha_j-\alpha_k)\\
&=r^2(1/q)(1/q)=(r/q)^2.
\end{align*}
As all the rows are independent, $\Pr(A\sigma_1 +b < \alpha_r\cdot \mathds{1} \wedge A\sigma_2 +b < \alpha_r\cdot \mathds{1})=(\frac{r}{q})^{2m}.$
\end{proof}

\begin{proof}[Proof of Lemma~\ref{lem:cru}]
 Consider the set of $\lfloor t^j\rfloor$ heaviest configuration
\begin{align*}
\Omega_j=\{\sigma_1,\ldots \sigma_{\lfloor t^j \rfloor} \} .
\end{align*}
Let $S_j(h_i) = |\{\sigma \in \Omega_j: A^i\sigma+b^i < \alpha_r\cdot\mathbf{1}\}|.$
By the uniformity property of the hash function,
\begin{align*}
\avg S_j(h_i)=\avg \sum_{\sigma \in \Omega_j}\mathds{1}[h_i(\sigma) < \alpha_r\cdot \mathbf{1}]=\frac{\lfloor t^j\rfloor}{t^i}.
\end{align*}
For each configuration $\sigma$ let us denote the random variable $\bar{Z^i_{\sigma}}=\mathds{1}[\{h_i(\sigma) <  \alpha_r\cdot \mathbf{1} \}]-\frac{1}{t^i}$. By our design $\avg \bar{Z^i_{\sigma}}=0.$ Note that, $S_j(h_i) - \avg S_j(h_i) = \sum_{\sigma\in \Omega_j}\bar{Z^i_{\sigma}}$.
Also, from Lemma \ref{lem:pair}, the random variables $\bar{Z^i_{\sigma}}$s are pairwise independent.
Therefore,
\begin{align*}
\Var S_j(h_i) =\Var(\sum_{\sigma\in \Omega_j}\bar{Z^i_{\sigma}})=\sum_{\sigma\in \Omega_j}\avg \bar{Z^i_{\sigma}}^2=\frac{\lfloor t^j\rfloor}{t^i}(1-\frac{1}{t^i}).
\end{align*}
Now,
for any $1 \le k \le \ell$,
\begin{align*}
\Pr(w_i^{(k)}\geq \beta_j)=\Pr(w_i^{(t)} \geq w(\sigma_{\lfloor t^j\rfloor}))\geq \Pr(S_j(h_i)\ge1)
 \geq 1- \Pr(S_j(h_i)\le 0).
\end{align*}
Let $j =i+1$. Then, using Chebyshev inequality,
\begin{align*}
\Pr(S_j(h_i)\le 0) = \Pr\Big(S_j(h_i) - \avg S_j(h_i) \le - \frac{\lfloor t^j\rfloor}{t^i}\Big)
\le  \frac{\Var S_j(h_i)}{(\frac{\lfloor t^j\rfloor}{t^i})^2} \le \frac{t^{i}(1-1/t^i)}{\lfloor t^{j}\rfloor} < \frac{t^i-1}{t^{i+1}-1} \le \frac{1}{t} = \frac{r}{q}.
\end{align*}

Therefore,
\begin{align*}
\Pr(w_i^{(k)}\geq \beta_{i+1}) \ge 1 - \frac{r}{q}.
\end{align*}

Also,
$
\Pr(w_i^{(k)}\leq \beta_{i-1}) = \Pr(w_i^{(k)} \leq w(\sigma_{\lfloor t^{i-1}\rfloor})) \ge \Pr(S_{i-1}(h_i) =0).
$
Notice that the last inequality is satisfied because $S_{i-1}(h_i) =0$ implies  $w_i^{(k)} \leq w(\sigma_{\lfloor t^{i-1}\rfloor})$.
Now, continuing the chain of inequalities, using Markov inequality,
\begin{align*}
\Pr(w_i^{(k)}\leq \beta_{i-1})  \ge 1 - \Pr(S_{i-1}(h_i) \ge 1) \ge 1 - \avg S_{i-1}(h_i)
= 1-\frac{\lfloor t^{i-1}\rfloor}{t^i}\ge 1-\frac1t = 1-\frac{r}{q}.
\end{align*}

Now just by using Chernoff bound\footnote{We use the fact that if $X$ is a sum of iid $\{0,1\}$ random variables then $\Pr(X \ge \avg X(1+\delta)) \le \exp(-\avg X \delta^2/3)$.},
\begin{align*}
\Pr(M_i \le \beta_{i+1})\le \exp\Big(-\frac{\ell q}{3r}\Big(\frac12-\frac{r}{q}\Big)^2 \Big),

\quad \text{ and } \quad
\Pr(M_i \ge \beta_{i-1})\le \exp\Big(-\frac{\ell q}{3r}\Big(\frac12-\frac{r}{q}\Big)^2 \Big).
\end{align*}
This proves the lemma.

\end{proof}

\vspace{0.1in}

\begin{proof}[Proof of Theorem \ref{thm:mb}]
From Lemma \ref{lem:cru}, we have,
$\Pr \bigg[ \bigcap_{i=1}^{n} M_i \in [\beta_{\min(i+1,n')},\beta_{\max(i-1,0)}] \bigg] \ge 1-2n\exp(-\gamma \ell) = 1-\delta$
 for $\ell = \frac{1}{\gamma}\ln\frac{2n}{\delta}$ and by definition $M_{0}=\beta_{0}.$

The algorithm outputs $M_0+(t-1)\sum_{i=0}^{n'-1}M_{i+1}t^i$ which lies in the range $[L',U']$ with probability at least $1-\delta$ where
\begin{align*}
L'=\beta_0+(t-1)\sum_{i=0}^{n'-1}\beta_{\min\{i+2,n'\}}t^i  \quad \text{ and } \quad  U'=\beta_0+(t-1)\sum_{i=0}^{n'-1}\beta_it^i.\end{align*}

Now notice that, as $ \beta_0 \ge \beta_1$, we have
\begin{align*}
U' &=\beta_0+(t-1)\sum_{i=0}^{n'-1}\beta_it^i
\\& = \beta_0 + (t-1) (\beta_0 +\beta_1t)+\sum_{i=2}^{n'-1}\beta_it^i \\
& \leq t^2\beta_0+ t^2.(t-1). \sum_{i=2}^{n'-1}\beta_it^{i-2}\\
&\leq t^2(\beta_0+(t-1)\sum_{i=0}^{n'-1}\beta_{\min\{i+2,n'\}}t^i)=t^2L'.
\end{align*}

The only thing that remains to be proved is that $L' \le S_w(\Omega) \le U'.$ However that is true, by just following an argument similar to \eqref{eq:bound}. Indeed,
\begin{align*}
\sum_{i=0}^{n'-1}\beta_{i+1}(t^{i+1} -  t^{i}) \le
S_\Omega(w) \le\sum_{i=0}^{n'-1}\beta_{i}( t^{i+1} - t^{i} )\end{align*}
which implies $L' \le S_\Omega(w) \le U'.$

Therefore Algorithm \ref{algo:mbwish} provides a $t^2$-approximation to $S_\Omega(w)$. The total number of calls to the MAX-oracle is
$n'\ell +1 = O(n\log(n/\delta))$.
\end{proof}

\section{\texttt{MB-WISH} with unconstrained optimization oracle}
\label{sec:mba}
In this section, we provide an algorithm that uses unconstrained optimizations for the oracle, as in the case of Achlioptas and Jiang \cite{achlioptas2015stochastic}. We call this algorithm \texttt{Unconstrained MB-WISH}.

Let us assume $\Omega = \ff_q^n$ where $q$ is a prime-power. As before,  let us also fix an ordering among the elements of $\ff_q\equiv  \{\alpha_0, \alpha_1,\dots , \alpha_{q-1}\}$ and write $\alpha_0< \alpha_1 < \dots <\alpha_{q-1}$. Recall that, here the symbol `$<$' just signifies a fixed  ordering and has no real meaning over the finite field.

The MAX-oracle for \texttt{Unconstrained MB-WISH} performs an unconstrained optimization of the following form, given $A\in \ff_q^{m\times n}, b \in \ff_q^n$ and a set $B \subseteq \ff_q^m$:
\begin{align}
\label{eq:umb}
\max_{\sigma \in B}w(A\sigma +b).
\end{align}

The \texttt{Unconstrained MB-WISH} algorithm is presented as Algorithm~\ref{algo:mbawish}. The main result of this section is the following.
\begin{theorem}\label{thm:mba}
Let $q$ be a power of prime and  $r=\lfloor{\frac{q-1}{2}}\rfloor$. Let $\Omega = \ff_q^n$. For any  $\delta> 0$ and a positive constant $\gamma = \frac{q}{3r}(\frac{1}{2}-\frac{r}{q})^{2}$, Algorithm \ref{algo:mbawish} makes $\Theta(n \log \frac{n}{\delta})$ calls to the MAX-oracle (cf.~\eqref{eq:umb}) and, with probability at least $1-\delta$ outputs a $(\frac{q}{r})^{2}$-approximation of $S_w(\Omega)$.
\end{theorem}

\begin{algorithm}
\caption{\texttt{Unconstrained MB-WISH} algorithm for $\Omega = \ff_q^n$\label{algo:mbawish}}

\begin{algorithmic}
      \REQUIRE $\ell \rightarrow \lceil \frac{1}{\gamma}\ln \frac{2n}{\delta} \rceil$, $r,n'=\lceil n \log_{q/r} q\rceil$
      \STATE $M_0 \equiv \max_{\sigma\in \ff_q^n} w(\sigma)$
       \FOR{$i\in\{1,2,\dots,n\}$}
         \FOR{$k \in \{1,\dots,\ell\}$}
           \STATE Sample $n$ linearly independent vectors spanning $\ff_{q}^{n}$ and construct matrices $A$ and $R$ by taking the first $n-i$ columns and the last $i$ columns respectively
           \STATE $w_{i}^{(k)}=\max_{\substack{x \in \ff_{q}^{n-i} \\ y \in \{\alpha_0,\alpha_1,\dots,\alpha_{r-1}\}^{i}}} w(Ax+Ry+b)$
        \ENDFOR
         \STATE $M_{i}={\rm Median}(w_{i}^{(1)},w_{i}^{(2)},\dots,w_{i}^{(\ell)})$
       \ENDFOR
       \FOR{$i\in\{n+1,\dots,n'\}$}
         \FOR{$k \in \{1,\dots,\ell\}$}
           \STATE Sample full rank matrix $A \in \ff_{q}^{n \times n},b \in \ff_{q}^{n}$ uniformly at random. Set $\mathcal{S}_{i}$ as defined in Equation \ref{def:si}
           \STATE $w_{i}^{(k)}=\max_{y \in \mathcal{S}_i}w(Ay+b)$
        \ENDFOR
         \STATE $M_{i}={\rm Median}(w_{i}^{(1)},w_{i}^{(2)},\dots,w_{i}^{(\ell)})$
       \ENDFOR
       \STATE Return $M_{0}+(\frac{q}{r}-1)\sum_{i=0}^{n'-1}M_{i+1}\big(\frac{q}{r}\big)^{i} $
\end{algorithmic}
\end{algorithm}

To prove this theorem we borrow some ideas from coding theory. We define a linear $q$-ary code $C$ of dimension $n-m$ and length $n$ as the set of vectors $\{Ax: x\in \ff_q^{n-m}\}$ where $A$ is a full-rank matrix of size $n \times n-m$ and rank $n-m$. For a vector $a \in \ff_{q}^{n}$, we define the set $\{a+C\}$ as a coset of $C$. It is well known that  $\ff_{q}^{n}$ is partitioned by the $q^{m}$ distinct  cosets, each of size $q^{n-m}$. The main technique behind our algorithm is that for a random linear code $C$ of size $q^{n-m}$, we randomly sample $r^{m}$ distinct cosets of $C$. Subsequently, we find the maximum value $w(x)$ of an element among those $r^{m}$ cosets.

Let $E\equiv \{e_i\}_{i=1}^{n}$ be a set of $n$ linearly independent vectors in $\ff_{q}^{n}$ chosen randomly and uniformly (one can uniformly sample the columns of $E$ and resample if it belongs to the span of the already chosen columns). Let $A$ denote the matrix formed by the first $n-m$ vectors of $E$ as columns and let $R$ be the matrix formed by the remaining $m$ vectors as column. Also let $b$ be a vector sampled randomly and uniformly from $\ff_{q}^{n}$.
The MAX-oracle for \texttt{Unconstrained MB-WISH} is going to perform the following optimization   when $m \le n$:
\begin{align}
\max_{\sigma_1 \in \ff_q^{n-m}, \sigma_2 \in  \{\alpha_0,\alpha_1,\dots,\alpha_{r-1}\}^{n \times m}} w(A\sigma_1+R\sigma_2 +b ).
\end{align}

Analogous to Theorem~\ref{thm:mb}, here we are creating union of $r^m$ distinct random bins. If we can prove that, for any element of $\ff_q^n$, the probability that it belongs to one of these bins is $(\frac{r}{q})^m$ and for any pair of different elements from $\ff_q^n$, whether they belong to one of these bins are independent (pairwise independence), the rest of the proof of Theorem~\ref{thm:mba} will  just follow that of Theorem~\ref{thm:mb}.

In particular, we just have to prove the lemma that is analogous to Lemma~\ref{lem:pair}.
Define a set
\begin{align*}
S_{A,R,b}\equiv \{Ax+b+Ry \mid x \in \ff_q^{n-m},y \in \{\alpha_0,\alpha_1,\dots,\alpha_{r-1}\}^{m}\}.
\end{align*}
For each configuration $\sigma \in \ff_{q}^{n}$, associate an indicator random variable $Z_{\sigma}$ denoting whether $\sigma \in S_{A,R,b}$ .

\begin{lemma}\label{lem:ach1}
For each configuration $\sigma \in \ff_{q}^{n}$, we must have $\Pr(Z_{\sigma}=1)=\Big(\frac{r}{q}\Big)^{m}$ and moreover for any two distinct configurations $\sigma_1,\sigma_2 \in \ff_{q}^{n}$ the random variables $Z_{\sigma_1}$ and $Z_{\sigma_2}$ are independent.
\end{lemma}

\begin{proof}
Notice that $S_{A,R,b}$ is a union of distinct cosets and therefore,
\begin{align*}
S_{A,R,{\bf 0}} \equiv \bigcup_{y \in \{\alpha_0,\alpha_1,\dots,\alpha_{r-1}\}^{m}} S_{A,{\bf 0}}(y),
\end{align*}
where $S_{A,{\bf 0}}(y) \equiv  \{Ax+Ry \mid x \in \ff_q^{n-m}\}$ is defined as a particular coset with a fixed $y \in \{\alpha_0,\alpha_1,\dots,\alpha_{r-1}\}^{m}$. Hence $|S_{A,R,{\bf 0}}|=q^{n-m}r^{m}$ and since $S_{A,R,b}$ is simply a random affine shift of $S_{A,R,{\bf 0}}$, $|S_{A,R,b}|=q^{n-m}r^{m}$ as well. Now for a vector $\sigma \in \ff_q^{n}$, we must have
\begin{align*}
\Pr(Z_{\sigma})=\sum_{y \in S_{A,R,b}} \Pr(\sigma=y)=\frac{|S_{A,R,b}|}{q^{n}}=\Big(\frac{r}{q}\Big)^{m}.
\end{align*}
Next, for two configurations $\sigma_1,\sigma_2 \in \ff_{q}^{n}$, we have that
\begin{align*}
\Pr(Z_{\sigma_1}=1 \wedge Z_{\sigma_2}=1)&=\sum_{y_1,y_2} \Pr(\sigma_1 \in S_{A,b}(y_1) \wedge \sigma_2  \in S_{A,b}(y_2)) \\
&=\sum_{y_1,y_2} \Pr(\sigma_{1} \in S_{A,b}(y_1) \mid \sigma_{2} \in S_{A,b}(y_2))\Pr(\sigma_2 \in S_{A,b}(y_2)) \\
&=\sum_{y_1,y_2} \Pr(\sigma_{1}-\sigma_{2} \in S_{A,{\bf 0}}(y_{1}-y_{2}))\Pr(\sigma_2 \in S_{A,b}(y_2)).
\end{align*}
Therefore we just need to evaluate the probability of the event $\Pr(\tau \in S_{A,{\bf 0}}(z))$ for $\tau=\sigma_1-\sigma_2 \neq 0$ and $z=y_1-y_2$. Now, if $z={\bf 0}$, $\Pr(\tau \in  S_{A,{\bf 0}}(z))$ is simply equal to $\frac{1}{q^{m}}$ since the columns of $A$ are independent, i.e.,  $|\{Ax:x \in \ff_q^{n-m} \}|=q^{n-m}$. Now, since the linearly independent columns of $R$ are sampled randomly.
We have,
\begin{align*}
&\Pr(\tau \in S_{A,{\bf 0}}(z) \mid z \neq 0, \tau  \in \{Ax:x \in \ff_q^{n-m} \})=0 \\
&\Pr(\tau \in S_{A,{\bf 0}}(z) \mid z \neq 0, \tau \notin \{Ax:x \in \ff_q^{n-m} \})=\frac{1}{q^{m}-1}.
\end{align*}
Hence,
\begin{align*}
\Pr(\tau \in \tau \in S_{A,{\bf 0}}(z) \mid z \neq 0)&=
\Big(1-\frac{1}{q^{m}}\Big)\frac{1}{q^{m}-1}=\frac{1}{q^{m}}.
\end{align*}
Therefore, we have that
\begin{align*}
\Pr(Z_{\sigma_1}=1 \wedge Z_{\sigma_2}=1)&=\sum_{y_1,y_2}  \frac{1}{q^{2m}}=\Big( \frac{r}{q} \Big)^{2m}
\end{align*}
and hence we have the statement of the lemma.
\end{proof}

From Algorithm \ref{algo:mbawish} it is clear that Lemma \ref{lem:ach1} allows us to obtain the values of $M_i$ for $i \in \{1,2,\dots,n\}$.
Indeed, the MAX-oracle is not well defined when $m > n$.
 In order to obtain the values of $M_{i}$ for $i \in \{n+1,\dots,n'\}$, we propose the following technique.

Recall that the elements of $\ff_{q}^{n}$ can be represented as $n$ dimensional vectors where each element belongs to $\ff_{q}$. Moreover we defined an ordering over the elements of the finite field $\ff_{q}\equiv  \{\alpha_0, \alpha_1,\dots , \alpha_{q-1}\}$ so that $\alpha_{i} < \alpha_{j}$ for $i <j$.
Consider the lexicographic ordering of the elements (vectors) of $\ff_{q}^{n}$.

Let $s_m$ be the $\lceil \frac{r^{m}}{q^{m-n}} \rceil$ th element in this ordering of $\ff_{q}^{n}$. Define the set
\begin{align}\label{def:si}
\mathcal{S}_m=\{x \in \ff_{q}^{n} \mid x < s_m \}
\end{align}
for all $m>n$. Generate a random uniform full rank matrix $A \in \ff_{q}^{n \times n}$ and a uniform random vector $b \in \ff_{q}^{n}$. Subsequently, the MAX-Oracle for \texttt{MBA-WISH} solves the following optimization problem for $m>n$:
\begin{align*}
\max_{y\in \mathcal{S}_m} w(Ay+b).
\end{align*}
In order to analyze the statistical properties of this oracle, define the random set
\begin{align*}
T_{A,b,m} \equiv \{Ay+b \mid y \in \mathcal{S}_{m}\}.
\end{align*}
Again, for each configuration $\sigma \in \ff_{q}^{n}$, associate an indicator random variable $Z_{\sigma}$ denoting $\sigma \in T_{A,b,m}$.

\begin{lemma}
For each configuration $\sigma \in \ff_{q}^{n}$, we must have $\Pr(Z_{\sigma}=1)\approx  \Big(\frac{r}{q}\Big)^{m}$ and moreover for any two configurations $\sigma_1,\sigma_2 \in \ff_{q}^{n}$,  $\Pr(Z_{\sigma_1}=1 \wedge Z_{\sigma_2}=1) \le (\Pr(Z_{\sigma}=1))^2.$
\end{lemma}
\begin{proof}
We have,
\begin{align*}
\Pr(Z_{\sigma}=1)=\sum_{y \in T_{A,b,m}} \Pr(\sigma=y)=\frac{|T_{A,b,m}|}{q^{n}}=\frac{1}{q^n}\Big\lfloor \frac{r^{m}}{q^{m-n}} \Big\rfloor \approx \Big(\frac{r}{q}\Big)^{m}.
\end{align*}
Next, for two distinct configurations $\sigma_1,\sigma_2 \in \ff_{q}^{n}$, we have that
\begin{align*}
\Pr(Z_{\sigma_1}=1 \wedge Z_{\sigma_2}=1)&=\sum_{y_1,y_2 \in T_{A,b,m}} \Pr(\sigma_1=y_1 \wedge \sigma_2=y_2) \\
&=\sum_{y_1,y_2 \in \mathcal{S}_{m}} \Pr(\sigma_1=Ay_1+b \wedge \sigma_2=Ay_2+b) \\
&=\sum_{y_2 \in \mathcal{S}_{m}} \Pr(\sigma_2=Ay_2+b)\sum_{y_1 \in \mathcal{S}_m} \Pr(\sigma_1=Ay_1+b \mid \sigma_2=Ay_2+b) \\
&=\sum_{y_2 \in \mathcal{S}_{m}} \Pr(\sigma_2=Ay_2+b)\sum_{y_1 \in \mathcal{S}_m} \Pr(\sigma_1-\sigma_2=A(y_1-y_2)) \\
\end{align*}
Since $\sigma_1 \neq \sigma_2$, we must have that $\Pr(\sigma_1-\sigma_2=A(y_1-y_2) \mid y_1=y_2)=0$. For $y_1 \neq y_2$, every configuration $\sigma \in \ff_{q}^{n}, \sigma \neq 0$ is equally probable to be $A(y_1-y_2)$ since $A$ is uniformly and randomly sampled full rank matrix. Hence,
\begin{align*}
\Pr(Z_{\sigma_1}=1 \wedge Z_{\sigma_2}=1)&=\frac{1}{q^{n}(q^n-1)}\Big\lfloor \frac{r^{m}}{q^{m-n}} \Big\rfloor \Big(\Big\lfloor \frac{r^{m}}{q^{m-n}} \Big\rfloor-1\Big) \le (\Pr(Z_{\sigma}=1))^2.
\end{align*}

\end{proof}

The remainder of the proof of Theorem~\ref{thm:mba} follows that of Theorem \ref{thm:mb} in a straightforward manner.

\section{The \texttt{PERM-WISH} Algorithm}\label{sec:perm_}

\begin{algorithm}[H]
\caption{\texttt{PERM-WISH} for $\Omega=S_n$}                               \label{algo:permanent}
\begin{algorithmic}
      \REQUIRE $\ell \rightarrow \lceil \frac{1}{\gamma}\ln \frac{2n}{\delta} \rceil$, $q>n$, $r,n'=\lceil n \log_{q/r} q\rceil$
      \STATE $M_0 \equiv \max_{\sigma \in S_n} w(\sigma)$
       \FOR{$i\in\{1,2,\dots,n'\}$}
         \FOR{$k \in \{1,\dots,\ell\}$}
           \STATE Sample hash functions $h_i \equiv h_{A^{i},b^{i}}$ uniformly at random from $\Hc_{i,n}$ as defined in \eqref{eq:fam}
           \STATE $w_{i}^{(k)}=\max_{\sigma\in \ff_q^n} w(\sigma) $ such that $ A^{i}\sigma +b^i < \alpha_r\cdot \mathbf{1} ; \sigma < \alpha_{n-1}\cdot \mathbf{1} ;\sigma(k) \ne \sigma(l) \forall k \ne l$.
        \ENDFOR
         \STATE $M_{i}={\rm Median}(w_{i}^{(1)},w_{i}^{(2)},\dots,w_{i}^{(\ell)})$
       \ENDFOR
       \STATE Return $M_{0}+(\frac{q}{r}-1)\sum_{i=0}^{n'-1}M_{i+1}\big(\frac{q}{r}\big)^{i} $
\end{algorithmic}
\end{algorithm}

\section{Omitted Proofs from Section~\ref{sec:der}} \label{app:der}

\begin{proof}[Proof of Prop.~\ref{prop:aff}]
It is sufficient to show that $z$ can be obtained as a linear transform of $\nu$. Note that the product of $\phi(x)$ and $\phi(\nu)$ can be written as a convolution between $x$ and $\nu\equiv (\nu_1, \nu_2, \dots, \nu_n)$ (as we can view this as product between two polynomials). Let $\Gamma$ be the $(2n -1) \times n$ matrix,
$$ \Gamma=
\begin{bmatrix}
\nu_1 & 0 & 0 &  \dots & 0\\
\nu_2 & \nu_1 & 0 & \dots & 0\\
\nu_3& \nu_2 & \nu_1 & \dots & 0\\

\vdots & \vdots & \vdots &\vdots & \vdots\\
\nu_n & \nu_{n-1} & \nu_{n-2} & \dots & \nu_1\\
0 & \nu_n & \nu_{n-1} & \dots & \nu_2\\
\vdots & \vdots & \vdots &\vdots & \vdots\\
0 & 0 & 0 & \dots & \nu_n
\end{bmatrix}.
$$
The reduction modulo $f(\zeta)$ can also be written as a linear operation. Just consider the $n \times (2n-1)$ matrix
$P$ whose $i$th column contains the coefficients of the polynomial $\zeta^{i-1} \mod f(\zeta), 1\le i \le 2n-1.$ Note that the first $n$ columns of the matrix is simply the identity matrix.
We can write,
$
z= P\Gamma x.
$
\end{proof}

\vspace{0.2in}

\begin{proof}[Proof of Prop.~\ref{prop:pair}]
Suppose $\nu,b$ are randomly and uniformly chosen. For any $x_1,x_2\in \ff_2^n$ and $y_1,y_2 \in \ff_2^m$, first of all
$$
\Pr(h_{\nu,b}(x_1) = y_1) = \frac{1}{2^m},
$$
since $b$ is uniform. Now,
\begin{align*}
&\Pr(h_{\nu,b}(x_1) = y_1, h_{\nu,b}(x_2) = y_2) \\
&= \frac{1}{2^m}\Pr(h_{\nu,b}(x_2) = y_2 | h_{\nu,b}(x_1) = y_1)\\
&= \frac{1}{2^m}\Pr(h_{\nu,b}(x_2) -h_{\nu,b}(x_1)= y_2-y_1 | h_{\nu,b}(x_1) = y_1)\\
&= \frac{1}{2^m}\Pr(h_{\nu,b}(x_2) -h_{\nu,b}(x_1)= y_2-y_1).
\end{align*}
Now, since $\Pr((\phi(x_1)-\phi(x_2)) \cdot \phi(\nu) \mod f(\zeta) = u) = \frac1{2^n}$ for any $u$, we must have $\Pr(h_{\nu,b}(x_2) -h_{\nu,b}(x_1)= y_2-y_1)= \frac1{2^m}$.
Therefore the claim is proved.

\end{proof}

\vspace{0.2in}

\begin{proof}[Proof of Prop.~\ref{prop:top}]

First of all, the uniformity of the family is immediate since $b$ is uniformly chosen.
For any $x_1,x_2\in \ff_2^n$ and $y_1,y_2 \in \ff_2^m$,
$
\Pr(h_{A_T,b}(x_1) = y_1, h_{A_T,b}(x_2) = y_2)
= \frac{1}{2^m}\Pr(h_{A_T,b}(x_2) = y_2 | h_{A_T,b}(x_1) = y_1)

= \frac{1}{2^m}\Pr(A_T(x_1-x_2)= y_1-y_2).
$
It remains to prove that $\Pr(A_Tx =y) = \frac{1}{2^m}$ for any fixed $x,y$.  Let the $k$th coordinate of $x$ is the first to be in the support of $x$.
Now consider the inner product of the $j$th row of $A_T$ with $x$. This product will contain the entry $A_T(j,k)$, the  $(j,k)$th entry of $A_T$. Note that, this entry would not have appeared in any of the inner products of $i$th row of $A_T$ and $x$, for $i <j$. Therefore the probability that this inner product is any fixed value is exactly $\frac12$ given inner product of all previous rows with $x$. Therefore, $\Pr(A_Tx =y) = \frac{1}{2^m}$.

\end{proof}

\end{document}


TITLE: expected value of an exponential distribution
QUESTION [0 upvotes]: I've been stuck at this problem for hours now. So first I know that the expected value of $X\sim\operatorname{Exponential}(\lambda)$ is $1/\lambda.$ But I can't figure out what the expected value is given this information:
$X\sim\operatorname{Exponential}(3)$ and $Y=\exp(2x)$
So, how do I calculate the $\operatorname E[Y]$? A step by step would be helpful, I've been stuck for hours.

REPLY [1 votes]: First, if $X \sim Exp(\lambda)$, then it has a p.d.f. given by $f_\lambda(x) = \lambda \mathrm{e}^{-\lambda x}$. By the definition of expectation:
$$ \mathrm{E}[Y] = \mathrm{E}[\mathrm{e}^{2X}] = \int_0^\infty \mathrm e^{2x} f_3(x) \mathrm{d} x = 3 \int_0 ^\infty \mathrm{e}^{2x}\mathrm{e}^{ - 3 x}\mathrm d x = 3 \int_0^\infty \mathrm{e}^{-x}\mathrm{d}x = 3.$$

REPLY [1 votes]: Let $X\sim\mathsf{Expo}(\lambda)$ and $Y=e^{aX}$ with $\lambda>a>0$. Then for $y>1$ we have
\begin{align}
\mathbb P(Y\leqslant y) &= \mathbb P\left(e^{aX}\leqslant y\right)\\
&= \mathbb P\left(X\leqslant \frac{\log y}a\right)\\
&= 1 - e^{-\frac{\lambda}a\log y}\\
&= 1 - \left(e^{\log y}\right)^{-\frac\lambda a}\\
&= 1 - y^{-\frac\lambda a}.
\end{align}
It follows that
\begin{align}
\mathbb E[Y] &= \int_1^\infty \left(1-\mathbb P(Y\leqslant y)\right)\ \mathsf dy\\
&= \int_1^\infty y^{-\frac\lambda a}\ \mathsf dy\\
&= \frac\lambda{\lambda-a}.
\end{align}
If $\lambda=3$ and $a=2$ then this is equal to $3$.


TITLE: Weighted P-norm
QUESTION [2 upvotes]: It is known that $ l^p \subseteq l^r $ if $ 1 \le p \le r < \infty $. Also, $ l_w^p $ are spaces where we consider a collection of positive weights, $ w_n > 0 $ for each $ n \in \mathbb{N} $, that given sequence $ x = (x_n) $, the weighted $p$-norm is
 $$ ||x||_{w,p} = \big (\sum_{n=1}^\infty w_n |x_n|^p \big )^\frac{1}{p} \text{.} $$
 Prove that if $ \sum_{n=1}^\infty w_n < \infty $, then $ l^1_w \not \subseteq l^2_w $.One way to prove this is to find a subsequence $ (w_{n_k}), $ so that $ w_{n_k} \le \frac{1}{2^k} $ and consider elements which are nonzero only at the $ n_k $ coordinates.
the subsequence is probably very easy to think up, but I can just not figure it out... :(

REPLY [2 votes]: Let $(n_k)$ be the increasing sequence of integers mentioned in the opening post and define
$x_{n_k}=\left(w_{n_k}k^2\right)^{-1}$ and $x_l=0$ if $l$ is not equal to $n_k$ for some $k$. The series $\sum_{k}w_{k}\left\lvert x_{k}\right\rvert$ is convergent and $\sum_{k}w_{n_k}x_{n_k}^2$ is divergent since $w_{n_k}x_{n_k}^2\geqslant 2^k/k^4$.


TITLE: Embedding of $C(X)$ into $B(H)$ where $H$ is separable
QUESTION [6 upvotes]: I would like to ask a question which may look strange at the first sight nevertheless I find it interesting. Let $H$ be a separable Hilbert space: for any separable $C^*$-algebra $A$ one can embed $A$ into $B(H)$. Very often one can do even better and many nonseparable $C^*$-algebras still embed into $B(H)$, for example $\ell^{\infty}, L^{\infty}[0,1]$ or $B(H)$ itself. However it is not always possible even for quotients of such algebras, $\ell^{\infty}/c_0$ being an example. One can view $\ell^{\infty}$ as $C(\beta\mathbb{N})$ where $\beta\mathbb{N}$ is the Stone-Cech compactification of $\mathbb{N}$ (which is of cardinality $2^{\mathfrak{c}}$) and $\ell^{\infty}/c_0$ as $C(\beta\mathbb{N} \setminus \mathbb{N})$ (still of cardinality $2^{\mathfrak{c}}$). So here comes my question:

Let $X$ be a compact nonmetrizable Hausdorff topological space of cardinality $\mathfrak{c}$. Is it true that $C(X)$ embeds into $B(H)$ for separable $H$?

REPLY [7 votes]: No, not necessarily.  If $C(X)$ embeds into $B(H)$ then $X$ must be ccc, but there exist compact Hausdorff spaces of cardinality $\mathfrak{c}$ which are not ccc.  (One example is to take $\mathfrak{c}$ with its discrete topology and form the one-point compactification $\mathfrak{c} \cup \{\infty\}$; then all the singletons except $\infty$ are disjoint open sets.  Another is the ordinal $\mathfrak{c}+1$ with its order topology; all the successor singletons $\{\alpha+1\}$, $\alpha < \mathfrak{c}$, are disjoint open sets.)
To see this, suppose that $X$ is not ccc, so that there is an uncountable family $\{U_\alpha\}$ of disjoint nonempty open sets.  By the appropriate version of Urysohn's lemma, for each $\alpha$ there is a nonzero real-valued $f_\alpha \in C(X)$ supported inside $U_\alpha$; in particular, $f_\alpha f_\beta = 0$ for $\alpha \ne \beta$.  Now if there is an injective *-homomorphism $\Phi : C(X) \to B(H)$, then each $\Phi(f_\alpha)$ is nonzero, so we may find $h_\alpha \in H$ so that $\|\Phi(f_\alpha) h_\alpha\|=1$.  Also, $\Phi(f_\alpha)$ is self-adjoint, so for $\alpha \ne \beta$ we have
$$\langle \Phi(f_\alpha) h_\alpha, \Phi(f_\beta) h_\beta\rangle = \langle  h_\alpha, \Phi(f_\alpha)\Phi(f_\beta) h_\beta\rangle = \langle  h_\alpha, \Phi(f_\alpha f_\beta) h_\beta\rangle = 0$$
so that $\{\Phi(f_\alpha) h_\alpha\} \subset H$ is an uncountable orthonormal set, which is impossible if $H$ is separable.

This is fairly close to optimal because every separable compact Hausdorff $X$ does have $C(X)$ embeddable into $B(H)$ (recall that separable implies ccc but not conversely).  Indeed, if $\{x_n\}$ is a countable dense subset of $X$, then identifying $f \in C(X)$ with the sequence $(f(x_1), f(x_2), \dots)$ gives an isometric embedding of $C(X)$ into $\ell^\infty$, which as we know embeds in $B(\ell^2)$.  Presumably there is a necessary and sufficient condition somewhere in between, but I don't know what it is.


TITLE: Bijective proof for equal number of odd-sized/even-sized subsets of a finite set
QUESTION [2 upvotes]: Possible Duplicate:
Number of even and odd subsets

Let $X$ be a finite set of cardinality $n$. Let ${\cal E}$ ($\cal O$ )denote the set of all subsets of $X$ with an even (odd) number of elements. It is easy to see that $\cal E$ and $\cal O$ contain the same number of subsets ; indeed, Newton’s binomial formula gives
$$
0=0^n=(1-1)^n=\sum_{k=0}^n \binom{n}{k} (-1)^k=
\sum_{k \ \text{even}} \binom{n}{k} -\sum_{k \ \text{odd}} \binom{n}{k}=
|{\cal E}|-|{\cal O}|
$$
When $n$ is odd, there is a simpler, more combinatorial proof  : $Y \mapsto X \setminus Y$ defines a bijection between $\cal E$ and $\cal O$. Is there a similar, combinatorial, bijection between  $\cal E$ and $\cal O$ when $n$ is even ?

REPLY [7 votes]: First note that the claim only holds if $n>0$. There is a general combinatorial proof that works for $n$ even or odd: since $n>0$ there is at least one element in  $X$, call it $x_0$. Now define $f:E\to O$ by mapping $A\in E$ to $A\cup \{x_0\}$ if $x_0\notin A$ and to $A-\{x_0\}$ if $x_0\in A$. $f$ is easily seen to have an inverse (defined in a completely analogous way to $f$'s definition), and so is a bijection.


TITLE: Solving a system of equations using modular arithmetic
QUESTION [4 upvotes]: I am trying to implement a solver for the game lights out. You have a grid of lights, when you click on one of them the light you clicked and its four neighbours change colour, with the light starting over when it runs out of colours. The aim is to get all the lights to a particular colour.
Because the way the colours change is cyclic, I thought I could implement it as a system of equations and do all calculations mod n (n being the number of available colours).
This method worked for some puzzles but I got stuck in others.
I am representing the game as a system of equations in an augmented matrix and reducing it to reduced row echelon form using gaussian elimination. As I said in some cases this worked well. However there are some cases (example to follow) where I end up with a line which I cannot reduce completely, reason being that, since I'm using modular arithmetic, some values don't have a multiplicative inverse so I get stuck.
Here is an example:
The game shown here represents a 4x4 grid with 4 colours available. Here is the matrix as it started out:
1100100000000000 3
1110010000000000 3
0111001000000000 2
0011000100000000 3
1000110010000000 3
0100111001000000 3
0010011100100000 0
0001001100010000 0
0000100011001000 0
0000010011100100 0
0000001001110010 2
0000000100110001 1
0000000010001100 3
0000000001001110 1
0000000000100111 0
0000000000010011 0

and this is as far as I've managed to reduce it:
1000000000000333 2
0100000000003323 3
0010000000003233 0
0001000000003330 0
0000100000001232 2
0000010000002003 2
0000001000003002 3
0000000100002321 3
0000000010001320 1
0000000001003332 1
0000000000102333 0
0000000000010231 2
0000000000000220 2
0000000000002222 0
0000000000002222 0
0000000000000220 2

In this case the last line, for example, is ...220 2 and I cannot figure out how I can reduce it since I cannot simply divide by 2 (2 has no multiplicative inverse in z4). Whatever I tried, I always ended up with a leading 2 in a row. I am absolutely certain a solution exists for the game (I did solve it correctly myself) but I'm not sure if I'm doing missing something here, or just that the method simply does not work for all cases.
Thanks
Edit: Fixed the matrices in the example. There were inconsistencies at the end. I've figured out that any inconsistencies there mean there is no solution. In the updated example a solution does exist for sure, and this results in no inconsistencies at the end. However I still cannot solve it

REPLY [0 votes]: Here is a probably much more useful answer than my previous one. When solving a system of equations, one can use row operations freely (as usual modifying the right hand side correspondingly). One can however also allow column operations provided one keeps track of them. For instance, one could add twice the first column to the third column, which amounts to replacing the first unknown $x_1$ by a new unknown $x_1'=x_1-2x_3$, keeping all other unknowns $x_i$ (one has $ax_1+cx_3=ax'_1+(2a+c)x_3$ for any $a,c$, which explains the column operation). After solving the new system, one must not forget to recover $x_1=x'_1+2x_3$ from the solution in terms of the new variables. In matrix terms, if we modify the system $A\cdot x=b$ by column operations that amounts to right-multiplication by an invertible matrix $R$, then we get a new system $AR\cdot x'=b$ in terms of the variable $x'=R^{-1}\cdot x$, and at the end one then computes $x=R\cdot x'$. So each column operation applied to$~A$ is accompanied by the same column operation applied to$~R$ (which starts out as the identity), and at the end $R$ is applied (from the left) to the solution vector $x'$ found to obtain the true solution $x$ of the original system.
Given that one can mix row and column operations (recording the latter), one can now arrange that any pivot ultimately used divides all remaining entries (those to it bottom-right), and also the modulus$~n$. The way to achieve this is inspired by Euclid's algorithm, and very similar to the algorithm for reducing a matrix to its Smith normal form. Basically one continually reduces the size of the minimal remaining nonzero entry$~d$ by computing the remainder of the Euclidean division of another remaining entry by$~d$ as long as this is possible (i.e., while $d$ does not divide all others): if that other entry is in the same row or column as$~d$, a single row or column operation will do; if not, one performs a row operation to make a non-multiple of$~d$ appear in the column of$~d$, of which one the forms the remainder by a column operation. As a difference with the Smith normal form algorithm over$~\Bbb Z$, the entries are now all defined modulo$~n$, and so all arithmetic is done modulo$~n$. This in fact allows a significant extra operation to reduce$~d$: whenever $d$ does not divide$~n$, just form a multiple of its row so that$~d$ becomes smaller after reduction modulo$~n$; one easily sees that it can be reduced to $\gcd(d,n)$ in this manner. All in all, even though there are a large number of possibilities to be considered in general (the program is more complicated than for Gaussian elimination), in practice the reduction of$~d$ is very rapid, and most of the time (until the very end of the procedure) one will succeed to reduce it to$~1$.
With this normal form algorithm the system $A\cdot x=b\pmod n$ can be transformed into such a system in which the only nonzero entries of$~A'$ are on the main diagonal, and each one divides its successor. Once in this form, the solution of the system is obvious; to have any solution, each diagonal entry of$~A'$ should divide the corresponding entry of$~b$ (and in particular every zero row of$~A'$ should have the corresponding entry of$~b$ zero.


TITLE: Elementary set theory - Confused about A $\cup$ B = B
QUESTION [3 upvotes]: I am trying to solve:
Proof that  $\exists!A\subset U(\forall B \subset U(A \cup B = B)$.
But once I think about it, I can see not one but 2 solutions.
a) $ \varnothing \cup B = B $
b) $ B \cup B = B $
The solution is supposed to be a), what is wrong about b)?
Thanks!

REPLY [2 votes]: It needs to work for all possible sets $B$.  Not just one specific $B$.
For example: if $U =\mathbb R$ the we need to find an $A$ so that $A \cup B = B$ for ALL $B$.  So for example $A \cup [5, 7.3) = [5, 7.3)$ and $A \cup \mathbb Q = \mathbb Q$ and $A \cup \{-27, 0, \sqrt{\pi}, e^{-\frac 23}\}= \{-27, 0, \sqrt{\pi}, e^{-\frac 23}\}$.
It should be clear from those examples that the only possible set that will do that is $A = \emptyset$.
We can formally prove that.  But notice because $A\cup [5,7.30)  = [5,7.30)$ then $2 \not \in A\cup [5,7.3)$ so $2 \not \in A$ and $A \cup \{-27, 0, \sqrt{\pi}, e^{-\frac 23}\}= \{-27, 0, \sqrt{\pi}, e^{-\frac 23}\}$ so $5\frac 12 \not \in A \cup \{-27, 0, \sqrt{\pi}, e^{-\frac 23}\}= \{-27, 0, \sqrt{\pi}, e^{-\frac 23}\}$ so $5\frac 12 \not \in A$, etc.  Can anything be in $A$?  Well, no because.... if $x \in U$ then $A\cup \{x\}^c = \{x\}^c$ and $x \not \in A\cup \{x\}^c$ so $x \not \in A$.
=====
By the way.  Using your incorrect  interpretation of the question: If $B\subset U$ then show there exists a $A$ so that $A \cup B = B$. (a different question than what was asked.) Then the answer to that can be any subset of $B$.  Not just $\emptyset$ or $B$ but for any $A \subset B$ then $A \cup B = B$.
But that's an entirely different question.


TITLE: do numerical integration with fixed abscissas
QUESTION [1 upvotes]: Can I do an integral, possibly using gaussian quadrature, when the abscissas are fixed (for reasons that I don't want to get into right now), i.e. is it possible to calculate the weights for fixed abscissas that I don't get to choose?

REPLY [1 votes]: I have settled this question here; in brief, a FORTRAN implementation of algorithms for solving this problem have been published in the Collected Algorithms of the ACM.

At the risk of being redundant, I'd like to mention here some other things I mentioned in the (now expanded) answer of mine at that other site for completeness' sake.
Robert's warning of the Runge phenomenon happening is a good one, and it does happen if your abscissas are perversely distributed (relatedly, the underlying Vandermonde matrix is ill-conditioned); the equispaced case being among the worst-behaved point distributions. Abscissas that are "nicely distributed" (e.g. Legendre, Chebyshev, or any other point distribution which "clusters" near the ends) will generally ensure that you have a quadrature rule that behaves tamely even for large numbers of points. (As an aside, the chebfun project hinges on functions being nicely approximated by interpolating polynomials with abscissas that are transformed Chebyshev polynomial roots/extrema.)
Lastly, whatever you finally settle with, you will want to perform the sanity check of ensuring that all the weights of your quadrature rule are of the same sign; any change of sign in the weights can lead to subtractive cancellation when you employ the quadrature rule, and you wouldn't want that...


TITLE: Why do siamese magic squares have real eigenvalues, symmetric around zero?
QUESTION [20 upvotes]: There is a standard method to construct magic squares of odd size, known as the
 Siamese construction. I'll write $S_m$ for the $m \times m$ Siamese square. For example, here is $S_5$.
17    24     1     8    15
23     5     7    14    16
 4     6    13    20    22
10    12    19    21     3
11    18    25     2     9

As you can see, this matrix certainly isn't symmetric! It has one obvious eigenvector, the all ones vector, whose corresponding eigenvalue is $\frac{m^3+m}{2}$, and there is no obvious reason its other eigenvalues should have nice structure. Nonetheless, experimentation suggests:

(1) All eigenvalues of $S_m$ are real.
(2) If $\lambda$ is an eigenvalue of $S_m$ other than $\frac{m^3+m}{2}$, then $-\lambda$ is as well.

Can anyone explain why?
Remark In case anyone is curious how I found this, I am teaching myself MATLAB and came up with the task of listing characteristic polynomials of magic squares as a more or less random task which would involve learning some basic control structures and some mathematical tools. This piece of MATLAB documentation describes some alternate ways of thinking about the Siamese construction, which may be relevant.

REPLY [4 votes]: Here are some partial results.
This paper states that every regular magic square $\mathbf A\in\mathbb N^{n\times n}$ (that is, its entries satisfy $a_{i,j}+a_{n-i+1,n-j+1}=n^2+1$) is representable as a rank-$1$ correction of a skew-centrosymmetric matrix $\mathbf Z$ (that is, $\mathbf Z=-\mathbf J\mathbf Z\mathbf J$, where $\mathbf J$ is the exchange matrix, or the identity matrix with its columns reversed). In other words, the matrix you get when you subtract $\frac{n^2+1}{2}$ from all of the entries of $\mathbf A$ is skew-centrosymmetric. This is important since the eigenvalues of $\mathbf A$, apart from the magic eigenvalue $\frac{n^3+n}{2}$, are the same as the eigenvalues of $\mathbf Z$, which now has $0$ as an eigenvalue.
We can now use the results of this paper, which says that if $\lambda$ is an eigenvalue of $\mathbf Z$, then $-\bar{\lambda}$ must also be an eigenvalue, with the same algebraic multiplicity.
Thus, the remaining task is to prove that all the eigenvalues of a Siamese magic square are necessarily real.

A tantalizing result: Pressman, in this paper, shows the stronger result that if $\lambda$ is an eigenvalue of $\mathbf Z$, then $-\lambda$ is also an eigenvalue. There seems to be a criterion in the paper for saying when a skew-centrosymmetric matrix has all its eigenvalues real, but I have not yet figured how to apply it to this case.


TITLE: $n$ by $n$ Primally Magic Squares
QUESTION [5 upvotes]: (Again copied verbatim from a September 2009 thread I made.)
A Primally Magic Square (PMS) is exactly like a traditional magic square with a change of criteria. Where a traditional magic square is one where the sum of each row, column, and both diagonals are some number, usually 15 (uses 1-9). By contrast, a PMS is one where all rows, columns, and both diagonals form a prime number. Note: if a number is prime either forward or back, it counts.
Example:
$\displaystyle \begin{matrix} 1 & 3 & 7 \\ 3 & 9 & 7 \\ 7 & 7 & 9 \end{matrix}$
137, 397, 977, 797, and 199 are all prime, so this square is a PMS.
How many 3x3 PMS's are there?
How many 4x4 PMS's are there?
How many 5x5 PMS's are there?
How many 9x9 PMS's are there?
Side note: Feel free to calculate the number of PMS's for some other $n$.

REPLY [0 votes]: Following a comment, let's look at the $2×2$ case.  I shall consider the case where all digits are distinct and to avoid issues with convention, nonzero.
Within these limits we select a set of four digits such as $\{1,2,3,4\}$.  If each of the six two-element permutations generates a prime in either order, we have three primary magic squares apart from rotations and reflections (depending on which element is opposite a given one).
The set $\{1,2,3,4\}$ fails because $42$ and $24$ are both composite, and similarly there cannot be any two even digits, nor the digit $5$ (only two numbers in the 50s are prime and we would need three).  We also cannot have two digits with  residues $1,2\bmod 3$ or two multiples of $3$ (which lead to multiples of $3$), nor specifically $4$ together with $9$ (the odd number $49$ is $7×7$).
Here are the sets that work:
$\{1,3,4,7\}$
And that seems to be it!  Thus three squares (plus trivial operations):
$\begin{array}&1&3\\4&7\end{array}$
$\begin{array}&1&4\\7&3\end{array}$
$\begin{array}&1&7\\3&4\end{array}$


TITLE: $I_n=\int_{\pi}^{2\pi} \frac {\left|\sin(nx)\right|}{x}dx$ 12th grade, inequality
QUESTION [0 upvotes]: Let $I_n=\int_{\pi}^{2\pi} \frac {\left|\sin(nx)\right|}{x}dx$
($a$) Prove $I_n \leq \ln(2)$
($b$) Prove $I_n \geq \frac {2}{\pi}(\frac 1{n+1}+ \frac 1{n+2}+...+\frac 1{2n})$

For ($a$)
We have
$$\left|\sin(nx)\right|\leq 1 \forall x\in\mathbb{R}$$
Let's divide by $x$ and integrate from $\pi$ to $2\pi$
We get:
$$\int_{\pi}^{2\pi}\frac{\left|\sin(nx)\right|}{x}dx\leq \int_{\pi}^{2\pi}\frac 1xdx$$
Leading to
$$I_n \leq \ln(2)$$
For ($b$) however I don't know how to start. Actually I tried to apply substitution $u=nx$ and then do a sum of integrals but doesn't really work.

REPLY [3 votes]: Note that by letting $t=nx$ and then $s=t-(n\pi+ k\pi)$ we have that
$\left|\sin(t)\right|=\left|\sin(s)\right|$, and
$$\begin{align}I_n&=\int_{\pi}^{2\pi} \frac {\left|\sin(nx)\right|}{x}dx=
\int_{n\pi}^{2n\pi} \frac {\left|\sin(t)\right|}{t}dt\\
&
=\sum_{k=0}^{n-1}\int_{t=n\pi+k\pi}^{n\pi+(k+1)\pi}\frac {\left|\sin(t)\right|}{t}dt=
\sum_{k=0}^{n-1}\int_{s=0}^{\pi}\frac {\left|\sin(s)\right|}{n\pi+k\pi+s}ds\\&\geq \sum_{k=0}^{n-1}\frac {1}{n\pi+k\pi+\pi}\int_{s=0}^{\pi}\left|\sin(s)\right|ds=\frac{2}{\pi}\left(\frac 1{n+1}+ \frac 1{n+2}+\dots+\frac 1{2n}\right).\end{align}$$


TITLE: Do reversible processes exist in Nature?
QUESTION [2 upvotes]: Classical thermodynamics always discusses entropy in the light of reversible processes, and it lies at the heart of the definition of entropy. But do these reversible processes exist in Nature, or are they just a gedankenexperiment?
Since if they do not exist, then classical thermodynamics gives a lower bound on the increase of entropy, but says nothing about the difference in entropy between the reversible and irreversible process. Doesn't this make the classical definition useless for a quantitative approach? Or is there proof that the entropy difference between the reversible and irreversible process is small for most everyday cases?

REPLY [3 votes]: Of course, thermodynamics says things about differences in entropy between reversible and irreversible processes. In fact, we can analyse how irreversible a process was and what to do in order to, for instance, extract more work from it. That's a very important part of thermodynamics.
But, no, there are no truly reversible processes — at least for daily macroscopic, bulk phenomena; I might be unaware of some crazy quantum-relativistic-boson-fermion experiment :). In practice, all real processes are irreversible, even though some can be good approximations of reversible ones (e.g. if you go really slowly).
Also, your statement "classical thermodynamics gives a lower bound on the increase of entropy" is a bit weird. That "lower bound" would be zero, so I'd prefer to say that it states that global entropy never decreases. Apart from that, classical thermodynamics also goes on to explore quantities like availability, irreversibility and efficiency, and it can do a good work of optimizing the hell out of your irreversible processes.


TITLE: Flip a coin 10 times, probability that no face shows up three times in a row?
QUESTION [1 upvotes]: Given a fair coin, if we flip the coin $10$ times, what is the probability that the string of flips does not contain a $'HHH'$ or $'TTT'$ anywhere inside it?
That is, what is the probability that neither heads or tails shows up three times in a row? (As the title suggests)

Total number of outcomes: $2^{10}$
From seeing related posts and doing some googling, it seems a recurrence relation can be used to solve a problem of this sort, but I got stuck. I don't know if there is another way to solve this kind of problem. (Besides brute force of counting the outcomes we desire) The answer given to this question by N.F Taussig gave me a good start to setting up the problem, but their answer deals with the stipulation that only heads doesn't show up twice in a row.
Letting $a_n$ be the number of sequences of length $n$ that do not have $3$ consecutive heads or tails in them.
It is relatively straightforward to see that $a_1 = 2$, $a_2 = 4$, and $a_3 = 6$
If our sequence of length $n$ ends with $HH$, it must have a $T$ in the $n-2$ position.
If our sequence ends with $TT$, it must have a $H$ in the $n-2$ position.
I am unsure of how to go on from here.

I would also be interested in generalizing this process to $N$ flips of a fair coin where neither heads or tails shows up $k$ times in a row. ($k \ge 2$)
Thanks for any help given.

REPLY [2 votes]: One approach is coupled recurrence relations.  Let $b_n$ be the number of $n$ length sequences that end in one of the same letter, so they end in HT or TH (or are length $1$), and $c_n$ be the number of $n$ length sequences that end in two of the same, so HH or TT.  We then have $$a_n=b_n+c_n\\b_n=b_{n-1}+c_{n-1}\\c_n=b_{n-1}\\b_1=2\\b_2=2\\c_1=0\\c_2=2$$
where the $b_n$ comes because we can we can put the other flip after any series and the $c_n$ comes because we can put another on if there is only one at the end so far.  This gives $$b_n=b_{n-1}+b_{n-2}$$ so $b_n$ is just twice the Fibonacci numbers and $c_n$ is twice the Fibonacci numbers offset by $1$.
$$b_n=2F_n\\c_n=2F_{n-1}\\a_n=2F_n+2F_{n-1}=2F_{n+1}$$