A SIMPLE POLYNOMIAL TIME ALGORITHM TO APPROXIMATE THE PERMANENT WITHIN A SIMPLY EXPONENTIAL FACTOR

ALEXANDER BARVINOK

Abstract. We present a simple randomized polynomial time algorithm to approximate the mixed discriminant of n positive semidefinite n×n matrices within a factor 2^{O(n)}. Consequently, the algorithm allows us to approximate in randomized polynomial time the permanent of a given n×n non-negative matrix within a factor 2^{O(n)}. When applied to approximating the permanent, the algorithm turns out to be a simple modification of the well-known Godsil-Gutman estimator.

1. Introduction

In this paper we address the question of how to approximate the permanent of a non-negative matrix, and, more generally, the mixed discriminant of positive semidefinite matrices. Our main result is that a simple modification of the well-known Godsil-Gutman estimator ([10], see also Chapter 8 of [19]) yields a randomized polynomial time algorithm, which, given an n × n non-negative matrix, approximates its permanent within a 2^{O(n)} factor. It turns out that the ideas of the algorithm and the proof become more transparent when generalized to mixed discriminants. The first randomized polynomial time algorithm that approximates the permanent within a 2^{O(n)} factor was suggested by the author in [5]. The algorithm described in this paper has some advantages compared to the algorithm from [5]. It has a much more transparent structure, it is much easier to implement and it is easily parallelizable. Besides, it sheds some additional light on the properties of the Godsil-Gutman estimator, which was studied in several papers (see [9], [16] and Chapter 8 of [19]).

(1.1) Permanent. Let A = (a_{ij}) be an n × n matrix and let S_{n} be the symmetric group, that is the group of all permutations of the set {1, …, n}. The number

per A = ∑_{σ∈S_{n}} ∏_{i=1}^{n} a_{iσ(i)}

is called the permanent of A. If A is a 0-1 matrix, then per A can be interpreted as the number of perfect matchings in a bipartite graph G on 2n vertices v_{1}, …, v_{n} and u_{1}, …, u_{n}, where (v_{i}, u_{j}) is an edge of G if and only if a_{ij} = 1. To compute the permanent of a given 0-1 matrix is a #P-complete problem and even to estimate per A seems to be difficult. Polynomial time algorithms for computing per A are known when A has some special structure, for example, when A is sparse [11] or has small rank [4]. Polynomial time approximation schemes are known for dense 0-1 matrices [12], for “almost all” 0-1 matrices (see [12], [21] and [9]) and for some special 0-1 matrices, such as corresponding to lattice graphs, (see [13] for a survey on approximation algorithms). However, not much is known on how to approximate in polynomial time the permanent of an arbitrary 0-1 matrix (see [15] for the fastest known “mildly exponential” approximation scheme), let alone the permanent of an arbitrary non-negative matrix.

Key words and phrases. permanent, mixed discriminant, randomized algorithms, approximation algorithms.

Let t_{1}, …, t_{n} be real variables. Then the permanent per A can be expressed as the coefficient of t_{1} ⋯ t_{n} in the product of linear forms:

(1.1.1)  per A = \frac{∂^{n}}{∂t_{1} … ∂t_{n}} ∏_{i=1}^{n} ∑_{j=1}^{n} a_{ij}t_{j}.

(1.2) Mixed Discriminant. Let Q_{1}, …, Q_{n} be n × n symmetric matrices and let t_{1}, …, t_{n} be real variables. Then det(t_{1}Q_{1} +…+t_{n}Q_{n}) is a homogeneous polynomial of degree n in t_{1}, …, t_{n}. The number

D(Q_{1}, …, Q_{n}) = \frac{∂^{n}}{∂t_{1} … ∂t_{n}} det(t_{1}Q_{1} + … + t_{n}Q_{n})

is called the mixed discriminant of Q_{1}, …, Q_{n}. Sometimes the normalizing factor 1/n! is used (cf. [18]). The mixed discriminant D(Q_{1}, …, Q_{n}) is a polynomial in the entries of Q_{1}, …, Q_{n} with coefficients −1 and 1.

The mixed discriminant can be considered as a generalization of the permanent. Indeed, from (1.1.1) we deduce that for diagonal matrices Q_{1}, …, Q_{n} (1.2.1)

D(Q_{1}, …, Q_{n}) = per A, where Q_{i} = diag{a_{i1}, …, a_{in}} and A = (a_{ij}).

Mixed discriminants were introduced by A.D. Aleksandrov in his proof of the Aleksandrov - Fenchel inequality for mixed volumes ([2], see also [18]). The relation between the mixed discriminant and the permanent was used in the proof of the Van der Waerden conjecture for permanents of doubly stochastic matrices (see [6]).

The mixed discriminant is linear in each argument, that is

D(Q_{1}, …, Q_{i−1}, αQ_{i} + βQ′_{i}, Q_{i+1}, …, Q_{n})

= αD(Q_{1}, …, Q_{i−1}, Q_{i}, Q_{i+1}, …, Q_{n}) + βD(Q_{1}, …, Q_{i−1}, Q′_{i}, Q_{i+1}, …, Q_{n}), and D(Q_{1}, …, Q_{n}) ≥ 0 provided Q_{1}, …, Q_{n} are positive semidefinite (see [18]).

The paper is organized as follows. In Section 2, we describe the algorithms for computing mixed discriminants and permanents, state the main results, and discuss them. Section 3 addresses the behavior of quadratic forms on random vectors drawn from the Gaussian distribution in ℝ^{n} – our main tool for analyzing the algorithms. Section 4 contains proofs of the main results of the paper. In Section 5, we present an application of the mixed discriminant for counting problems and also make some remarks on “binary versions” of our algorithms.

2. Main Results

