Physics-inspired derivations of some algorithms for computing the permanent

TL;DR

This paper introduces physics-inspired methods for computing the matrix permanent, using Grassmann integrals and fermion problems, leading to new derivations and estimators that connect to existing algorithms.

Contribution

It formulates the permanent computation as a Grassmann integral, deriving known and new algorithms through physics-inspired techniques.

Findings

Derivation of the Godsil-Gutman and Karmarkar estimators from Grassmann integrals.

Connection of the permanent computation to interacting many-fermion problems.

Introduction of new estimators and formulas based on gauge invariance.

Abstract

We provide physics-inspired derivations of a number of algorithms for computing the permanent of a matrix. In particular we formulate the computation of the permanent as a Grassmann integral that may be viewed as an interacting many-fermion problem. Applying a discrete Hubbard-Stratonovich decoupling then gives approximation schemes that are equivalent to the familiar determinant Monte Carlo algorithm. This leads to elementary derivations of the well-known estimators of Godsil-Gutman and Karmarkar et al. Another straightfoward manipulation of the Grassmann integral, making use of gauge invariance, gives the efficient exact formula of Glynn. In addition to these known results we also give some additional estimators and formulas that are natural in our formulation.

Equations

$$\mathrm{per}A=\sum _{\mathcal{P}}\prod _{i=1}^n{A}_{i{\mathcal{P}}_i},$$

$$\int\bigl{[}d\xi^{*}d\xi\bigr{]}e^{-\sum_{i,j}\xi^{*}_{i}A_{ij}\xi_{j}}=\det A.$$

$$\operatorname{\text{per}}A=\int\bigl{[}d\phi^{*}d\phi\bigr{]}e^{\sum_{i,j}\phi^{*}_{i}A_{ij}\phi_{j}},$$

$$\int\bigl{[}d\phi^{*}d\phi\bigr{]}=\int\bigl{[}d\xi^{*}d\xi\bigr{]}\int\bigl{[}d\eta^{*}d\eta\bigr{]}.$$

$$e^{{\varphi }^{\ast }a\varphi }=e^{{\eta }^{\ast }{\xi }^{\ast }a\xi \eta }=\frac{1}{2}\sum _{s=\pm 1}{e}^{s\sqrt{a}{\xi }^{\ast }\xi }{e}^{s\sqrt{a}{\eta }^{\ast }\eta },$$

$$e^{{\sum }_{i,j}{\varphi }_i^{\ast }{A}_{ij}{\varphi }_j}=\frac{1}{2^m}\sum _{\{S\}}{e}^{{\sum }_{ij}{\xi }_i^{\ast }{G}_{ij}(S){\xi }_j}{e}^{{\sum }_{ij}{\eta }_i^{\ast }{G}_{ij}(S){\eta }_j},$$

$$G_{ij}(S)=\sqrt{A_{ij}}{s}_{ij},$$

$$\operatorname{\text{per}}A=\frac{1}{2^{m}}\sum_{\{S\}}\bigl{(}\det G(S)\bigr{)}^{2}.$$

$$e^{{\varphi }^{\ast }a\varphi }=\frac{1}{p}\sum _{q=1}^p{e}^{{\omega }^q\sqrt{a}{\xi }^{\ast }\xi +{\omega }^{-q}\sqrt{a}{\eta }^{\ast }\eta },$$

$$H_{ij}(Q)=\sqrt{A_{ij}}{\omega }^{q_{ij}},$$

$$\operatorname{\text{per}}A=\frac{1}{p^{m}}\sum_{\{Q\}}\bigl{|}\det H(Q)\bigr{|}^{2},$$

$$A=\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right),H(S)=\left(\begin{array}{cc}{s}_{11}& i{s}_{12}\\ s_{21}& s_{22}\end{array}\right).$$

$$e^{{\varphi }^{\ast }a\varphi }=\frac{1}{p}\sum _{q=1}^p{e}^{{\omega }^q\sqrt{a}{\varphi }^{\ast }}{e}^{{\omega }^{-q}\sqrt{a}\varphi }.$$