Our algorithms use two standard procedures from linear algebra: computing the determinant of a matrix and computing a decomposition Q = T T^{∗}, where Q is a positive semidefinite matrix and T^{∗} is the transpose of T. One can compute the determinant of a given n × n matrix using O(n^{3}) arithmetic operations (see, for example, Chapter 2, Section 7 of [8]). For a positive semidefinite n × n matrix Q one can compute a decomposition Q = T T^{∗}, where T is a lower triangular matrix, by using O(n^{3}) arithmetic operations and n times taking square root from a nonnegative number (see, for example, Chapter 2, Section 10 of [8]). The “random” part of the algorithms consists of sampling n^{2} random variables independently from the standard Gaussian distribution in ℝ with the density

ψ_{1}(x) = \frac{1}{\sqrt{2π}} e^{−x^{2}/2}.

Various ways to simulate this distribution from the uniform distribution on the interval [0, 1] are described in Section 3.4.1 (C) of [17]. In particular, this allows us to sample vectors x = (x_{1}, …, x_{n}) ∈ ℝ^{n} from the standard n-dimensional Gaussian distribution in ℝ^{n} with the density

ψ_{n}(x) = (2π)^{−n/2} exp{−∥x∥^{2}/2}, where ∥x∥^{2} = x_{1}^{2} + … + x_{n}^{2},

by sampling the coordinates x_{1}, …, x_{n} independently from the one-dimensional Gaussian distribution. We will write ψ(x) instead of ψ_{n}(x) if the choice of the ambient space ℝ^{n} is clear from the context.

To simplify our analysis, we assume that we operate with real numbers and that we can perform arithmetic operations, take square root and sample from the Gaussian distribution exactly. In Section 5.1 we briefly discuss how to adjust our algorithms to the “binary model” of computation, that is, if we allow only integers (rational numbers) and arithmetic operations on them. As a general source on algorithms and complexity we use [1].

Now we describe our main algorithm.

(2.1) Algorithm for computing the mixed discriminant. Input: Positive semidefinite n × n matrices Q_{1}, …, Q_{n}. Output: A number α approximating the mixed discriminant D(Q_{1}, …, Q_{n}).

Algorithm: For i = 1, …, n compute a decomposition Q_{i} = T_{i}T_{i}^{∗}. Sample independently n vectors u_{1}, …, u_{n} at random from the standard Gaussian distribution in ℝ^{n} with the density ψ(x) = (2π)^{−n/2} exp{−∥x∥^{2}/2}. Compute

α = (det[T_{1}u_{1}, …, T_{n}u_{n}])^{2},

the squared determinant of the matrix with the columns T_{1}u_{1}, …, T_{n}u_{n}. Output α.

(2.2) Theorem. (2.2.1) The expectation of α is the mixed discriminant D(Q_{1}, …, Q_{n}); (2.2.2) For any C > 1 the probability that

α ≥ C · D(Q_{1}, …, Q_{n})

does not exceed C^{−1}; (2.2.3) Let

c_{0} = exp\{\frac{4}{\sqrt{2π}} ∫_{0}^{+∞}(ln t)e^{−t^{2}/2} dt\} ≈ 0.2807297419.

Then for any 1 > ϵ > 0 the probability that

α ≤ (ϵc_{0})^{n}D(Q_{1}, …, Q_{n})

does not exceed \frac{8}{n ln^{2} ϵ}.

The algorithm becomes especially simple when we apply it to compute the mixed discriminant of diagonal matrices, that is the permanent of a non-negative matrix.

(2.3) Algorithm for computing the permanent of a non-negative matrix. Input: Non-negative n × n matrix A. Output: A number α approximating the permanent per A. Algorithm: Sample independently n^{2} numbers u_{ij} : i, j = 1, …, n at random

from the standard Gaussian distribution in ℝ with the density ψ(x) = \frac{1}{\sqrt{2π}} e^{−x^{2}/2}.

Compute the n × n matrix B = (b_{ij}), where b_{ij} = u_{ij}\sqrt{a_{ij}}. Compute α = (det B)^{2}. Output α.

(2.4) Theorem. (2.4.1) The expectation of α is the permanent per A; (2.4.2) For any C > 1 the probability that

α ≥ C · per A

does not exceed C^{−1}; (2.4.3) Let

c_{0} = exp\{\frac{4}{\sqrt{2π}} ∫_{0}^{+∞}(ln t)e^{−t^{2}/2} dt\} ≈ 0.2807297419.

Then for any 1 > ϵ > 0 the probability that

α ≤ (ϵc_{0})^{n}per A

does not exceed \frac{8}{n ln^{2} ϵ}.

As implied by (1.2.1), Theorem 2.4 is a straightforward corollary of Theorem 2.2. From (2.2.2) and (2.2.3) it follows that for all sufficiently large n, the output α of Algorithm 2.1 satisfies the inequalities

(0.28)^{n}D(Q_{1}, …, Q_{n}) ≤ α ≤ 3D(Q_{1}, …, Q_{n})

with probability at least 0.6. We can make the probability as high as 1−ϵ for any ϵ > 0 by running the algorithm independently O(ln(ϵ^{−1})) times and taking the median of the computed α’s (cf. [14]). Hence we get a randomized polynomial time algorithm approximating the mixed discriminant of positive semidefinite matrices (and hence the permanent of a non-negative matrix) within a simply exponential factor 2^{O(n)}.

(2.5) Relation to the Godsil-Gutman estimator. It is seen immediately that Algorithm 2.3 is a modification of the Godsil-Gutman estimator (see [10]) and Chapter 8 of [19]). Indeed, in the Godsil-Gutman estimator we sample u_{ij} from the binary distribution:

u_{ij} = \{ \begin{matrix}1 & with probability 1/2, \\ −1 & with probability 1/2. \end{matrix}

Furthermore, parts (2.4.1) and (2.4.2) of Theorem 2.4 remain true as long as we sample u_{ij} independently from some distribution with the expectation 0 and variance 1. However, (2.4.3) is not true for the binary distribution. Indeed, let A = (a_{ij}) be the following n × n matrix

a_{ij} = \{ \begin{matrix}1 & if i = j or \{i, j\} = \{1, 2\}, \\ 0 & elsewhere. \end{matrix}

Then per A = 2. If u_{ij} ∈ {−1, 1} are sampled from the binary distribution, then α = (u_{11}u_{22} − u_{12}u_{21})^{2}. So, we get that α = 0 with probability 1/2 and α = 4 with probability 1/2. Apparently, sampling from continuous distributions allows us to approximate better.

Another (discrete) distribution for the Godsil-Gutman estimator was studied in [16] and [9]. It was shown that for a “typical” 0-1 matrix we get a polynomial approximation scheme for the permanent [9], whereas for a “worst case” 0-1 matrix it allows us to construct an approximation scheme [16] whose complexity, while still exponential, is significantly better than the complexity of known exact methods (cf. Chapter 7 of [20]).

(2.6) How well can we approximate the permanent in polynomial time? Suppose we have a polynomial time (probabilistic or deterministic) algorithm that for any given n × n (non-negative or 0-1) matrix A computes a number α such that ϕ(n)per A ≤ α ≤ per A. What sort of function might ϕ(n) be? For an n × n matrix A and k > 0 let us construct the nk × nk block-diagonal matrix A ⊗ I_{k}, having k diagonal copies of A. We observe that A⊗I_{k} is non-negative if A is non-negative and A⊗I_{k} is 0-1 if A is 0-1, and that per A = (per A⊗I_{k})^{1/k}. Applying our algorithm to A⊗I_{k} and taking the root we get an approximation of per A within a factor ϕ^{1/k}(nk). So, if we suppose that ϕ is the best possible, we may assume that ϕ(n) ≥ ϕ^{1/k}(nk). There are few reasonable choices for such functions ϕ.

a) ϕ(n) ≡ 1. This doesn’t look likely, given that the problem is #P-hard.

b) For ϵ > 0 one can choose ϕ_{ϵ}(n) = 1 − ϵ and the algorithm is polynomial in ϵ^{−1}. In the author’s opinion this conjecture is overly optimistic, especially for arbitrary non-negative matrices.

c) For ϵ > 0 one can choose ϕ_{ϵ}(n) = (1 − ϵ)^{n}, but the algorithm is not polynomial in ϵ^{−1}. This type of approximation was conjectured by V.D. Milman.

d) ϕ(n) = c^{n} for some fixed constant c. This is the type of a bound achieved by Algorithm 2.3.

We note that functions ϕ(n), decreasing faster than a simply exponential function of n (for example, ϕ(n) = 1/n!) are not interesting since they are beaten by c^{n} of Algorithm 2.3. The author does not know if the constant c ≈ 0.28 from Theorem 2.4 can be improved.

3. The Gaussian Measure and Quadratic Forms

(3.1) Notation. Let us fix the standard orthonormal basis in ℝ^{n}. Thus we can identify n × n matrices with linear operators on ℝ^{n}.

Let ⟨·, ·⟩ be the standard inner product on ℝ^{n}, that is

⟨x, y⟩ = x_{1}y_{1} + … + x_{n}y_{n}, where x = (x_{1}, …, x_{n}) and y = (y_{1}, …, y_{n}).

We denote the corresponding Euclidean norm by ∥ · ∥.

Let μ be the standard Gaussian measure on ℝ^{n} with the density

ψ(x) = (2π)^{−n/2} exp{−∥x∥^{2}/2}.

For a μ-integrable function f : ℝ^{n} ⟶ ℝ we define its expectation

E(f) = ∫_{ℝ^{n}} f dμ = ∫_{ℝ^{n}} f(x)ψ(x) dx.

Furthermore, if F : ℝ^{n} ⟶ ℝ^{m}, F (x) = (f_{1}(x), …, f_{m}(x)), we let E(F) = (E(f_{1}), …, E(f_{m})) ∈ ℝ^{m}.

If x ∈ ℝ^{n} is a vector then we denote by x ⊗ x the n × n matrix whose (i, j)-th entry is x_{i} · x_{j}. So, x ⊗ x is a positive semidefinite matrix of rank 1, provided x ≠ 0.

Let Q be a symmetric n × n matrix and q(x) = ⟨x, Qx⟩ be the corresponding quadratic form. We define the trace of q:

tr(q) = ∑_{i=1}^{n} Q_{ii}.

We start with a simple lemma.

(3.2) Lemma. (3.2.1) Let q : ℝ^{n} ⟶ ℝ be a quadratic form. Then

E(q) = tr(q).

(3.2.2) Let us fix an n × n matrix T. Then

E(T u ⊗ T u) = T T^{∗},

where u ∈ ℝ^{n} is sampled from the standard Gaussian distribution μ.

Proof. Since E(x_{i}^{2}) = 1 for i = 1, … n and E(x_{i}x_{j}) = 0, the proof of (3.2.1) follows. Therefore, E(u ⊗ u) = I, the identity matrix, and hence

E(T u ⊗ T u) = E(T (u ⊗ u)T^{∗}) = T E(u ⊗ u)T^{∗} = T T^{∗},

and (3.2.2) follows.

Now we prove the main result of this section.

(3.3) Theorem. Let q : ℝ^{n} ⟶ ℝ be a positive semidefinite quadratic form, such that E(q) = 1. Let

C_{0} = \frac{4}{\sqrt{2π}} ∫_{0}^{+∞}(ln t)e^{−t^{2}/2} dt ≈ −1.270362845.

Then

(3.3.1)  C_{0} ≤ E(ln q) ≤ 0

and

(3.3.2)  0 ≤ E(ln^{2} q) ≤ 8.