$$\operatorname{\text{per}}A=\frac{1}{p^{m}}\sum_{\{Q\}}\\ \times\prod_{i}\Bigl{(}\sum_{k_{i}}\sqrt{A_{k_{i}i}}\omega^{-q_{k_{i}i}}\Bigr{)}\Bigl{(}\sum_{j_{i}}\sqrt{A_{ij_{i}}}\omega^{q_{ij_{i}}}\Bigr{)}.$$

$$\operatorname{\text{per}}A=\frac{1}{2^{m}}\sum_{\{S\}}\prod_{i}\Bigl{(}\sum_{k_{i}}s_{k_{i}i}\Bigr{)}\Bigl{(}\sum_{j_{i}}s_{ij_{i}}\Bigr{)}.$$

$$e^{\sum_{i,j}\phi^{*}_{i}A_{ij}\phi_{j}}=\prod_{i}\Bigl{(}1+\phi^{*}_{i}\sum_{j_{i}}A_{ij_{i}}\phi_{j_{i}}\Bigr{)}.$$

$$\operatorname{\text{per}}A=\int\bigl{[}d\phi\bigr{]}\prod_{i}\Bigl{(}\sum_{j_{i}}A_{ij_{i}}\phi_{j_{i}}\Bigr{)}.$$

$$\operatorname{\text{per}}A=\frac{1}{2^{n}}\sum_{\{S\}}\Bigl{(}\prod_{k}s_{k}\Bigr{)}\int\bigl{[}d\phi\bigr{]}\prod_{i}\Bigl{(}\sum_{j_{i}}A_{ij_{i}}s_{j_{i}}\phi_{j_{i}}\Bigr{)}.$$

$$\operatorname{\text{per}}A=\frac{1}{2^{n}}\sum_{\{S\}}\Bigl{(}\prod_{k}s_{k}\Bigr{)}\prod_{i}\Bigl{(}\sum_{j_{i}}A_{ij_{i}}s_{j_{i}}\Bigr{)}.$$

$$\displaystyle\operatorname{\text{per}}A=\frac{1}{2^{n-1}}{\sum_{\{S\}}}^{\prime}\Bigl{(}\prod_{k}s_{k}\Bigr{)}\prod_{i}\Bigl{(}\sum_{j_{i}}A_{ij_{i}}s_{j_{i}}\Bigr{)},$$

$$\operatorname{\text{per}}A=\frac{1}{p^{n}}\sum_{\{Q\}}\Bigl{(}\prod_{k}\omega^{-q_{k}}\Bigr{)}\prod_{i}\Bigl{(}\sum_{j_{i}}A_{ij_{i}}\omega^{q_{j_{i}}}\Bigr{)}.$$

$$e^{\sum_{i,j,k}\phi_{i}^{*}L_{ij}U_{jk}\phi_{k}}=\int\bigl{[}d\varphi^{*}d\varphi\bigr{]}e^{\sum_{i,j}(\phi_{i}^{*}L_{ij}\varphi_{j}+\varphi_{j}^{*}U_{ji}\phi_{i})},$$

$$\int\bigl{[}d\varphi^{*}d\varphi\bigr{]}=\int\prod_{j}\frac{d\varphi_{j}^{*}d\varphi_{j}}{2\pi i}e^{-|\varphi_{j}|^{2}}.$$

$$\operatorname{\text{per}}A=\int\bigl{[}d\varphi^{*}d\varphi\bigr{]}\prod_{i}\Bigl{(}\sum_{k_{i}}\varphi^{*}_{k_{i}}U_{k_{i}i}\Bigr{)}\Bigl{(}\sum_{j_{i}}L_{ij_{i}}\varphi_{j_{i}}\Bigr{)}.$$

$$\operatorname{\text{per}}A=\int\bigl{[}d\varphi^{*}d\varphi\bigr{]}\\ \times\prod_{i}\Bigl{(}\sum_{j_{i}}U_{ij_{i}}\sqrt{\sigma_{j_{i}}}\varphi_{j_{i}}\Bigr{)}\Bigl{(}\sum_{k_{i}}V_{ik_{i}}\sqrt{\sigma_{k_{i}}}\varphi^{*}_{k_{i}}\Bigr{)}.$$

$$\int\bigl{[}d\varphi^{*}d\varphi\bigr{]}\varphi^{*i}\varphi^{j}=\delta_{ij}j!\qquad\text{for }0\leq i,j\leq n.$$

$$e^{{\varphi }^{\ast }a\varphi }=\frac{1}{2}\sum _{s=\pm 1}{e}^{s\sqrt{a}{\mu }^{\ast }\mu }{e}^{s\sqrt{a}{\nu }^{\ast }\nu },$$

$$\operatorname{\text{per}}A=\frac{1}{2^{m}}\sum_{\{S\}}\bigl{(}\operatorname{\text{per}}G(S)\bigr{)}^{2}.$$

Full text

The permanent of an $n\times n$ matrix $A$ is defined by

[TABLE]

where the sum is over all permutations $\mathcal{P}$ of $1,2,\ldots,n$. The permanent is therefore somewhat similar to the determinant, but the signature of the permutation is absent in the formula. In physics the permanent appears for example when calculating the overlap of two bosonic many-body wave functions Negele and Orland (1988). It is widely appreciated that it is much more difficult to evaluate the permanent than the determinant. This was put on firm ground in the context of complexity theory when Valiant proved that even the computation of the permanent of a matrix with $0\,$-$1$ valued entries is in the class of $\#P$-complete problems Valiant (1979). On a less sophisticated level one important reason for this complexity is that the permanent is not invariant under similarity transformations, whereas the determinant is.

Grassmann integrals are commonly used as a device in many-body quantum mechanics to study interacting fermions in the path integral formalism. They were introduced by Berezin Berezin (1965), and are since long textbook material, see e.g. Negele and Orland (1988). The Gaussian Grassmann integral is of particular importance and results in a determinant

[TABLE]

The permanent of a matrix can also be expressed as a Gaussian Grassmann integral

[TABLE]

over even Grassmann numbers that satisfy $\phi_{i}\phi_{j}=\phi_{j}\phi_{i}$ and $\phi_{i}^{2}=0$. This simple formula is one of crucial importance for this letter. It has appeared before in the physics literature Palumbo (1994a, b), but is apparently not well-known. The even Grassmann numbers have also been introduced in the applied mathematics and computer science communities where they go under the name of “zeons”, see e.g. Schott and Staples (2012). The proof of (3) is easily obtained by direct expansion sup . Since different $\phi_{i}$ commute no sign is generated upon rearranging the $\phi$’s into canonical order after expanding the exponential. The final result is therefore the permanent instead of the determinant. For our purposes it is important to note that these Grassmann numbers can also be viewed as composite objects Palumbo (1994a, b); sup , defined as products of two ordinary anti-commuting Grassmann numbers $\xi_{i}$ and $\eta_{i}$: $\phi_{i}=\xi_{i}\eta_{i}$. We will also use that the integration measure for $\phi$ may be separated into a product of two independent ones over $\xi$ and $\eta$

[TABLE]

In the language of physics integrals of the type (3) can describe an interacting fermion theory of doublons in the Hubbard model, and lattice models involving spin-1/2 objects and hard-core bosons nil . In this case the matrix $A$ is of large dimension but has a specific block structure and is very sparse.

.1 Formulas from discrete Hubbard-Stratonovich transformations

The Hubbard-Stratonovich (HS) transformation is another standard method used in many-body physics to reformulate an interacting theory in terms of a weighted sum of non-interacting ones Stratonovich (1957); Hubbard (1959). This transformation is important both in analytical approaches as well as in the construction of the determinant quantum Monte Carlo formalism Blankenbecler et al. (1981). For our Grassmann numbers we will first use a discrete variant of this transformation which is the simple identity

[TABLE]

which is valid for any real or complex number $a$. Similar discrete HS transformations can also be introduced at the operator level Hirsch (1983), but the Grassmann version is almost trivial in comparison. Introducing a sign $s_{ij}$ for each of the $m$ non-vanishing matrix elements of $A_{ij}$ we may rewrite the exponential using

[TABLE]

where the matrix elements of $G(S)$ are

[TABLE]

and $S$ is a shorthand for the particular configuration of all of the $s_{ij}$. Using the behavior of the measure (4) we can now perform the Gaussian Grassmann integrals over $\xi$ and $\eta$ independently to get

[TABLE]

This is the unbiased Godsil-Gutman estimator for the permanent Godsil and Gutman (1978). This estimator unfortunately has a very large variance in the worst case, and does therefore not work well for all matrices.