Proof. Since ln is a concave function, we have E(ln q) ≤ ln(E(q)) = 0. Let us decompose q into a non-negative linear combination q = λ_{1}q_{1} + … + λ_{n}q_{n} of positive semidefinite forms q_{i} of rank 1. We can scale q_{i} so that E(q_{i}) = 1 for i = 1, …, n and then we have λ_{1} + … + λ_{n} = 1. In fact, one can choose λ_{i} to be the eigenvalues of q and q_{i} = ⟨x, u_{i}⟩^{2}, where u_{i} is the corresponding unit eigenvector. Since ln is a concave function we have ln(λ_{1}q_{1} + … + λ_{n}q_{n}) ≥ λ_{1} ln q_{1} + … + λ_{n} ln q_{n}. Furthermore, if q_{i} is a positive semidefinite form of rank 1 such that E(q_{i}) = 1, then by an orthogonal transformation of the coordinates it can be brought into the form q_{i}(x) = x_{1}^{2} (see (3.2.1)). Therefore, E(ln q_{i}) = E(ln x_{1}^{2}) = C_{0} and

E(ln q) ≥ λ_{1}E(ln q_{1}) + … + λ_{n}E(ln q_{n}) = (λ_{1} + … + λ_{n})C_{0} = C_{0},

so (3.3.1) is proven (we note that this reasoning proves that E(ln q) is well-defined). Let X = {x ∈ ℝ^{n} : q(x) ≤ 1} and Y = ℝ^{n} \ X. Then

E(ln^{2} q) = ∫_{X} ψ(x) ln^{2} q(x) dx + ∫_{Y} ψ(x) ln^{2} q(x) dx.

Let us estimate the first integral. Decomposing q = λ_{1}q_{1} + … + λ_{n}q_{n} as above, we get ln q ≥ λ_{1} ln q_{1} + … + λ_{n} ln q_{n}. Since ln q(x) ≤ 0 for x ∈ X, we get that

ln^{2} q(x) ≤ ∑_{i,j=1}^{n} λ_{i}λ_{j}(ln q_{i}(x))(ln q_{j}(x))

for x ∈ X. Therefore,

∫_{X} ψ(x) ln^{2} q(x) dx ≤ ∑_{i,j=1}^{n} λ_{i}λ_{j} ∫_{X} ψ(x)(ln q_{i}(x))(ln q_{j}(x)) dx

≤ ∑_{i,j=1}^{n} λ_{i}λ_{j}(∫_{X} ψ(x) ln^{2} q_{i}(x) dx)^{1/2}(∫_{X} ψ(x) ln^{2} q_{j}(x) dx)^{1/2}

(we applied the Cauchy-Schwartz inequality)

≤ ∑_{i,j=1}^{n} λ_{i}λ_{j}(E(ln^{2} q_{i}))^{1/2}(E(ln^{2} q_{j}))^{1/2}.

Now, as in the proof of (3.3.1) we have

E(ln^{2} q_{i}) = E(ln^{2} x_{1}^{2}) = \frac{8}{\sqrt{2π}} ∫_{0}^{+∞}(ln^{2} t)e^{−t^{2}/2} dt ≈ 6.548623960 ≤ 7.

Summarizing, we get

∫_{X} ψ(x) ln^{2} q(x) dx ≤ ∑_{i,j=1}^{n} λ_{i}λ_{j}(E(ln^{2} q_{i}))^{1/2}(E(ln^{2} q_{j}))^{1/2} ≤ 7 ∑_{i,j=1}^{n} λ_{i}λ_{j} = 7.

Since for 0 ≤ ln t ≤ \sqrt{t} for t ≥ 1 we have

∫_{Y} ψ(x) ln^{2} q(x) dx ≤ ∫_{Y} q(x)ψ(x) dx ≤ E(q) = 1.

Therefore, E(ln^{2} q) ≤ 7 + 1 = 8 and (3.3.2) is proven.

Finally, we will need the following simple result.

(3.4) Lemma. Let u_{1}, …, u_{n} be vectors from ℝ^{n}. Then

D(u_{1} ⊗ u_{1}, …, u_{n} ⊗ u_{n}) = (det[u_{1}, …, u_{n}])^{2},

the squared determinant of the matrix with the columns u_{1}, …, u_{n}.

Proof. Let e_{1}, …, e_{n} be the standard orthonormal basis of ℝ^{n}. Let G be the operator such that G(e_{i}) = u_{i} for i = 1, …, n. Then u_{i} ⊗ u_{i} = G(e_{i} ⊗ e_{i})G^{∗} and from the definition of Section 1.2 we get

D(u_{1} ⊗ u_{1}, …, u_{n} ⊗ u_{n}) = \frac{∂^{n}}{∂t_{1} … ∂t_{n}} det(t_{1}u_{1} ⊗ u_{1} + … + t_{n}u_{n} ⊗ e_{n})

= \frac{∂^{n}}{∂t_{1} … ∂t_{n}} det(G(t_{1}e_{1} ⊗ e_{1} + … + t_{n}e_{n} ⊗ e_{n})G^{∗})

= det(GG^{∗})\frac{∂^{n}}{∂t_{1} … ∂t_{n}} det(t_{1}e_{1} ⊗ e_{1} + … + t_{n}e_{n} ⊗ e_{n}) = (det G)^{2}.

4. Proof of the Main Results

(4.1) Conditional expectations. In this subsection, we summarize some general results on measures and integration, which we exploit later. As a general source we use [3].

Suppose that we have k copies of the Euclidean space ℝ^{n}, each with the standard Gaussian probability measure μ. We will consider functions f : ℝ^{n} × … × ℝ^{n} ⟶ ℝ that are defined almost everywhere and integrable with respect to the measure ν_{k} = μ × … × μ. Let f(u_{1}, …, u_{k}) be such a function. Then for almost all (k − 1)-tuples (u_{1}, …, u_{k−1}) with u_{i} ∈ ℝ^{n}, the function f(u_{1}, …, u_{k−1}, ·) is integrable (Fubini’s Theorem) and we can define the conditional expectation g(u_{1}, …, u_{k−1}) = E_{k}(f) by letting

E_{k}(f)(u_{1}, …, u_{k−1}) = ∫_{ℝ^{n}} f(u_{1}, …, u_{k−1}, u_{k})ψ(u_{k}) du_{k}.

Fubini’s theorem implies that

E(f) = E_{1} … E_{k}(f),

where E is the expectation with respect to the product measure ν_{k}. Tonelli’s Theorem states that if f is ν_{k}-measurable and non-negative almost surely with respect to ν_{k} then f is ν_{k}-integrable, provided E_{1} … E_{k}(f) < +∞.

If f(u_{1}, …, u_{i}) is a function of i < k arguments, we may formally extend it to ℝ^{n} × … × ℝ^{n} (k times) by letting f(u_{1}, …, u_{k}) = f(u_{1}, …, u_{i}). The distribution of values of f(u_{1}, …, u_{i}) with respect to ν_{i} is the same as the distribution of values of f(u_{1}, …, u_{k}) with respect to ν_{k}. In particular, if f(u_{1}, …, u_{i}) is ν_{i}-integrable then f(u_{1}, …, u_{k}) is ν_{k}-integrable.

We note the following useful facts:

(4.1.1) The linear operator E_{k} is monotone, that is, if f(u_{1}, …, u_{k}) ≤ g(u_{1}, …, u_{k}) almost surely with respect to ν_{k}, then E_{k}(f) ≤ E_{k}(g) almost surely with respect to ν_{k−1}; (4.1.2) If f(u_{1}, …, u_{k}) is integrable and g(u_{1}, …, u_{i}), i < k is a function, then E_{k}(gf) = gE_{k}(f); (4.1.3) If f = a is a constant almost surely with respect to ν_{k}, then E_{k}(f) = a almost surely with respect to ν_{k−1}.

First, we prove the following technical lemma (a martingale inequality).

(4.2) Lemma. Let f_{k}(u_{1}, …, u_{k}), k = 1, …, n be integrable functions on ℝ^{n} × …×ℝ^{n} (k times), and let ν_{n} = μ×…×μ (n times). Suppose that for some numbers a and b we have

a ≤ E_{k}(f_{k}) and E_{k}(f_{k}^{2}) ≤ b almost surely with respect to ν_{k−1}

for k = 1, …, n.

Then for any δ > 0 we have

ν_{n}\{(u_{1}, …, u_{n}) : \frac{1}{n} ∑_{k=1}^{n} f_{k}(u_{1}, …, u_{k}) ≤ a − δ\} ≤ \frac{b}{nδ^{2}} .

Proof. Let g_{k} = E_{k}(f_{k}) and h_{k} = f_{k} − g_{k}. Since g_{k} does not depend on u_{k}, using (4.1.2) we have E_{k}(h_{k}^{2}) = E_{k}(f_{k}^{2}) − 2E_{k}(g_{k}f_{k}) + E_{k}(g_{k}^{2}) = E_{k}(f_{k}^{2}) − g_{k}^{2}. Summarizing, we may write

f_{k} = g_{k} + h_{k}, where E_{k}(h_{k}) = 0, g_{k} ≥ a, and E_{k}(h_{k}^{2}) ≤ b

almost surely with respect to ν_{k−1}. Let

H(u_{1}, …, u_{n}) = \frac{1}{n} ∑_{k=1}^{n} h_{k}(u_{1}, …, u_{k}).

Now for U = (u_{1}, …, u_{n}) we have

ν_{n}\{U : \frac{1}{n} ∑_{k=1}^{n} f_{k}(u_{1}, …, u_{k}) ≤ a − δ\}

= ν_{n}\{U : H(U) + \frac{1}{n} ∑_{k=1}^{n} g_{k}(u_{1}, …, u_{k−1}) ≤ a − δ\}

≤ ν_{n}\{U : H(U) ≤ −δ\} ≤ \frac{E(H^{2})}{δ^{2}} = \frac{1}{δn^{2}} ∑_{k=1}^{n} E(h_{k}^{2}) + \frac{2}{δn^{2}} ∑_{1≤i<j≤n} E(h_{i}h_{j})

(we used Chebyshev’s inequality).

We note that it is legitimate to pass to global expectations E here. Indeed, since h_{k}^{2} is non-negative and E_{k}h_{k}^{2} ≤ E_{k}f_{k}^{2} ≤ b it follows by (4.1.1), (4.1.3) and Tonelli’s Theorem that h_{k}^{2} is ν_{n}-integrable. Since |h_{i}h_{j}| ≤ (h_{i}^{2} + h_{j}^{2})/2, the products h_{i}h_{j} are also ν_{n}-integrable. Therefore, H^{2} is ν_{n}-integrable.

Since h_{k} does not depend on u_{k+1}, …, u_{n}, using (4.1.2) we have E(h_{k}^{2}) = E_{1} … E_{n}h_{k}^{2} = E_{1} … E_{k}h_{k}^{2} and since E_{k}h_{k}^{2} ≤ b almost surely on ν_{k−1}, by (4.1.1) and (4.1.3) we get that E(h_{k}^{2}) ≤ b for each k = 1, …, n. Furthermore, since h_{i} does not depend on u_{i+1}, …, u_{n} for j > i, using (4.1.2) and (4.1.3) we have

E(h_{i}h_{j}) = E_{1} … E_{n}(h_{i}h_{j}) = E_{1} … E_{j}(h_{i}h_{j}) = E_{i} … E_{j−1}h_{i}E_{j}(h_{j}) = 0.

The proof now follows.

Now we are ready to prove the main results of the paper.

(4.3) Proof of Theorem 2.2. The output α of Algorithm 2.1 is a function of vectors u_{1}, …, u_{n} ∈ ℝ^{n}, which are drawn at random from the standard Gaussian distribution μ in ℝ^{n}. We consider the distribution of α(u_{1}, …, u_{n}) with respect to the product measure ν_{n} = μ × … × μ (n times).

Using Lemma 3.4, we may write