The Godsil-Gutman estimator has been generalized in many different ways to reduce its variance. In the physics language some such generalizations may be generated by using other HS transformations than the $\mathbb{Z}_{2}$ one of (5). In particular we could instead use a $\mathbb{Z}_{p}$ decoupling

[TABLE]

where $\omega=e^{2\pi i/p}$ is a $p\,$-th root of unity. Note that (5) is a special case of this corresponding to $p=2$. Clearly this can again be done independently on each non-vanishing element of $A_{ij}$. Introducing $q_{ij}$ and $H(Q)$ in analogy with $s_{ij}$ and $G(S)$ above with

[TABLE]

we may perform the Grassmann integrals with the result

[TABLE]

when $A$ has only positive semi-definite elements. This includes (setting $p=3$) the KKLLL esimator of Karmarkar et al. Karmarkar et al. (1993). The bounds on the variance of this estimator is much better than that of Godsil-Gutman. The variance is however still extremely large in the worst case Karmarkar et al. (1993). On the other hand it has been proven that the KKLLL estimator is very efficient for sufficiently dense matrices Frieze and Jerrum (1995).

In the physics context it is well-known that it is possible to decouple interaction terms in different channels Negele and Orland (1988). This leads to additional formulas for the permanent. So far we have only considered decouplings that preserves “spin rotational symmetry” around one axis, meaning that we have only made use of bilinears of the types $\xi^{*}\xi$ and $\eta^{*}\eta$. In (5) we only use the density channel, meaning that the bilinears are of the form $\xi^{*}\xi+\eta^{*}\eta$. It is also possible to do decouplings in the spin channel, the simplest of these involves $i(\xi^{*}\xi-\eta^{*}\eta)$, which do appear in (9). We may also decouple using the bilinears $i(\xi^{*}\eta+\eta^{*}\xi)$ and $(\xi^{*}\eta-\eta^{*}\xi)$, which breaks the spin rotational symmetry. These are all examples of time-reversal invariant decoupling schemes Wu and Zhang (2005), which have the appealing property that the weights are always positive semi-definite when the elements of $A$ are. In fact it is also possible to use different forms of the decoupling on different matrix elements. A simple example demonstrates that using this freedom may be extremely fruitful. Consider the following matrix $A$ and one particular decoupled version $H(S)$

[TABLE]

The corresponding estimator has zero variance! Of course it is a very difficult problem to pick the “best” decoupling scheme for a given large matrix. In the physics context the choice of decoupling is supposed to be “guided by the physics of the problem” Negele and Orland (1988), something that seems quite difficult in the abstract setting of a generic matrix $A$.

It is also possible to perform a decoupling in the pairing channel, which is equivalent to a decoupling in $\phi^{*}$ and $\phi$. Let us use a general $\mathbb{Z}_{p}$ decoupling ($p\geq 2)$

[TABLE]

Implementing this on each non-vanishing term in $A$ as above and performing the Grassmann integral we obtain

[TABLE]

Specializing to matrices with $0\,$-$1$ valued entries and taking $p=2$ we get

[TABLE]

This formula treats rows and columns symmetrically and has a simple interpretation: all non-zero elements of $A$ are substituted with $\pm 1$ with equal probability, and the contribution from a configuration is the product of all row sums and columns sums.

.2 Efficient exact formulas

Efficient exact formulas may also be easily obtained from manipulations of the Grassmann integral. Let us first note that because $\phi_{i}^{*2}=0$ the Gaussian exponent may be expanded as

[TABLE]

Sticking this into the formula for the permanent and performing the integrals over the $\phi_{i}^{*}$’s we get

[TABLE]

Now let us attach a local $\mathbb{Z}_{2}$ gauge freedom to each $\phi_{j}$. We implement this by making a change of variables inside the integral, taking $\phi_{j}\rightarrow s_{j}\phi_{j}$ with $s_{j}=\pm 1$. The Jacobian generated by going to the new measure is simply $\prod_{k}s_{k}$. Since nothing can depend on the variable change, and hence the $s_{i}$’s, we may average over all possibilities

[TABLE]

In this formula the role of the Grassmann integral, i.e., to pick out the terms for which each $\phi_{i}$ appears exactly once, is superfluous since terms in the integrand without this property are anyway set to zero by the sum over $S$. Getting rid of the integral we therefore get