α = α(u_{1}, …, u_{n}) = (det[T_{1}u_{1}, …, T_{n}u_{n}])^{2} = D(T_{1}u_{1} ⊗ T_{1}u_{1}, …, T_{n}u_{n} ⊗ T_{n}u_{n}).

For k = 0, …, n let

α_{k}(u_{1}, …, u_{k}) = D(T_{1}u_{1} ⊗ T_{1}u_{1}, …, T_{k}u_{k} ⊗ T_{k}u_{k}, Q_{k+1}, …, Q_{n}).

In particular, α_{0} = D(Q_{1}, …, Q_{n}) and α_{n} = α. Since the mixed discriminant of positive semidefinite matrices is non-negative (Section 1.2), we deduce that α_{k}(u_{1}, …, u_{k}) is non-negative. By (3.2.2) we have E_{k}(T_{k}u_{k} ⊗ T_{k}u_{k}) = T_{k}T_{k}^{∗} = Q_{k}. Since α_{k} is a polynomial in the coordinates of u_{1}, …, u_{k} and the mixed discriminant is linear in every argument (see Section 1.2) we can interchange D and E_{k}:

E_{k}(α_{k})

= D(T_{1}u_{1} ⊗ T_{1}u_{1}, …, T_{k−1}u_{k−1} ⊗ T_{k−1}u_{k−1}, E_{k}(T_{k}u_{k} ⊗ T_{k}u_{k}), Q_{k+1}, …, Q_{n})
= D(T_{1}u_{1} ⊗ T_{1}u_{1}, …, T_{k−1}u_{k−1} ⊗ T_{k−1}u_{k−1}, Q_{k}, Q_{k+1}, …, Q_{n}) = α_{k−1}.

Since α_{k} is non-negative and ν_{n}-measurable, applying theorems of Tonelli and Fubini, we have

(4.3.1)  E(α) = E_{1}E_{2} … E_{n}(α_{n}) = E_{1} … E_{k}(α_{k}) = α_{0} = D(Q_{1}, …, Q_{n}).

and (2.2.1) is proven. Since α(u_{1}, …, u_{n}) is non-negative, (2.2.2) follows by Chebyshev’s inequality.

If D(Q_{1}, …, Q_{n}) = 0 then by (2.2.1) and non-negativity of α it follows that α is identically 0 and (2.2.3) would follow. Therefore, without loss of generality, we may assume that D(Q_{1}, …, Q_{n}) > 0. Since α_{k}(u_{1}, …, u_{k}) is a non-negative polynomial, by (4.3.1) we conclude that α_{k}(u_{1}, …, u_{k}) > 0 almost surely with respect to ν_{k}. For k = 1, …, n let

g_{k}(u_{1}, …, u_{k}) = \frac{α_{k}(u_{1}, …, u_{k})}{α_{k−1}(u_{1}, …, u_{k−1})} = \frac{α_{k}}{E_{k}(α_{k})},

Hence we may write

\frac{α(u_{1}, …, u_{n})}{D(Q_{1}, …, Q_{n})} = ∏_{k=1}^{n} g_{k}(u_{1}, …, u_{k})

almost surely with respect to ν_{n}. Since the mixed discriminant is linear in each argument, g_{k}(u_{1}, …, u_{k}) is a quadratic form in u_{k} for any fixed u_{1}, …, u_{k−1}, such that α_{k−1}(u_{1}, …, u_{k−1}) > 0. Furthermore, since the mixed discriminant is nonnegative for positive semidefinite arguments, we conclude that g_{k}(u_{1}, …, u_{k}) is a positive semidefinite quadratic form in u_{k} for every such choice of u_{1}, …, u_{k−1}. Since E_{k}(g_{k}) = 1 almost surely with respect to ν_{k−1}, by Theorem 3.3 we conclude that

C_{0} ≤ E_{k}(ln g_{k}) ≤ 0 with C_{0} = \frac{4}{\sqrt{2π}} ∫_{0}^{+∞}(ln t)e^{−t^{2}/2} dt

and that |E_{k}(ln^{2} g_{k})| ≤ 8 almost surely with respect to ν_{k−1}. In particular, since ln^{2} g_{k} is non-negative almost surely with respect to ν_{k}, we deduce that ln^{2} g_{k} is ν_{n}- integrable, and since | ln g_{k}| ≤ 1 + ln^{2} g_{k} we deduce that ln g_{k} is ν_{n}-integrable.

Now for U = (u_{1}, …, u_{n}) we have

ν_{n}\{U : \frac{α(u_{1}, …, u_{n})}{D(Q_{1}, …, Q_{n})} ≤ (ϵc_{0})^{n}\} = ν_{n}\{U : \frac{1}{n} ln(\frac{α(u_{1}, …, u_{n})}{D(Q_{1}, …, Q_{n})}) ≤ C_{0} + ln ϵ\}

= ν_{n}\{(u_{1}, …, u_{n}) : \frac{1}{n} ∑_{k=1}^{n} ln g_{k}(u_{1}, …, u_{k}) ≤ C_{0} + ln ϵ\}.

To complete the proof of (2.2.3), we use Lemma 4.2 with f_{k} = ln g_{k}, a = C_{0}, b = 8 and δ = − ln ϵ.

(4.4) Proof of Theorem 2.4. For i = 1, …, n let us define the diagonal matrices Q_{i} =diag{a_{i1}, …, a_{in}}. Algorithm 2.3 with the input A is Algorithm 2.1 with the input Q_{1}, …, Q_{n}. By (1.2.1) per A = D(Q_{1}, …, Q_{n}), and the proof follows by Theorem 2.2.

5. Remarks

(5.1) The algorithms in the binary model. Suppose we want to operate in the standard binary (Turing) model of computation. That is, we allow arithmetic operations with integral (rational) numbers stored as bit strings (see [1]). In the probabilistic setting, we suppose also that we can generate a random bit, that is we can sample a random variable x ∈ {0, 1}, which assumes either value with probability 1/2. Algorithms 2.1 and 2.3 can be transformed into randomized polynomial time algorithms in this model, which approximate the mixed discriminant, resp. the permanent within a 2^{O(n)} factor. This transformation is relatively straightforward for Algorithm 2.3 and first we briefly sketch it here, omitting tedious details (although it may not be the most efficient binary version).