[TABLE]

This is an exact formula that is a sum of $2^{n}$ terms. A more efficient albeit less symmetric formula involving a sum of $2^{n-1}$ terms may easily be obtained from this one by noting that each term is invariant upon inverting all signs. This implies that for each configuration with say $s_{1}=-1$ there is another one with $s_{1}=1$ with equal weight. We may therefore write

[TABLE]

where the primed sum is over all sign configurations but with $s_{1}=1$ always. This is the formula of Glynn Glynn (2010). Glynn has also formulated this in terms of polarization identities Glynn (2013), this makes clear the connection to the Ryser formula Ryser (1963), which is of comparable efficiency. It is interesting to note that in the physics language (19) is a result of local $\mathbb{Z}_{2}$ gauge invariance, and the reduction to (20) a consequence of global $\mathbb{Z}_{2}$ gauge invariance.

.3 Additional formulas

The derivations above motivates us to consider a few generalizations that gives some additional formulas for the permanent.

It is clear that the formulas (19) and (20) may also be considered as unbiased estimators for the permanent when the signs are treated as random variables Aaronson and Hance (2014). We now note that in the discrete HS case the variance was reduced in going from Godsil-Gutman to KKLLL, which may be viewed as changing the decoupling variable from $\mathbb{Z}_{2}$ to $\mathbb{Z}_{p}$ in our language. In the derivation of the exact formula above we may easily implement the same idea, i.e., we make the gauge transformation (variable change) $\phi_{i}\rightarrow\omega^{q_{i}}\phi_{i}$. This leads to the formula

[TABLE]

Viewed as an exact formula this is obviously less efficient than (19) since the number of terms is $p^{n}$ instead of $2^{n}$. From the point of view of an unbiased estimator the variance is however substantially reduced for certain matrices. A common example is the $n\times n$ block-diagonal matrix with the matrix $A$ in (12) repeated on the diagonal. For such a matrix the second moment is reduced by a factor of $(3/4)^{n/2}$ in going from (19) to (21) for all $p\geq 3$. This is similar to the reduction in going from Godsil-Gutman to KKLLL Karmarkar et al. (1993); Chien et al. (2003).

It is also possible to consider various continuous decoupling schemes. Let us decompose our matrix as $A=LU\negthinspace P$ where $P$ is a permutation matrix and $L$ ($U$) are lower (upper) triangular matrices (such a decomposition is always possible). $P$ does not affect the permanent and we only need to consider the permanent of the matrix $LU$. Now we may perform a conventional continuous HS transformation on the exponent as follows

[TABLE]

with a normalized Gaussian integration measure

[TABLE]

Sticking this into the formula for the permanent it is easy to perform the Grassmann integral with the result

[TABLE]

This formula is a $2n$-dimensional integral representation of the permanent. Similar formulas also works for other types of matrix decompositions. Consider for example a singular value decomposition $A=U\Sigma V^{T}$, with $U$ and $V$ orthogonal and $\Sigma$ diagonal. Let us denote the diagonal elements of $\Sigma$ by $\sigma_{i}$, then the same manipulations gives

[TABLE]

These $2n$-dimensional integral representations involves an unbounded integration region. It is also possible to reformulate these on bounded integration regions, what is needed is a weight function that satisfy

[TABLE]

The integrals (24) and (25) may be estimated using stochastic Monte Carlo methods. The last version is particularly efficient for matrices of low rank since the number of necessary $\varphi$’s is equal to the rank of the matrix.

Let us finally mention another possibility that may be obtained from a slight extension of our formalism. Suppose that we represent $\phi=\mu\nu$, with $\mu$ and $\nu$ a pair of even Grassmann numbers. We may then decouple a term in the exponent as before

[TABLE]

and the measure still factorizes in analogy with (4). Doing this for each of the $m$ non-vanishing elements of $A$ like we did to get to (8) and using (3) we arrive at the formula

[TABLE]

This formula should be compared with (8) and can be used recursively to generate additional formulas.

Acknowledgements.

.4 Acknowledgements

Funding from the Knut and Alice Wallenberg Foundation and the Swedish research council Vetenskapsrådet is gratefully acknowledged.