We are given a non-negative integer matrix A and we want to approximate per A. To compute the matrix B from Algorithm 2.3 we need to compute square roots \sqrt{a_{ij}} and to sample variables u_{ij} from the standard Gaussian distribution in ℝ. It is known that the square root of a positive rational number can be computed within any given error ϵ > 0 in time that is polynomial in ln ϵ^{−1}. We observe that if we choose u_{ij} : 1 ≤ i, j ≤ n independently from the standard Gaussian distribution in ℝ, then we will have |u_{ij}| ≤ n for all i, j with the probability that goes to 1 exponentially fast as n ⟶ +∞. So, we compute \sqrt{a_{ij}} for a_{ij} ≥ 1 with such a precision ϵ that for any choice of u_{ij} : |u_{ij}| ≤ n the value of the output α = (det B)^{2} gets computed with an error at most 10^{−n} (we recall that to compute the determinant we need arithmetic operations only). One can show that the bit size of ϵ can be bounded by a polynomial in the input size. We note that if A is a 0-1 matrix, we don’t need square roots.

The next step is to approximate sampling from the standard Gaussian distribution in ℝ. Let us choose a variation of the “polar method” (see Section 3.4.1 (C) of [17]). It can be shown that if x and y are independent random variables, uniformly distributed on [0, 1], then u = sin(2πy)\sqrt{−2 ln x} has the standard Gaussian distribution in ℝ. Let (x_{ij}, y_{ij}) : 1 ≤ i, j ≤ n be the coordinates in the 2n^{2}-dimensional unit cube [0, 1]^{n^{2}} × [0, 1]^{n^{2}} = [0, 1]^{2n^{2}} and u_{ij} be the coordinates in ℝ^{n^{2}}. This allows us to construct a map

Φ : [0, 1]^{n^{2}} × [0, 1]^{n^{2}} ⟶ ℝ^{n^{2}}, where Φ_{ij}(x, y) = sin(2πy_{ij})\sqrt{−2 ln x_{ij}},

such that the push-forward measure of the standard Lebesgue measure λ on [0, 1]^{2n^{2}} is the standard Gaussian measure μ in ℝ^{n^{2}}. The output α(u_{ij}) is a function on ℝ^{n^{2}}. The map Φ is singular on the part of the boundary where x_{ij} = 0 or x_{ij} = 1, so we find a parallelepiped Π ⊂ [0, 1]^{2n^{2}}, such that λ(Π) ≥ 0.99 for all sufficiently large n and 1−1/n^{3} ≥ x_{ij} ≥ 1/n^{3} for any point in Π. The map Φ restricted to Π is Lipschitz, so we are able to find a rational δ > 0 such that if ∥z_{1} − z_{2}∥ ≤ δ for z_{1}, z_{2} ∈ Π then |α(Φ(z_{1})) − α(Φ(z_{2}))| ≤ 10^{−n}. It can be shown that the size of δ can be bounded by a polynomial in the input size. Next, we have to approximate sampling from the Lebesgue measure on Π by the binary sampling. We do it by sampling each coordinate x_{ij} and y_{ij} independently. To sample a coordinate x we sample N random bits b_{1}, …, b_{N} and let

x = ∑_{k=1}^{N} 2^{−k}b_{k}.

The 2^{2n^{2}N} points that can be obtained in such a way form a lattice grid in [0, 1]^{2n^{2}}. We choose N so large that this grid forms a δ-net in [0, 1]^{2n^{2}}. One can choose an N that is bounded by a polynomial in n and ln δ. Now we can approximate sampling from the standard Gaussian measure in ℝ^{n^{2}}: we sample a point z ∈ [0, 1]^{2n^{2}} as above, accept it if z ∈ Π and compute u = Φ(a) approximately with a sufficiently high precision, so that the corresponding value of α(u) = (det B)^{2} gets computed with an error at most 2 · 10^{−n} (taking into account the error from approximate computation of \sqrt{a_{ij}}). Again, this can be done in polynomial time.

Finally, the output α ≤ 3 · 10^{−n} is rounded to 0.

It is interesting to note, that even to approximate the permanent of a 0-1 matrix (a problem, which sounds purely combinatorial) we seem to have to deal with approximate computation of such “non-combinatorial” functions as sin and ln.

Algorithm 2.3 is modified similarly, except that on the first step we perturb the matrices Q_{i} ⟼ Q_{i} + ϵI, where I is the identity matrix and ϵ > 0 is such that the value of D(Q_{1}, …, Q_{n}) changes by not more than 10^{−n} (the bit size of ϵ can be bounded by a polynomial in the size of the input). Now Q_{i} are strictly positive definite and this makes computation of the decomposition Q_{i} = T_{i}T_{i}^{∗} stable, so that we can compute T_{i} with the desired precision.

(5.2) An application of the mixed discriminant to counting. Suppose we are given a rectangular n × m matrix A with the columns u_{1}, …, u_{m}, which we interpret as vectors from ℝ^{n}. Suppose that for any subset I ⊂ {1, …, m}, I = {i_{1}, …, i_{n}}, the determinant of the submatrix A_{I} = [u_{i_{1}}, …, u_{i_{n}}] is either 0, −1 or 1. Such an A represents a unimodular matroid on the set {1, …, m} and a subset I with det A_{I} ≠ 0 are called a base of the matroid (see [22]).

Suppose now, that the columns of A are colored into n different colors, which induces a partition {1, …, m} = J_{1} ∪ … ∪ J_{n}. We are interested in the number of bases that have precisely one index of each color. Let us define the positive semidefinite matrices Q_{1}, …, Q_{n} as follows:

Q_{k} = ∑_{i∈J_{k}} u_{i} ⊗ u_{i}, k = 1, …, n.

Then the number of such bases can be expressed as D(Q_{1}, …, Q_{n}). Indeed, using the linearity of the mixed discriminant (Section 1.2) and Lemma 3.4, we have

D(Q_{1}, …, Q_{n}) = ∑_{I={i_{1},... ,i_{n}}} D(u_{i_{1}} ⊗ u_{i_{1}}, …, u_{i_{n}} ⊗ u_{i_{n}})

= ∑_{I={i_{1},... ,i_{n}}}(det[u_{i_{1}}, …, u_{i_{n}}])^{2},

where the sums are taken over all n-subsets I, having precisely one element of each color and the proof follows.

(5.2.1) Example. Trees in a graph. Suppose we have a connected graph G with n vertices and m edges. Suppose further, that the edges of G are colored into n − 1 different colors. We are interested in spanning trees T in G such that all edges of T have different colors. Let us number the vertices of G by 1, …, n and the edges of G by 1, …, m. Let us make G an oriented graph by orienting its edges arbitrarily. We consider the truncated incidence matrix (with the last row removed) A = (a_{ij}) for 1 ≤ i ≤ n − 1 and 1 ≤ j ≤ m as an (n − 1) × m matrix such that

a_{ij} = \{ \begin{matrix}1 & if i is the beginning of j \\ −1 & if i is the end of j \\ 0 & otherwise. \end{matrix}

The spanning trees of G are in one-to-one correspondence with non-degenerate (n − 1) × (n − 1) submatrices of A and the determinant of a such a submatrix is either 1 or −1 (see, for example, Chapter 4 of [7]). Hence counting colored trees reduces to computation of the mixed discriminant of some positive semidefinite matrices, computed from the incidence matrix of the graph.

Acknowledgment

This research was partially supported by Alfred P. Sloan Research Fellowship, by NSF grant DMS 9501129 and by the Mathematical Sciences Research Institute, Berkeley, CA through NSF grant DMS 9022140.

References

[1] A. Aho, J. Hopcroft and J. Ullman, The Design and Analysis of Computer Algorithms, Addison-Wesley, Reading, 1974.

[2] A.D. Aleksandrov, On the theory of mixed volumes of convex bodies, IV, Mixed discriminants and mixed volumes (in Russian), Math. Sb., N.S., 3(1938), 227-251.

[3] R.G. Bartle, The Elements of Integration, John Wiley & Sons, New York, 1966.

[4] A.I. Barvinok, Two algorithmic results for the Traveling Salesman Problem, Mathematics of Operations Research, 21(1996), 65-84.

[5] A. I. Barvinok, Computing mixed discriminants, mixed volumes and permanents, Discrete & Computational Geometry, to appear.

[6] G.P. Egorychev, The solution of van der Waerden’s problem for permanents. Advances in Mathematics, 42 (1981), 299-305.

[7] V.A. Emelichev, M.M. Kovalev and M.K. Kravtsov, Polytopes, Graphs and Optimization, Cambridge University Press, New York, 1984.

[8] Faddeeva V.N., Computational Methods of Linear Algebra, Dover, 1975.

[9] A. Frieze and M. Jerrum, An analysis of a Monte Carlo algorithm for estimating the permanent, Combinatorica, 15 (1995), 67-83.

[10] C.D. Godsil and I. Gutman, On the matching polynomial of a graph, in.: Algebraic Methods in Graph Theory, Vol. I, II, (Szeged, 1978), North-Holland, Amsterdam - New York, 1981, 241-249.

[11] D.Yu. Grigoriev and M. Karpinski, The matching problem for bipartite graphs with polynomially bounded permanents is in NC, Proc. Twenty-eight Annual IEEE Sympos. Foundations of Computer Sci., IEEE Computer Soc. Press, Washington D.C., 1987, 162-172.

[12] M. Jerrum and A. Sinclair, Approximating the permanent, SIAM Journal on Computing, 18(1989), 1149-1178.

[13] M. Jerrum and A. Sinclair, The Markov chain Monte Carlo method: an approach to approximate counting and integration, in: Approximation Algorithms for NP-hard Problems, D.S. Hochbaum, ed., PWS Publishing, Boston, MA, 1996, 482-520.

[14] M.R. Jerrum, L.G. Valiant and V.V. Vazirani, Random generation of combinatorial structures from a uniform distribution, Theoretical Computer Science, 43 (1986), 169-188.

[15] M. Jerrum and U. Vazirani, A mildly exponential approximation algorithm for the permanent, Algorithmica, 16 (1996), 392–401.

[16] N. Karmarkar, R. Karp, R. Lipton, L. Lovász and M. Luby, A Monte Carlo algorithm for estimating the permanent, SIAM Journal on Computing, 22(1993), 284-293.

[17] D.E. Knuth, The Art of Computer Programming. Vol. 2: Seminumerical Algorithms, second edition, Addison-Wesley, MA 1981.

[18] K. Leichtweiß, Convexity and Differential Geometry, in: Handbook of Convex Geometry, ed. P.M. Gruber and J.M. Wills, vol. B, Chapter 4.1, North-Holland, 1993, 1045-1080.

[19] L. Lovász and M.D. Plummer, Matching Theory, North - Holland, Amsterdam - New York and Akadémiai Kiadó, Budapest, 1986.

[20] H. Minc, Permanents, Addison-Wesley, Reading, MA, 1978.

[21] L.E. Rasmussen, Approximating the permanent: a simple approach, Random Structures & Algorithms, 5(1994), 349-361.

[22] N. White, Unimodular matroids, in: Combinatorial Geometries, Cambridge Univ. Press, Cambridge, 1987, 40-52.

Alexander Barvinok, Department of Mathematics, University of Michigan, Ann Arbor, MI 48109-1109, USA

E-mail address: barvinok@math.lsa.umich.edu