#entropy zero g/t | Explore Tumblr posts and blogs

theresawelchy · 6 years ago

Text

tmp

Lecturer: Aram Harrow

Scribes: Sinho Chewi, William Moses, Tasha Schoenstein, Ary Swaminathan

November 9, 2018

Outline

Sampling from thermal states was one of the first and (initially) most important uses of computers. In this blog post, we will discuss both classical and quantum Gibbs distributions, also known as thermal equilibrium states. We will then discuss Markov chains that have Gibbs distributions as stationary distributions. This leads into a discussion of the equivalence of mixing in time (i.e. the Markov chain quickly equilibrates over time) and mixing in space (i.e. sites that are far apart have small correlation). For the classical case, this equivalence is known. After discussing what is known classically, we will discuss difficulties that arise in the quantum case, including (approximate) Quantum Markov states and the equivalence of mixing in the quantum case.

Gibbs distributions

We have already learned about phase transitions in a previous blog post, but they are important, so we will review them again. The Gibbs or thermal distribution is defined as follows: Suppose that we have an energy function E:{0,1}n→RE:{0,1}n→R

E : {\{0,1\}}^n \to {\mathbb R}, which takes nn

n-bit strings to real numbers. Usually, E=∑mi=1EiE=∑i=1mEi

E = \sum_{i=1}^m E_i, where each EiEi

E_i term depends only on a few bits. For example, the energy might be the number of unsatisfied clauses in a 3-SAT formula, or it may arise from the Ising model. The Gibbs distribution is

p(x)=exp{−E(x)/T}Zp(x)=exp⁡{−E(x)/T}Z

p(x) = \frac{\exp\{-E(x)/T\}}{Z}

where the normalization factor in the denominator, also called the partition function, is Z=∑x∈{0,1}nexp{−E(x)/T}Z=∑x∈{0,1}nexp⁡{−E(x)/T}

Z = \sum_{x \in {\{0,1\}}^n} \exp\{-E(x)/T\}. Another, perhaps more operational, way to define the Gibbs distribution is:

<span id="MathJax-Element-8-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" style="font-size: 109%; text-align: center; position: relative;" data-mathml="p=argmaxq∈P({0,1}n)⁡H(q) subject to the constraint ⟨q,E⟩=E¯.” role=”presentation”>p=argmaxq∈P({0,1}n)H(q) subject to the constraint ⟨q,E⟩=¯E.p=argmaxq∈P({0,1}n)⁡H(q) subject to the constraint ⟨q,E⟩=E¯.

p = \operatorname*{arg\,max}_{q \in {\mathcal{P}}({\{0,1\}}^n)} H(q)~\text{subject to the constraint}~\braket{q,E} = \bar{E}.

In this expression, P({0,1}n)P({0,1}n)

{\mathcal{P}}({\{0,1\}}^n) is the set of probability distributions on {0,1}n{0,1}n

{\{0,1\}}^n, HH

H is the Shannon entropy, and ¯EE¯

\bar E is a constant representing the average energy. We are thinking of probability distributions and EE

E as vectors of size 2n2n

2^n. It turns out that if we solve this optimization problem, then the Gibbs distribution is the unique solution.

Uses of Gibbs distributions

Why is it useful to work with Gibbs distributions?

Gibbs distributions arise naturally in statistical physics systems, such as constraint satisfaction problems (CSPs), the Ising model, and spin glasses. One approach to deal with Gibbs distributions is through belief propagation (BP), which yields exact inference on tree graphical models and sometimes phase transition predictions on loopy graphs. Instead, we will focus on a different approach, namely, sampling from the Gibbs distribution.

If we want to minimize EE

E (say, to find a 3-SAT solution), we can use simulated annealing. The idea of annealing is that we want to produce a crystal; a crystal is the lowest energy configuration of molecules. If we heat up the substance to a liquid and then cool it quickly, we will not get a nice crystal, because little bits of the material will point in different directions. In order to form a crystal, we need to cool the system slowly.

In computer science terms, we take a sample from a high temperature because sampling is generally easier at a higher temperature than at a lower temperature. We then use that sample as the starting point for an equilibration process at a slightly lower temperature, and repeat this procedure. If we reach zero temperature, then we are sampling from the minimizers of EE

E. In practice, the system will usually stop mixing before we get to zero temperature, but this is a good heuristic. You can think of this process as gradient descent, with some additional randomness.

Gibbs distributions are used to simulate physical systems.

Gibbs distributions are used in Bayesian inference due to the Hammersley-Clifford theorem, which will be discussed next.

Gibbs distributions are also connected to multiplicative weights for linear programming (not discussed in this blog post).

Bayesian inference & the Hammersley-Clifford theorem

In order to present the Hammersley-Clifford theorem, we must first discuss Markov networks. For this part, we will generalize our setup to a finite alphabet ΣΣ

\Sigma, so the energy function is now a function Σn→RΣn→R

\Sigma^n \to \mathbb R.

Markov chains

First, let us recall the idea of a Markov chain with variables X1X1

X_1, X2X2

X_2, X3X3

X_3.

The random variables X1X1

X_1, X2X2

X_2, X3X3

X_3 form a Markov chain if their joint distribution can be written in a factored way: p(x1,x2,x3)=p1,2(x1,x2)p3∣2(x3∣x2)p(x1,x2,x3)=p1,2(x1,x2)p3∣2(x3∣x2)

p(x_1,x_2,x_3) = p_{1,2}(x_1,x_2)p_{3 \mid 2}(x_3 \mid x_2). For example, imagine that X1X1

X_1, X2X2

X_2, X3X3

X_3 represent the weather on Monday, Tuesday, and Wednesday respectively. These random variables form a Markov chain if, conditioned on the weather on Tuesday, we have all of the information we need to forecast the weather on Wednesday. Another way to say this is that conditioned on the weather on Tuesday, then the weather on Monday and the weather on Wednesday are conditionally independent. Note that the weather on Monday and the weather on Wednesday are not independent; there can be correlations, but these correlations are mediated through the weather on Tuesday. It is important to note that the definition of a Markov chain is symmetric with respect to going forwards or backwards in time, so we can also write the conditional independence condition as p(x1,x2,x3)=p2,3(x2,x3)p1∣2(x1∣x2)p(x1,x2,x3)=p2,3(x2,x3)p1∣2(x1∣x2)

p(x_1,x_2,x_3) = p_{2,3}(x_2,x_3) p_{1 \mid 2}(x_1 \mid x_2).

The conditional independence condition can also be written as p1,3∣2(x1,x3∣x2)=p1∣2(x1∣x2)p3∣2(x3∣x2).p1,3∣2(x1,x3∣x2)=p1∣2(x1∣x2)p3∣2(x3∣x2).

p_{1,3 \mid 2}(x_1, x_3 \mid x_2) = p_{1 \mid 2}(x_1 \mid x_2) p_{3 \mid 2}(x_3 \mid x_2). Recall that for two random variables X1X1

X_1 and X2X2

X_2 with joint distribution pp

p, they are independent, i.e., p(x1,x2)=p1(x1)p2(x2)p(x1,x2)=p1(x1)p2(x2)

p(x_1,x_2) = p_1(x_1) p_2(x_2), if and only if I(X1;X2)=0I(X1;X2)=0

I(X_1; X_2) = 0, where II

I here denotes the mutual information. Similarly, conditional independence is equivalent to the conditional mutual information I(X1;X3∣X2)I(X1;X3∣X2)

I(X_1; X_3 \mid X_2) equaling zero. This quantity is defined as I(X1;X3∣X2)=H(X1∣X2)+H(X3∣X2)−H(X1,X3∣X2)I(X1;X3∣X2)=H(X1∣X2)+H(X3∣X2)−H(X1,X3∣X2)

I(X_1;X_3 \mid X_2) = H(X_1 \mid X_2) + H(X_3 \mid X_2) – H(X_1, X_3 \mid X_2).

Keep in mind that conditional independence is characterized in two equivalent ways: via an algebraic condition on the distributions, and via mutual information.

Markov networks

A Markov network is like a Markov chain, but with more random variables and a more interesting structure. Imagine that we have a graph, where each node is associated with a random variable and the edges encode possible correlations. A Markov network has the property that if we take any disjoint collection of nodes AA

A, BB

B, and CC

C such that AA

A and CC

C are fully separated by BB

B (that is, any path from AA

A to CC

C must go through BB

B, or alternatively, removing BB

B leaves AA

A and CC

C disconnected), then I(XA;XC∣XB)=0I(XA;XC∣XB)=0

I(X_A; X_C \mid X_B) = 0. The notation XAXA

X_A here means the collection of random variables associated with the nodes in AA

For example:

Here, if A={1,5,6}A={1,5,6}

A=\{1,5,6\}, B={2,7}B={2,7}

B=\{2,7\}, and C={3,4}C={3,4}

C=\{3,4\}, then BB

B separates AA

A and CC

A Markov network is also called a graphical model or a Markov random field; and yet another name for them is Gibbs distribution, which is the content of the following theorem:

Theorem 1 (Hammersley-Clifford Theorem): Let pp

p be a strictly positive distribution on ΣnΣn

\Sigma^n. Then, pp

p can be represented as a Markov network with respect to a graph GG

G if and only if pp

p can be expressed as a Gibbs distribution p(x)∝exp{−∑C∈C(G)EC(xC)}p(x)∝exp⁡{−∑C∈C(G)EC(xC)}

p(x) \propto \exp\{-\sum_{C \in {\mathcal{C}}(G)} E_C(x_C)\}, where C(G)C(G)

{\mathcal{C}}(G) is the set of cliques (fully connected subsets) of GG

This theorem says that Markov networks are the same as Gibbs states, with the same notion of locality.

The Hammersley-Clifford theorem implies an area law for mutual information; we will explain what this is and sketch why this is true. Divide a system into two disjoint pieces AA

A and BB

B. We want to know about the mutual information between AA

A and BB

B, I(A;B)I(A;B)

I(A;B). The Hammersley-Clifford theorem gives us a bound which depends only on the size of the boundary ∂∂

\partial between these sets. For simplicity, assume ∂⊆B∂⊆B

\partial \subseteq B. Also, assume that the interactions have bounded range; then, the Hammersley-Clifford theorem tells us that I(A;B∣∂)=0I(A;B∣∂)=0

I(A; B \mid \partial) = 0.

Now, we will use the fact I(A;B∣∂)=I(A;B,∂)−I(A;∂)I(A;B∣∂)=I(A;B,∂)−I(A;∂)

I(A; B \mid \partial) = I(A; B,\partial) – I(A; \partial). We can see this by writing out the expressions, but the intuition is that the term on the left asks about how much AA

A knows about BB

B, having already known about ∂∂

\partial. This equals how much AA

A knows about BB

B and ∂∂

\partial combined, minus how much AA

A knows about ∂∂

\partial alone. In this case, since we said ∂⊆B∂⊆B

\partial \subseteq B, then I(A;B)I(A;B)

I(A; B) is the same as I(A;B,∂)I(A;B,∂)

I(A; B, \partial). In general, however, we have an upper bound:

<span id="MathJax-Element-88-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" style="font-size: 109%; text-align: center; position: relative;" data-mathml="I(A;B)≤I(A;B,∂)=I(A;∂)+I(A;B∣∂)≤H(∂)≤|∂|log⁡|Σ|” role=”presentation”>I(A;B)≤I(A;B,∂)=I(A;∂)+I(A;B∣∂) ≤H(∂)≤|∂|log|Σ|I(A;B)≤I(A;B,∂)=I(A;∂)+I(A;B∣∂)≤H(∂)≤|∂|log⁡|Σ|

I(A;B) \le I(A; B, \partial) = I(A; \partial) + \cancel{I(A;B \mid \partial)} \le H(\partial) \le |\partial| \log |\Sigma|

In this calculation, we have used I(A;∂)=H(∂)−H(∂∣A)I(A;∂)=H(∂)−H(∂∣A)

I(A; \partial) = H(\partial) – H(\partial \mid A) (the information between AA

A and ∂∂

\partial is the amount by which the entropy of ∂∂

\partial gets reduced once we know AA

A) and H(∂∣A)≥0H(∂∣A)≥0

H(\partial \mid A) \ge 0 (which is true classically).

Since the mutual information only scales with the surface area of the boundary and not with the area of the two regions AA

A and BB

B, this is known as an area law [1].

Relationship to Bayesian inference

In Bayesian inference, we have a model for a system which can be very complicated. The model represents our assumptions on how parts of the system are causally related to the rest of the system. We have some observations, and we want to sample from a distribution conditionally on the fixed observations. Sampling from a conditional distribution is not the same as sampling from the original distribution, but we can still formally represent the conditional distribution as a Markov network. Therefore, sampling from Markov networks is a broadly useful task.

As an example of a complicated Bayesian model, consider a hierarchical Bayesian model [2]. Bayesian statistics requires choosing a prior distribution, and when there is a natural parameterized family of priors that a statistician can use, it may make sense to introduce a distribution over the priors; this is known as introducing a hyperparameter, and inference in the resulting hierarchical model (including computation of the posterior distribution) is frequently intractable. However, it is still desirable to work with these models because they are often more accurate than models in which the prior is handpicked by a statistician.

Sampling from Gibbs distributions

The task of sampling from an arbitrary Gibbs distribution is MA-complete [3], and it is not hard to see that at low enough temperatures this problem is at least NP-hard. So, how do we sample from these distributions?

This section will discuss Monte Carlo Markov chain (MCMC) methods, namely the Metropolis-Hastings algorithm and Glauber dynamics. Readers familiar with these methods may wish to skip to the discussion of mixing in time. For readers who wish to build more intuition about Markov chains before proceeding, see the Appendix, where the simple example of the random walk on a cycle is treated in detail.

Monte Carlo Markov chain (MCMC) methods

The general approach is to use a Markov chain. Let Ω=ΣnΩ=Σn

\Omega=\Sigma^n be the possible states of the system. Effectively, a Markov chain is a way of doing a random walk over ΩΩ

\Omega.

The transition probabilities of the Markov chain are1 P{X(t+1)=y∣X(t)=x}=Ty,x.P{X(t+1)=y∣X(t)=x}=Ty,x.

{\mathbb P}\{X(t+1) = y \mid X(t) = x\} = T_{y,x}. Here, TT

T is the transition probability matrix. The column at index xx

x of TT

T is the probability distribution of the next state of the Markov chain, if the current state is xx

x. The row at index yy

y is a row of probability values which give the probabilities of jumping into state yy

y from every other state. It has the properties that its entries are non-negative and for every xx

x, ∑y∈ΩTy,x=1∑y∈ΩTy,x=1

\sum_{y \in \Omega} T_{y,x} = 1. These properties say that TT

T is a (column) stochastic matrix.

Suppose we start at a state x(0)x(0)

x(0); or, more generally, we will start with a distribution pp

p over ΩΩ

\Omega. If we move according to the chain once, the distribution will be TpTp

Tp. If we move agian, the distribution will be T2pT2p

T^2 p. In general, after tt

t movements, the distribution is TtpTtp

T^t p. So, we can express the dynamics of the chain as matrix-vector multiplication.

It is worth mentioning that if we are simulating the chain on a computer and we are manipulating nn

n-bit numbers, then these probability vectors are of size 2n2n

2^n so it becomes impractical to store the entire probability distributions.

The justification for our algorithms is the following theorem.

Theorem 2 (Perron-Frobenius Theorem): If TT

T is a stochastic aperiodic matrix, then one of the eigenvalues is 11

1, and all other eigenvalues have magnitude strictly less than 11

1. There is a unique probability distribution ππ

\pi such that Tπ=πTπ=π

T\pi = \pi.

The theorem implies that TtpTtp

T^t p will converge to the stationary distribution ππ

\pi as t→∞t→∞

t\to\infty. So, if we want to sample from a distribution, this provides a method of doing so: cook up a Markov chain that equilibrates to the desired distribution, and then run the Markov chain until convergence. A priori, it is not obvious how we can design the Markov chain. At first, our problem was to sample from a probability distribution (a vector), and now we have changed the problem to designing an entire matrix, which does not appear to make our task easier.

Now, the question becomes: how does one come up with Markov chains that give you the desired stationary distribution?

Metropolis-Hastings algorithm

The first algorithm we will introduce is the Metropolis-Hastings algorithm. One more desirable feature of a Markov chain is that it satisfies detailed balance, which says πxTy,x=πyTx,yπxTy,x=πyTx,y

\pi_x T_{y,x} = \pi_y T_{x,y} for all xx

x and yy

y. This condition says that if we pick a point with probability according to the stationary distribution and transition, the probability of picking xx

x and then moving to yy

y should be the same as picking yy

y and then moving to xx

For a Markov chain in equilibrium, the total amount of probability flowing out of xx

x must equal the total amount of probability flowing into xx

x. For example, the United States might export products to Europe and import from China. Detailed balance says that the flow along each edge must balance, which is a more demanding condition. In the example with country trade deficits, we are requiring that all bilateral trade deficits must be zero.

Mathematically, detailed balance implies that TT

T can be transformed, via similarity transformations, into a symmetric matrix. The Metropolis-Hastings algorithm says that we should choose TT

T with the property Tx,yTy,x=πxπy.Tx,yTy,x=πxπy.

\frac{T_{x,y}}{T_{y,x}} = \frac{\pi_x}{\pi_y}. Suppose that we have an underlying graph on our state space, and suppose that we are at a state xx

x. The algorithm chooses a random neighbor, say yy

y, and then accepts or rejects this move with some probability. If the move is accepted, then we move to yy

y and continue the algorithm from there. Otherwise, if the move is rejected, then we stay at xx

x. We are free to choose any underlying graph (as long as it is connected and has a self-loop), and then we will tune the acceptance probability so that detailed balance holds.

Look at the trial move x→yx→y

x\to y. One way we can accomplish detailed balance is by looking at the ratio πy/πxπy/πx

\pi_y/\pi_x. If πy/πx≥1πy/πx≥1

\pi_y/\pi_x \ge 1, then always accept the move. If πy/πx<1πy/πx<1

% <![CDATA[ \pi_y/\pi_x , then accept the move with probability πx/πyπx/πy

\pi_x/\pi_y.

To get an idea for how the algorithm works, suppose that our underlying graph is dd

d-regular. Then, for neighbors xx

x and yy

<span id="MathJax-Element-150-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" style="font-size: 109%; text-align: center; position: relative;" data-mathml="Ty,x=min{1,πyπx}1d,Tx,y=min{1,πxπy}1d” role=”presentation”>Ty,x=min{1,πyπx}1d,Tx,y=min{1,πxπy}1dTy,x=min{1,πyπx}1d,Tx,y=min{1,πxπy}1d

Claim: Ty,xπx=1dmin{πx,πy},Ty,xπx=1dmin{πx,πy},

T_{y,x} \pi_x = \frac{1}{d} \min\{\pi_x,\pi_y\}, which is manifestly symmetric in xx

x and yy

y; thus, we have reversibility. This is the basic idea of the Metropolis-Hastings algorithm.

How does it work for a Gibbs distribution πx=exp{−E(x)/T}/Zπx=exp⁡{−E(x)/T}/Z

\pi_x = \exp\{-E(x)/T\}/Z, where the energy function might, for example, count the number of violated clauses in a 3-SAT formula? In this case, we might be a little worried. The numerator of πxπx

\pi_x is pretty easy to compute (we can count how many violated constraints there are), but the denominator is hard to compute. In general, it is #P-hard to compute the denominator, because as TT

T drops to 00

0, the partition function in this case approaches the number of 3-SAT solutions. So, how do we calculate the ratios πy/πxπy/πx

\pi_y/\pi_x that the algorithm requires? We’re able to do this because the ratio does not depend on ZZ

πyπx=expE(x)−E(y)T.πyπx=exp⁡E(x)−E(y)T.

\frac{\pi_y}{\pi_x} = \exp \frac{E(x)-E(y)}{T}.

Suppose that the energy is a sum of local terms, and the underlying graph corresponds to modifying one site at at a time. What this means is that the graph is Ω={0,1}nΩ={0,1}n

\Omega = {\{0,1\}}^n and the edges in the graph correspond to flipping exactly one bit. In this case, it becomes very easy to evaluate the computations needed for the algorithm; in fact, we can even do them in parallel.

How do we choose the underlying graph? The key idea is that we do not want the majority of our moves to be rejected. A good example to keep in mind is the Ising model, where the configurations are x∈{0,1}nx∈{0,1}n

x \in {\{0,1\}}^n and the energy is E(x)=−∑ni,j=1Ji,jxixjE(x)=−∑i,j=1nJi,jxixj

E(x) = -\sum_{i,j=1}^n J_{i,j} x_i x_j. If Ji,j≥0Ji,j≥0

J_{i,j} \ge 0 for all ii

i, jj

j, then we say that the model is ferromagnetic (we obtain lower energy by making the sites agree with each other). Of course, an antiferromagnetic model is just the opposite of this.

Assume that the bits are laid out in a square and Ji,j=JJi,j=J

J_{i,j} = J if ii

i and jj

j are neighbors on the square, and Ji,j=0Ji,j=0

J_{i,j} = 0 if they are not. As we vary the quantity J/TJ/T

J/T, we observe a phase transition. If J/TJ/T

J/T is small, then the coupling between the random variables is weak and the different parts of the system are almost independent; we call this the disordered phase. If J/TJ/T

J/T is large, then the spins want to align in the same direction and the Gibbs distribution will look almost like the following: with probability 1/21/2

1/2, all spins are ++

+, and with probability 1/21/2

1/2, all spins are −−

-; we call this the ordered phase.

In the disordered phase, when the spins do not need to align so closely, the Metropolis-Hastings algorithm will work well. In the ordered phase, the algorithm is doomed. Indeed, suppose that most of the spins are ++

+. As time proceeds, any −−

-s will switch to ++

+. There may be islands of −−

– spins initially, but it will be energetically favorable for these islands to shrink over time. Therefore, there will be an exponentially small chance for the system to switch to a configuration with mostly −−

-’s, and thus the chain takes exponentially long to mix. Here, people are interested in understanding the autocorrelation time, because the goal is to run the chain for some time, get one sample, run the chain for some more time, get another sample, etc.

Glauber dynamics

This next method (Glauber dynamics) is essentially the same as Metropolis-Hastings, but this is not immediately obvious. We are at a state x=(x1,…,xn)∈Σnx=(x1,…,xn)∈Σn

x = (x_1,\dotsc,x_n) \in \Sigma^n. (For the Metropolis-Hastings algorithm, we could be walking on a state space without a product structure. However, Glauber dynamics requires a product structure.) Then, we update xx

x to (x1,…,xi−1,x′i,xi+1,…,xn)(x1,…,xi−1,xi′,xi+1,…,xn)

(x_1,\dotsc,x_{i-1},x_i’,x_{i+1},\dotsc,x_n) with chance πi∣−i(x′i∣x1,…,xi−1,xi+1,…,xn)πi∣−i(xi′∣x1,…,xi−1,xi+1,…,xn)

\pi_{i\mid -i}(x_i’ \mid x_1,\dotsc,x_{i-1},x_{i+1},\dotsc,x_n). In other words, we hold all other bits fixed, and conditioned on those other bits, we resample the ii

ith bit. Like Metropolis-Hastings, ππ

\pi is stationary for this chain.

It is not obvious that these conditional distributions can be computed efficiently, but it is possible since normalizing the conditional distribution only requires summing over the possible configurations for a single random variable. On a Markov network, the conditional probability is πi∣N(i)(x′i∣xN(i))πi∣N(i)(xi′∣xN(i))

\pi_{i \mid N(i)}(x_i’ \mid x_{N(i)}), where N(i)N(i)

N(i) denotes the set of neighbors of ii

i. This makes the computation a constant-sized calculation (i.e., does not depend on the size of the system).

For example, in the Ising model, suppose we are at state x∈{±1}nx∈{±1}n

x \in {\{\pm 1\}}^n. In Glauber dynamics, we pick a vertex i∈[n]i∈[n]

i \in [n] u.a.r. and update it to ++

+ with probability <span id="MathJax-Element-195-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" data-mathml="pi∣N(i)(+∣xN(i))=exp⁡(T−1∑j∈N(i)xj)exp⁡(−T−1∑j∈N(i)xj)+exp⁡(T−1∑j∈N(i)xj).” role=”presentation” style=”font-size: 109%; position: relative;”>pi∣N(i)(+∣xN(i))=exp(T−1∑j∈N(i)xj)exp(−T−1∑j∈N(i)xj)+exp(T−1∑j∈N(i)xj).pi∣N(i)(+∣xN(i))=exp⁡(T−1∑j∈N(i)xj)exp⁡(−T−1∑j∈N(i)xj)+exp⁡(T−1∑j∈N(i)xj).

p_{i \mid N(i)}(+ \mid x_{N(i)}) = \frac{\exp(T^{-1}\sum_{j\in N(i)} x_j)}{\exp(-T^{-1} \sum_{j\in N(i)} x_j) + \exp(T^{-1} \sum_{j\in N(i)} x_j)}.

Mixing in time

Mixing in time means that the dynamics will equilibrate rapidly. It turns out that this is equivalent to mixing in space, which means that ππ

\pi itself has decaying correlations. For example, the Ising model at low temperature has a lot of long-range correlations, but at high temperature it does not. For the high temperature regime, we can prove that mixing in time occurs. We will prove this for the ferromagnetic Ising model. The result is known more generally, but the proofs are much easier for the Ising model.

People have known about the Metropolis-Hastings algorithm since the 1950s, but only recently have researchers been able to prove convergence guarantees for the 2D Ising model. There is a large gap between theory and practice, but in some situations we can prove that the algorithm works.

Sampling from the distribution is roughly equivalent to estimating the partition function (sampling-counting equivalence). There have been many papers addressing tasks such as estimating the non-negative permanent, the number of colorings of a graph, etc.2 A dominant way of accomplishing these tasks is proving that the Metropolis-Hastings algorithm converges for these problems. It is easy to find algorithms for these problems that converge to Gibbs distributions, but the convergence may take exponential time.

We will look at the situation when the energy function looks like the Ising model, in the sense that the interactions are local and reflect the structure of some underlying space. Also, assume that the interactions are of size O(1)O(1)

O(1) and that the scaling comes from the size of the system. When can we expect that our algorithms work? There are two main cases when we can argue that there should be rapid mixing.

High temperature regime: The system is very disordered, and in the limit as the temperature approaches infinity, we get the uniform distribution.

One-dimension: In 1D, we can exactly compute the partition function using dynamic programming. Before, we mentioned that if there are a sea of ++

+s and an island of −−

-s, then it is energetically favorable for the island to shrink; note that this is no longer true in 1D. In a way, 1D systems are more “boring” because they cannot exhibit arbitrarily long-range correlations.

In this part of the blog post, we will try to be more proof-oriented. We will start by explaining why it is plausible that high temperature means that the chain will mix rapidly in time.

Coupling method

One method of proving rates of convergence for Markov chains is by analzying the spectral gap. Another method is the coupling method.

The idea behind the coupling method is to start with two configurations X(0),Y(0)∈ΩX(0),Y(0)∈Ω

X(0),Y(0) \in \Omega. We want each one to evolve under the Markov chain.

The key part is that there is still some freedom with respect to what the dynamics looks like. In particular, we are allowed to correlate the XX

X and YY

Y processes. Thus, we are defining a joint transition probability P{X(1)=x(1),Y(1)=y(1)∣X(0),Y(0)}P{X(1)=x(1),Y(1)=y(1)∣X(0),Y(0)}

{\mathbb P}\{X(1)=x(1),Y(1)=y(1) \mid X(0),Y(0)\}. We want to design the process such that X(1)X(1)

X(1) and Y(1)Y(1)

Y(1) are closer together than X(0)X(0)

X(0) and Y(0)Y(0)

Y(0). Imagine that we have two particles bouncing around. Each particle follows the dynamics of TT

T, but they are correlated so that they drift together, and once they meet, they stick together. It turns out that the mixing time can be upper bounded by the time it takes for the particles to meet each other.

Assume we have some sort of distance function distdist

\operatorname{dist} on the underlying space and we can prove that Edist(X(1),Y(1))≤exp(−α)dist(X(0),Y(0))E⁡dist⁡(X(1),Y(1))≤exp⁡(−α)dist⁡(X(0),Y(0))

\operatorname{\mathbb E}\operatorname{dist}(X(1),Y(1)) \le \exp(-\alpha) \operatorname{dist}(X(0),Y(0)). Then, it turns out that the mixing time tmix(ϵ)tmix(ϵ)

t_{\rm mix}(\epsilon), i.e. the time required to get within ϵϵ

\epsilon of the stationary distribution, is upper bounded as

tmix(ϵ)≤log{(diamΩ)/ϵ}αtmix(ϵ)≤log⁡{(diam⁡Ω)/ϵ}α

t_{\rm mix}(\epsilon) \le \frac{\log\{(\operatorname{diam}\Omega)/\epsilon\}}{\alpha}

Initially, the two particles can be diamΩdiam⁡Ω

\operatorname{diam}\Omega apart, but the expected distance is exponentially shrinking as we run the coupling, so the mixing time is logarithmic in the diameter.

The distance between probability distributions is defined as follows. Let pp

p and qq

q be two probability distributions on ΩΩ

\Omega. Then, the metric is:3

<span id="MathJax-Element-218-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" style="font-size: 110%; text-align: center; position: relative;" data-mathml="12‖p−q‖1=12∑x∈Ω|p(x)−q(x)|=min(X,Y)∼r∈P(Ω×Ω)r1=pr2=qPr{X≠Y}.” role=”presentation”>12∥p−q∥1=12∑x∈Ω|p(x)−q(x)|=min(X,Y)∼r∈P(Ω×Ω)r1=pr2=qPr{X≠Y}.12‖p−q‖1=12∑x∈Ω|p(x)−q(x)|=min(X,Y)∼r∈P(Ω×Ω)r1=pr2=qPr{X≠Y}.

\frac{1}{2} \|p-q\|_1 = \frac{1}{2}\sum_{x \in \Omega}|p(x)-q(x)| = \min_{\substack{(X,Y) \sim r \in {\mathcal{P}}(\Omega \times \Omega) \\ r_1 = p \\ r_2 = q}} {\mathbb P}_r\{X \ne Y\}.

In this expression, r1r1

r_1 and r2r2

r_2 denote the first and second marginals of rr

r respectively. The minimum is taken over all couplings of pp

p and qq

q. This is the correct way to measure the distance between distributions. To give some intuition for this quantity, the quantity on the right represents the best test to distinguish the two distributions. If pp

p and qq

q are the same, we can take a coupling in which X∼pX∼p

X \sim p and Y∼qY∼q

Y \sim q are always identical. If pp

p and qq

q have disjoint supports, then no matter what coupling we use, XX

X and YY

Y will never be equal.

It suffices to consider when X(0)X(0)

X(0) and Y(0)Y(0)

Y(0) are neighbors, i.e. at distance 11

1 apart. This is because if we have X(0)X(0)

X(0) and Y(0)Y(0)

Y(0) far apart, then we could look at the path between them and reduce to the case when they are neighbors. Formally, this is known as path coupling. The formal statement is in Theorem 12.3 of [4]:

Theorem 3: Let ΓΓ

\Gamma be a connected weighted graph on the state space, where no edge has weight less than dmindmin

d_{\min}. Let d(C,C′)d(C,C′)

d(C,C’) be the length of the shortest path from CC

C to C′C′

C’ in ΓΓ

\Gamma and let dmax=maxC,C′∈Ωd(C,C′)dmax=maxC,C′∈Ωd(C,C′)

d_{\max} = \max_{C,C’ \in \Omega} d(C,C’) be the diameter of ΓΓ

\Gamma. Suppose there is a coupling such that for some 0″ role=”presentation” style=”font-size: 110%; position: relative;”>δ>0δ>0

\delta > 0,

<span id="MathJax-Element-246-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" style="font-size: 110%; text-align: center; position: relative;" data-mathml="E⁡[d(X(1),Y(1))|(X(0),Y(0))=(C,C′)]≤(1−δ)d(C,C′)” role=”presentation”>E[d(X(1),Y(1))∣∣(X(0),Y(0))=(C,C′)]≤(1−δ)d(C,C′)E⁡[d(X(1),Y(1))|(X(0),Y(0))=(C,C′)]≤(1−δ)d(C,C′)

\operatorname{\mathbb E}\bigl[d\bigl(X(1),Y(1)\bigr) \bigm\vert \bigl(X(0),Y(0)\bigr) = (C,C’)\bigr] \le (1-\delta)d(C,C’)

for all neighboring pairs CC

C, C′C′

C’, i.e., those pairs connected by an edge in ΓΓ

\Gamma. Then, the mixing time is bounded by

<span id="MathJax-Element-250-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" style="font-size: 110%; text-align: center; position: relative;" data-mathml="tmix(ϵ)≤log⁡(ϵ−1dmax/dmin)δ.” role=”presentation”>tmix(ϵ)≤log(ϵ−1dmax/dmin)δ.tmix(ϵ)≤log⁡(ϵ−1dmax/dmin)δ.

t_{\rm mix}(\epsilon) \le \frac{\log(\epsilon^{-1}d_{\max}/d_{\min})}{\delta}.

Glauber dynamics at high temperature

Recall that in Glauber dynamics, we pick a site ii

i randomly and then update the site conditioned on its neighbors. The first way we will couple together X(1)X(1)

X(1) and Y(1)Y(1)

Y(1) is by picking the same site for both of them.

Pick a random i∈[n]i∈[n]

i \in [n].

If X(0)N(i)=Y(0)N(i)X(0)N(i)=Y(0)N(i)

{X(0)}_{N(i)} = {Y(0)}_{N(i)}, then set X(1)i=Y(1)iX(1)i=Y(1)i

{X(1)}_i = {Y(1)}_i (if the neighborhoods of the two points agree, then update them the same way). Otherwise, update them using the best possible coupling, i.e., pick a coupling for (X(1)i,Y(1)i)(X(1)i,Y(1)i)

({X(1)}_i, {Y(1)}_i) which minimizes P{X(1)i≠Y(1)i}P{X(1)i≠Y(1)i}

{\mathbb P}\{ {X(1)}_i \ne {Y(1)}_i \}.

So if X(0)=Y(0)X(0)=Y(0)

X(0) = Y(0), then the points will never drift apart. The reason why analyzing this coupling is non-trivial is because there is a chance that the distance between the two points can increase.

Assume that the degree of the graph is ΔΔ

\Delta. Suppose that dist(X(0),Y(0))=1dist⁡(X(0),Y(0))=1

\operatorname{dist}(X(0),Y(0)) = 1, that is, there is a single aa

a such that X(0)a≠Y(0)aX(0)a≠Y(0)a

{X(0)}_a \ne {Y(0)}_a. What will happen to X(1)X(1)

X(1) and Y(1)Y(1)

Y(1)? We start by picking a random i∈[n]i∈[n]

i \in [n]. There are three cases:

i∉({a}∪N(a))i∉({a}∪N(a))

i \notin (\{a\} \cup N(a)) (with probability 1−(Δ+1)/n1−(Δ+1)/n

1 – (\Delta + 1)/n): Nothing changes; X(0)X(0)

X(0) and Y(0)Y(0)

Y(0) agree at ii

i, and X(1)X(1)

X(1) and Y(1)Y(1)

Y(1) will also agree at ii

i. The distance remains at 11

i=ai=a

i = a (with probability 1/n1/n

1/n): We picked the one spot in which the two configurations differ. The neighborhoods of aa

a are the same for X(0)X(0)

X(0) and Y(0)Y(0)

Y(0), so we update in the same way for both processes, and the distance drops to 00

i∈N(a)i∈N(a)

i \in N(a) (with probability Δ/nΔ/n

\Delta/n): We could have different updates. Here, we have to use the high temperature assumption, which says that if we change one bit, the probability of a configuration cannot change too much.

In the Ising model, E(x)=∑ni,j=1Ji,jxixjE(x)=∑i,j=1nJi,jxixj

E(x) = \sum_{i,j=1}^n J_{i,j} x_i x_j. Changing aa

a can bias the energy by at most ΔmaxiJi,aΔmaxiJi,a

\Delta\max_i J_{i,a}, so the expected distance afterwards is 1+O(maxni,j=1Ji,j/T)1+O(maxi,j=1nJi,j/T)

1 + O(\max_{i,j=1}^n J_{i,j}/T).

Adding these cases up to get the overall expected distance gives

<span id="MathJax-Element-288-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" style="font-size: 110%; text-align: center; position: relative;" data-mathml="E⁡dist⁡(X(1),Y(1))=1−1n+O(ΔJmaxT)⏟≤11n=1−cn” role=”presentation”>Edist(X(1),Y(1))=1−1n+O(ΔJmaxT)≤11n=1−cnE⁡dist⁡(X(1),Y(1))=1−1n+O(ΔJmaxT)⏟≤11n=1−cn

\operatorname{\mathbb E}\operatorname{dist}\bigl(X(1), Y(1)\bigr) = 1-\frac{1}{n} + \underbrace{O\Bigl(\frac{\Delta J_{\max}}{T}\Bigr)}_{\le 1}\frac{1}{n} = 1 – \frac{c}{n}

for TT

T large enough, so the expected distance will shrink. This argument also tells us how large the temperature must be, which is important for applications. This gives us tmix(ϵ)=O(nlognϵ).tmix(ϵ)=O(nlog⁡nϵ).

t_{\rm mix}(\epsilon) = O\Bigl(n\log\frac{n}{\epsilon}\Bigr). Notice that this is the same dependence as the coupon collector problem. Therefore, in the high temperature regime, the system behaves qualitatively as if there are no correlations.

Temporal and spatial mixing equivalence

The analysis of Glauber dynamics at high temperature is already a version of the equivalence between mixing in time and mixing in space. It says that if the correlations even with the immediate neighbors of a node are weak, then Glauber dynamics rapidly mixes.

Now, we want to consider the situation in which there can be strong correlations between immediate neighbors, but weak correlation with far away sites. We want to show that spatial mixing implies temporal mixing.

We will give a few definitions of correlation decay. (Note: The definitions of correlation decay below are not exactly the ones from Aram’s lecture. These definitions are from [5] and [6].)

For non-empty W⊆VW⊆V

W \subseteq V and τ∈ΣV∖Wτ∈ΣV∖W

\tau \in \Sigma^{V\setminus W}, let μτWμWτ

\mu_W^\tau be the distribution of the spins in WW

W conditional on the spins in V∖WV∖W

V \setminus W being fixed to ττ

\tau. For Δ⊆WΔ⊆W

\Delta \subseteq W, let μτW,ΔμW,Δτ

\mu_{W,\Delta}^\tau be the marginal of μτWμWτ

\mu_W^\tau on the spins in ΔΔ

\Delta. We will assume that the interactions between the spins have finite range 0″ role=”presentation” style=”font-size: 110%; position: relative;”>r>0r>0

r > 0, and ∂rW∂rW

\partial_r W denotes the rr

r-boundary of WW

W, i.e., {v∈V∖W:dist(v,W)≤r}{v∈V∖W:dist⁡(v,W)≤r}

\{v \in V \setminus W : \operatorname{dist}(v,W) \le r\}.

(Weak decay of correlations) Weak spatial mixing holds for WW

W if there exist constants 0″ role=”presentation” style=”font-size: 110%; position: relative;”>C,ξ>0C,ξ>0

C, \xi > 0 such that for any subset Δ⊆WΔ⊆W

\Delta \subseteq W, <span id="MathJax-Element-309-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" data-mathml="supτ,τ′∈ΣV∖W‖μW,Δτ−μW,Δτ′‖1≤C∑x∈Δ,y∈∂rWexp⁡(−dist⁡(x,y)ξ).” role=”presentation” style=”font-size: 110%; position: relative;”>supτ,τ′∈ΣV∖W∥μτW,Δ−μτ′W,Δ∥1≤C∑x∈Δ,y∈∂rWexp(−dist(x,y)ξ).supτ,τ′∈ΣV∖W‖μW,Δτ−μW,Δτ′‖1≤C∑x∈Δ,y∈∂rWexp⁡(−dist⁡(x,y)ξ).

\sup_{\tau,\tau’ \in \Sigma^{V\setminus W}}\|\mu_{W,\Delta}^\tau – \mu_{W,\Delta}^{\tau’}\|_1 \le C\sum_{x\in\Delta, \; y \in \partial_r W} \exp\Bigl(- \frac{\operatorname{dist}(x,y)}{\xi}\Bigr).

(Strong decay of correlations) Strong spatial mixing holds for WW

W if there exist constants 0″ role=”presentation” style=”font-size: 110%; position: relative;”>C,ξ>0C,ξ>0

C,\xi > 0 such that for every Δ⊆WΔ⊆W

\Delta \subseteq W and every τ,τ′∈ΣV∖Wτ,τ′∈ΣV∖W

\tau,\tau’ \in \Sigma^{V\setminus W} differing only at site y∈V∖Wy∈V∖W

y \in V\setminus W, ∥μτW,Δ−μτ′W,Δ∥1≤Cexp(−dist(y,Δ)ξ).‖μW,Δτ−μW,Δτ′‖1≤Cexp⁡(−dist⁡(y,Δ)ξ).

\|\mu_{W,\Delta}^\tau – \mu_{W,\Delta}^{\tau’}\|_1 \le C\exp\Bigl(-\frac{\operatorname{dist}(y,\Delta)}{\xi}\Bigr).

(Strong decay of correlations) Strong spatial mixing in the truncated sense holds for VV

V if there exist 0″ role=”presentation” style=”font-size: 110%; position: relative;”>n,ξ>0n,ξ>0

n, \xi > 0 such that for all functions f,g:Ω→Rf,g:Ω→R

f, g : \Omega \to {\mathbb R} which depend only on the sites at ΛfΛf

\Lambda_f and ΛgΛg

\Lambda_g respectively and such that dist(Λf,Λg)≥ndist⁡(Λf,Λg)≥n

\operatorname{dist}(\Lambda_f,\Lambda_g) \ge n, <span id="MathJax-Element-322-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" data-mathml="supτ∈ΣV∖WcovμWτ⁡(f,g)≤|Λf||Λg|‖f‖∞‖g‖∞exp⁡(−d(Λf,Λg)ξ).” role=”presentation” style=”font-size: 110%; position: relative;”>supτ∈ΣV∖WcovμτW(f,g)≤|Λf||Λg|∥f∥∞∥g∥∞exp(−d(Λf,Λg)ξ).supτ∈ΣV∖WcovμWτ⁡(f,g)≤|Λf||Λg|‖f‖∞‖g‖∞exp⁡(−d(Λf,Λg)ξ).

\sup_{\tau \in \Sigma^{V\setminus W}} \operatorname{cov}_{\mu_W^\tau}(f, g) \le |\Lambda_f||\Lambda_g|\|f\|_\infty \|g\|_\infty \exp\Bigl(-\frac{d(\Lambda_f,\Lambda_g)}{\xi}\Bigr).

Here, ξξ

\xi is the correlation length (in physics, it is the characteristic length scale of a system). In the disordered phase, the correlation length is a constant independent of system size. For our purposes, the main consequence of these definitions is that the effective interaction range of each spin is O(ξ)O(ξ)

O(\xi). For the Ising model, there is a key simplification due to monotonicity. Namely, the ferromagnetic Ising model has the nice property (which is not true for other models) that if we flip a sign from −−

– to ++

+, this only makes ++

+ more likely everywhere. This is because the spins want to agree. There are a lot of boundary conditions to consider, but here, due to monotonicity, we only need to consider two: all of the spins are −−

-, and all of the spins are ++

+. All ++

+ spins will give the highest probability of a ++

+ spin, and all −−

– spin will give the lowest probability of a ++

+ spin. This monotonicity property is generally not required for time-space mixing equivalence to hold, but it greatly simplifies proofs.

It is a very non-obvious fact that all of these notions of spatial mixing are equivalent. We will sketch a proof that strong correlation decay implies that tmix=O(nlogn)tmix=O(nlog⁡n)

t_{\rm mix} = O(n\log n).

The idea is to use another coupling argument. Let X(0),Y(0)∈{±1}VX(0),Y(0)∈{±1}V

X(0), Y(0) \in {\{\pm 1\}}^V differ in one coordinate, i.e., X(0)a≠Y(0)aX(0)a≠Y(0)a

{X(0)}_a \ne {Y(0)}_a and X(0)i=Y(0)iX(0)i=Y(0)i

{X(0)}_i = {Y(0)}_i for i≠ai≠a

i \ne a. We want to argue that the expected distance between the processes will decrease. The proof uses a generalization of Glauber dynamics called block Glauber dynamics. In Glauber dynamics, we take a single spin and resample it conditioned on its neighbors. In block Glauber dynamics, we take an L×LL×L

L\times L box and resample it conditioned on its neighbors. There is an argument, called canonical paths, which can be used to show that if block Glauber dynamics mixes, then regular Glauber dynamics also mixes (slightly more slowly; we lose a poly(L)poly⁡(L)

\operatorname{poly}(L) factor, but anyway LL

L will be a large constant) so analyzing block Glauber dynamics is fine.

If aa

a lies in the box, then the expected change in distance is −L2/n−L2/n

-L^2/n. If aa

a is far away from the box, then there is no change. If aa

a is in the boundary of the box, then it is possible for the distance to increase. However, strong spatial mixing allows us to control the influence of a single site, so the expected change in distance is bounded by O(Lξ2/n)O(Lξ2/n)

O(L\xi^2/n). Now, since ξξ

\xi is a constant, if we choose LL

L sufficiently large, then we will have the same situation as in the high temperature case: the expected distance will exponentially shrink over time.

Quantum systems

The quantum version of Markov chains has many more difficulties. The first difficulty is that the Hammersley-Clifford theorem (which we have been relying on throughout this blog post) fails.

Notation

To properly discuss what we mean, let’s set up some notation. Readers already familiar with density matrices, quantum entropy, and quantum mutual information may wish to skip to the next subsection. Most of the time we discuss quantum objects here, we’ll be using density matricies, often denoted ρρ

\rho. A density matrix can be thought of as an extension to regular quantum states |ψ⟩|ψ⟩

\ket{\psi}, where there is some classical source of uncertainty.

A density matrix is a positive semidefinite matrix with trace 11

1. This extends the notion of a classical probability distribution; in the quantum setting, a classical probability distribution pp

p (thought of as a vector whose entries sum to 11

1) is represented as the density matrix diag(p)diag⁡(p)

\operatorname{diag}(p).

For example, we can consider a situation in which there is a 1/21/2

1/2 probability that we started with the quantum state |ψ⟩|ψ⟩

\ket{\psi} and a 1/21/2

1/2 probability that we started with the quantum state |ϕ⟩|ϕ⟩

\ket{\phi}. This would be denoted as follows:

ρ=12|ψ⟩⟨ψ|+12|ϕ⟩⟨ϕ|ρ=12|ψ⟩⟨ψ|+12|ϕ⟩⟨ϕ|

\rho = \frac{1}{2} \ket{\psi} \bra{\psi} + \frac{1}{2} \ket{\phi} \bra{\phi}

Density matricies are generally useful for a lot of tasks, but for our purposes a density matrix will be used to discuss both the classical and quantum “uncertainty” we have about what state we have.

Now let’s also talk about a second important piece of notation: the tensor product. Often when discussing quantum states, it is important to discuss multiple quantum states simultaneously. For example, Alice has one system AA

A and Bob has another system BB

B. However, these systems might be entangled, meaning that the results of the two systems are correlated.

For instance, let us consider the following state:

|ψ⟩=1√2(|+⟩A|+⟩B+|−⟩A|−⟩B)|ψ⟩=12(|+⟩A|+⟩B+|−⟩A|−⟩B)

\ket{\psi} = \frac{1}{\sqrt{2}}\left(\ket{+}_A \ket{+}_B + \ket{-}_A \ket{-}_B \right)

This particular state has the property that Alice and Bob will always both measure ++

+ or they will both measure −−

-. The notation for tensors is often ambiguous in the literature as there are many ways of specifying tensors. For instance, above we used subscripts to explicitly denote which particle was in system AA

A and which was in system BB

B. One may also choose to simply use the index of the system as below. The symbol ⊗⊗

\otimes is used to denote a tensor between states (where it is assumed that the first state is system AA

A and the second, system BB

B). Gradually folks may shorten the notation as follows:

<span id="MathJax-Element-370-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" style="font-size: 110%; text-align: center; position: relative;" data-mathml="|ψ⟩=12(|+⟩|+⟩+|−⟩|−⟩)|ψ⟩=12(|++⟩+|−−⟩)|ψ⟩=12(10)⊗(01)” role=”presentation”>|ψ⟩=1√2(|+⟩|+⟩+|−⟩|−⟩)|ψ⟩=1√2(|++⟩+|−−⟩)|ψ⟩=1√2(10)⊗(01)|ψ⟩=12(|+⟩|+⟩+|−⟩|−⟩)|ψ⟩=12(|++⟩+|−−⟩)|ψ⟩=12(10)⊗(01)

These are all notations for the same state. Let’s now talk about this state in the context of a density matrix. The density matrix of this state is as follows:

<span id="MathJax-Element-371-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" style="font-size: 110%; text-align: center; position: relative;" data-mathml="ρA,B=12(|++⟩+|−−⟩)(⟨++|+⟨−−|)ρA,B=12(|++⟩⟨++|+|−−⟩⟨++|+|++⟩⟨−−|+|−−⟩⟨−−|)ρA,B=12(1001000000001001)” role=”presentation”>ρA,B=12(|++⟩+|−−⟩)(⟨++|+⟨−−|)ρA,B=12(|++⟩⟨++|+|−−⟩⟨++|+|++⟩⟨−−|+|−−⟩⟨−−|)ρA,B=12⎛⎜ ⎜ ⎜⎝1001000000001001⎞⎟ ⎟ ⎟⎠ρA,B=12(|++⟩+|−−⟩)(⟨++|+⟨−−|)ρA,B=12(|++⟩⟨++|+|−−⟩⟨++|+|++⟩⟨−−|+|−−⟩⟨−−|)ρA,B=12(1001000000001001)

Writing the density matrix ρρ

\rho as ρA,BρA,B

\rho_{A,B} makes explicit that this is the density matrix over systems AA

A and BB

A crucial operation that one will often perform using density matricies is the partial trace. The partial trace is a way of allowing us to consider only a smaller part of the larger part of the system, while taking into account the influence of the larger system around it.

Here’s an example: Suppose Bob wants to know what his state is. However, Bob really doesn’t care about Alice’s system and just wants to know what the density matrix for his system is. Bob’s density matrix is simply the following density matrix (a 50% chance of being in |+⟩|+⟩

\ket{+} and a 50% chance of being in |−⟩|−⟩

\ket{-}).

ρB=12(|+⟩⟨+|+|−⟩⟨−|)ρB=12(|+⟩⟨+|+|−⟩⟨−|)

More explicitly, we could write the following:

ρB=12(|+⟩B⟨+|B+|−⟩B⟨−|B)ρB=12(|+⟩B⟨+|B+|−⟩B⟨−|B)

The partial trace is an operation that will let us take our original density matrix ρA,BρA,B

\rho_{A,B} and generates a new density matrix ρBρB

\rho_B that ignores system AA

A. This is specifically called the partial trace over AA

A, or trAtrA

\operatorname{tr}_A.

So how do we do this? We simply sum over the state AA

A (effectively taking a trace, but only along one axis):

<span id="MathJax-Element-386-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" style="font-size: 110%; text-align: center; position: relative;" data-mathml="trA⁡ρA,B=∑i⟨i|A12(|++⟩⟨++|+|−−⟩⟨++|+|++⟩⟨−−|+|−−⟩⟨−−|)|i⟩A” role=”presentation”>trAρA,B=∑i⟨i|A12(|++⟩⟨++|+|−−⟩⟨++|+|++⟩⟨−−|+|−−⟩⟨−−|)|i⟩AtrA⁡ρA,B=∑i⟨i|A12(|++⟩⟨++|+|−−⟩⟨++|+|++⟩⟨−−|+|−−⟩⟨−−|)|i⟩A

This is easier to evaluate using certain choices of notation:

<span id="MathJax-Element-387-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" style="font-size: 110%; text-align: center; position: relative;" data-mathml="trA⁡ρA,B=⟨+|A12(|++⟩⟨++|+|−−⟩⟨++|+|++⟩⟨−−|+|−−⟩⟨−−|)|+⟩A+⟨−|A12(|++⟩⟨++|+|−−⟩⟨++|+|++⟩⟨−−|+|−−⟩⟨−−|)|−⟩A=12(|+⟩B⟨++|+|+⟩B⟨−−|)|+⟩A+12(|−⟩B⟨++|+|−⟩B⟨−−|)|−⟩A=12(|+⟩B⟨+|B)+12(|−⟩B⟨−|B)=12(|+⟩B⟨+|B+|−⟩B⟨−|B)=ρB” role=”presentation”>trAρA,B=⟨+|A12(|++⟩⟨++|+|−−⟩⟨++|+|++⟩⟨−−|+|−−⟩⟨−−|)|+⟩A+⟨−|A12(|++⟩⟨++|+|−−⟩⟨++|+|++⟩⟨−−|+|−−⟩⟨−−|)|−⟩A=12(|+⟩B⟨++|+|+⟩B⟨−−|)|+⟩A+12(|−⟩B⟨++|+|−⟩B⟨−−|)|−⟩A=12(|+⟩B⟨+|B)+12(|−⟩B⟨−|B)=12(|+⟩B⟨+|B+|−⟩B⟨−|B)=ρBtrA⁡ρA,B=⟨+|A12(|++⟩⟨++|+|−−⟩⟨++|+|++⟩⟨−−|+|−−⟩⟨−−|)|+⟩A+⟨−|A12(|++⟩⟨++|+|−−⟩⟨++|+|++⟩⟨−−|+|−−⟩⟨−−|)|−⟩A=12(|+⟩B⟨++|+|+⟩B⟨−−|)|+⟩A+12(|−⟩B⟨++|+|−⟩B⟨−−|)|−⟩A=12(|+⟩B⟨+|B)+12(|−⟩B⟨−|B)=12(|+⟩B⟨+|B+|−⟩B⟨−|B)=ρB

This gives us the answer that we had expected.

We now have all of the tools we need to talk about quantum entropy. Intuitively, entropy can be thought of as the amount of uncertainty we have for our system, or equivalently the amount of information it takes to define our system. The entropy for a quantum system ρρ

\rho is defined as follows:

H(ρ)=−tr(ρlog2ρ)H(ρ)=−tr⁡(ρlog2⁡ρ)

Note that here we use the shorthand ρBρB

\rho_B to denote trAρA,BtrA⁡ρA,B

\operatorname{tr}_A \rho_{A,B}. Here, writing trtr

\operatorname{tr} without the subscript indicates that this is the full or normal trace that one might expect (or equivalently performing the partial trace over all systems). We can now define the conditional entropy of a system as follows:

H(A∣B)ρ=H(ρA,B)−H(ρB)H(A∣B)ρ=H(ρA,B)−H(ρB)

This definition intuitively makes sense since we can think of conditional entropy as the amount of information it takes to describe our joint system (A,B)(A,B)

(A,B), given that we already know what BB

B is.

We can now discuss quantum mutual information, the amount of information that measuring system AA

A will provide you about system BB

B. Like the classical case, this is defined as follows:

<span id="MathJax-Element-398-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" style="font-size: 110%; text-align: center; position: relative;" data-mathml="I(A;B)ρ=H(A,B)ρ−H(A∣B)ρ−H(B∣A)ρ” role=”presentation”>I(A;B)ρ=H(A,B)ρ−H(A∣B)ρ−H(B∣A)ρI(A;B)ρ=H(A,B)ρ−H(A∣B)ρ−H(B∣A)ρ

We can now finally discuss quantum mutual information (QCMI), defined as follows: I(A;B∣C)ρ=I(A;B,C)ρ−I(A;C)ρI(A;B∣C)ρ=I(A;B,C)ρ−I(A;C)ρ

{I(A;B \mid C)}_\rho = {I(A;B,C)}_\rho – {I(A;C)}_\rho. With some algebraic simplifications, one can arrive at the expression:

<span id="MathJax-Element-400-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" style="font-size: 110%; text-align: center; position: relative;" data-mathml="I(A;B∣C)ρ=H(A,C)ρ+H(B,C)ρ−H(A,B,C)ρ−H(C)ρ.” role=”presentation”>I(A;B∣C)ρ=H(A,C)ρ+H(B,C)ρ−H(A,B,C)ρ−H(C)ρ.I(A;B∣C)ρ=H(A,C)ρ+H(B,C)ρ−H(A,B,C)ρ−H(C)ρ.

The QCMI equals 00

0 if and only if ρρ

\rho is a quantum Markov state. Classically, the entropic characterization of conditional independence corresponds to an algebraic characterization.

Recovery Maps

Here, the algebraic characterization is more grueling. We have

<span id="MathJax-Element-403-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" style="font-size: 110%; text-align: center; position: relative;" data-mathml="ρABC=exp⁡(log⁡ρAB+log⁡ρBC−log⁡ρB)” role=”presentation”>ρABC=exp(logρAB+logρBC−logρB)ρABC=exp⁡(log⁡ρAB+log⁡ρBC−log⁡ρB)

\rho_{ABC} = \exp(\log \rho_{AB} + \log \rho_{BC} – \log \rho_B)

Equivalently,

<span id="MathJax-Element-404-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" style="font-size: 110%; text-align: center; position: relative;" data-mathml="ρABC=ρAB1/2ρB−1/2ρBCρB−1/2ρAB1/2=RB→AB(ρBC)” role=”presentation”>ρABC=ρ1/2ABρ−1/2BρBCρ−1/2Bρ1/2AB=RB→AB(ρBC)ρABC=ρAB1/2ρB−1/2ρBCρB−1/2ρAB1/2=RB→AB(ρBC)

\rho_{ABC} = \rho_{AB}^{1/2} \rho_B^{-1/2} \rho_{BC}\rho_B^{-1/2} \rho_{AB}^{1/2} = R_{B\to AB}(\rho_{BC})

Here, RB→ABRB→AB

R_{B\to AB} is called the Petz recovery map,4 ρB→AB(X)=ρ1/2ABρ−1/2BXρ−1/2Bρ1/2ABρB→AB(X)=ρAB1/2ρB−1/2XρB−1/2ρAB1/2

\rho_{B\to AB}(X) = \rho_{AB}^{1/2} \rho_B^{-1/2} X\rho_B^{-1/2} \rho_{AB}^{1/2}. One can think of a recovery may as a way that we can reconstruct the entire system A,BA,B

A, B using just system BB

B. It is not obvious that this is a quantum channel, but it is.

Suppose ρρ

\rho is a probability distribution, so ρ=diag(p)ρ=diag⁡(p)

\rho =\operatorname{diag}(p) for some vector pp

p. Then, all of the density matrices are diagonal and commuting. Then, the recovery map means that we divide by pBpB

p_B and multiply by pABpAB

p_{AB}, i.e., multiply by pA∣BpA∣B

p_{A \mid B}. This is the natural thing to do if we lost our information about AA

A and were trying to figure out what AA

A was based on our knowledge of BB

B. This is why RB→A,BRB→A,B

R_{B\to A,B} is known as a recovery map, and it is used to discuss conditional distributions in the quantum setting. In the classical case, if we start with B,CB,C

B, C, look only at BB

B, and use this to reconstruct AA

A, then we would have the whole state in a Markov chain. That is why this is a plausible quantum version of being a Markov chain.

However, quantum Gibbs states are not, in general, quantum Markov chains. The failure of this statement to hold is related to topological order, which is similar to the degrees of freedom that show up in error correcting codes.

Quantum Markov Networks

Here, we will formally define a quantum Markov network. The reference for this is [7].

Let G=(V,E)G=(V,E)

G = (V, E) be a finite graph. We associate with each vertex v∈Vv∈V

v \in V a Hilbert space HvHv

{\mathcal{H}}_v and we consider a density matrix ρVρV

\rho_V acting on ⨂v∈VHv⨂v∈VHv

\bigotimes_{v\in V} {\mathcal{H}}_v. Then, (G,ρV)(G,ρV)

(G, \rho_V) is a quantum Markov network if for all U⊆VU⊆V

U\subseteq V, UU

U is conditionally independent of V∖(U∪∂U)V∖(U∪∂U)

V \setminus (U \cup \partial U) given ∂U∂U

\partial U, where the conditional independence statement is w.r.t. ρVρV

\rho_V and means that the corresponding QCMI satisfies I(U;V∖(U∪∂U)∣∂U)ρV=0I(U;V∖(U∪∂U)∣∂U)ρV=0

{I(U; V\setminus (U \cup \partial U) \mid \partial U)}_{\rho_V} = 0.

A quantum Markov network is called positive if ρVρV

\rho_V has full rank. (Recall that in the statement of the Hammersley-Clifford Theorem, , it is assumed that the distribution is strictly positive.)

Now, consider the following example. First, we introduce the Pauli matrices

<span id="MathJax-Element-435-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" style="font-size: 110%; text-align: center; position: relative;" data-mathml="σx:=[0110],σz:=[100−1],σy:=[0−ii0].” role=”presentation”>σx:=[0110],σz:=[100−1],σy:=[0−ii0].σx:=[0110],σz:=[100−1],σy:=[0−ii0].

We define a Hamiltonian on three qubits AA

A, BB

B, CC

C by

<span id="MathJax-Element-439-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" style="font-size: 110%; text-align: center; position: relative;" data-mathml="H:=(σAxσBx+σAyσBy+σAzσBz)IC+IA(σBxσCx+σByσCy+σBzσCz)” role=”presentation”>H:=(σxAσxB+σyAσyB+σzAσzB)IC+IA(σxBσxC+σyBσyC+σzBσzC)H:=(σAxσBx+σAyσBy+σAzσBz)IC+IA(σBxσCx+σByσCy+σBzσCz)

H := (\sigma_A^x \sigma_B^x + \sigma_A^y \sigma_B^y + \sigma_A^z \sigma_B^z) I_C + I_A (\sigma_B^x \sigma_C^x + \sigma_B^y \sigma_C^y + \sigma_B^z \sigma_C^z)

(Juxtaposition in the above expression signifies the tensor product as discussed before.) Finally, for 0″ role=”presentation” style=”font-size: 110%; position: relative;”>β>0β>0

\beta > 0, we define the Gibbs state

ρA,B,C(β):=1Z(β)exp(−βH)ρA,B,C(β):=1Z(β)exp⁡(−βH)

\rho_{A,B,C}(\beta) := \frac{1}{Z(\beta)} \exp(-\beta H)

The Hamiltonian here has local terms which correspond to interactions (A,B)(A,B)

(A,B), (B,C)(B,C)

(B, C). However, it can be shown that the QCMI between AA

A and CC

C conditioned on BB

B w.r.t. ρA,B,CρA,B,C

\rho_{A,B,C} is non-zero, which means that this is not a quantum Markov network w.r.t. the line graph <span id="MathJax-Element-448-Frame" class="mjx-chtml MathJax_CHTML" tabindex="0" data-mathml="A

C” role=”presentation” style=”font-size: 110%; position: relative;”>A↔B↔CA↔B↔C

A \leftrightarrow B \leftrightarrow C. This demonstrates the failure of the Hammersley-Clifford Theorem in the quantum setting.

Important Results

We will briefly discuss the results of two papers.

[8] This paper shows that mixing in space implies mixing in time in the quantum case. However, the result of the paper only applies to commuting Hamiltonians. For commuting Hamiltonians, it turns out that quantum Gibbs states are quantum Markov networks. They use a version of Glauber dynamics, which can be simulated on a quantum computer but are also plausible dynamics for a physical system in nature. This is a difficult paper to read, but it is worth digesting if you want to work in the field.

[9] This second paper is much easier and more general, covering non-commuting Hamiltonians, but it requires more conditions. They give a method of preparing the Gibbs state which can run on a quantum computer, but the dynamics are not plausible as a physical system because they are too complicated. The more complicated dynamics allows them to make the proof work. The paper also uses QCMI.

They have two assumptions. The first assumption looks like mixing in space (weak correlation decay). The second assumption is that the state looks approximately like a quantum Markov network (this is definitely not met in general). A very important paper in this space is a recent breakthrough ([10]) which characterizes quantum Markov chains. They show that if the QCMI is bounded by ϵϵ

\epsilon, then the recovery map RB→A,B(ρBC)RB→A,B(ρBC)

R_{B\to A,B}(\rho_{BC}) is ϵ′ϵ′

\epsilon’-close to ρABCρABC

\rho_{ABC}, i.e., low QCMI implies that the recovery map works well. This is trivial to prove classically, but very difficult in the quantum world.

The algorithm in [9] is very elegant. Essentially, we take the entire system and punch out constant-sized boxes. If we can reconstruct the region outside of the boxes, then we can use the recovery maps to reconstruct the regions inside of the boxes, and the boxes are far apart enough so they are almost independent. For this argument, we must assume that the QCMI decays exponentially. Whenever we have exponential decay, we get a correlation decay that sets the size of the boxes. It is very difficult to condition on quantum states, but recovery maps provide a sense in which it is meaningful to do so. The paper gives an efficient method of preparing Gibbs states and simulating quantum systems on quantum computers.

Additional reading

The standard treatment of information theory is [11]. This book contains definitions and properties of entropy, conditional entropy, mutual information, and conditional mutual information.

To see a treatment of the subject of Markov chains from the perspective of probability theory, see [12] or the mathematically more sophisticated counterpart [13]. An introduction to coupling can be found in [14], as well as [4] (the latter also contains an exposition to spatial mixing). The connection between Markov chain mixing and the so-called logarithmic Sobolev inequality is described in [15].

Appendix: Intuition for Markov chains

Random walk on the cycle

We have nn

n points on the cycle, 0,1,…,n−10,1,…,n−1

0,1,\dotsc,n-1. At each step, we move left or right with probability 1/21/2

1/2. We can write the transition matrix as

T=S+S−12T=S+S−12

T = \frac{S + S^{-1}}{2}

where SS

S is the shift operator S|x⟩=|x+1modn⟩S|x⟩=|x+1modn⟩

S\ket{x} = \ket{x+1 \bmod n}. The matrix SS

S is diagonalized by the Fourier transform. Define, for k=0,1,…,n−1k=0,1,…,n−1

k=0,1,\dotsc,n-1,

\ket{\tilde k} = \frac{1}{\sqrt n} \sum_{x=0}^{n-1} \exp\Bigl( \frac{2\pi i k x}{n} \Bigr) \ket{x}

We have the same amount of amplitude at every point, but there is a varying phase which depends on kk

k. If k=0k=0

k = 0, we get the all-ones vector. If kk

k is small, then the phase is slowly varying. If kk

k is large, then the phase is rapidly varying. Look at what happens after we apply the shift operator:

After the shift, we pick up an additional phase based on how rapidly the phase is varying. From this, we get:

The eigenvalues are

λk=cos2πkn,k=0,1,…,n−1.λk=cos⁡2πkn,k=0,1,…,n−1.

\lambda_k = \cos \frac{2\pi k}{n}, \qquad k=0,1,\dotsc,n-1.

Only k=0k=0

k = 0 will give me an eigenvalue of 11

How do we analyze Tt|p⟩Tt|p⟩

T^t \ket{p}? We should Fourier transform the distribution.

\begin{aligned} T^t \ket{p} = T^t \sum_{k=0}^{n-1} p_k \ket{\tilde{k}} = \sum_{k=0}^{n-1} p_k \lambda_k^t \ket{\tilde k}.\end{aligned}

If nn

n is odd, then as t→∞t→∞

t\rightarrow\infty, λtk→0λkt→0

\lambda_k^t \to 0 for all k=1,…,n−1k=1,…,n−1

k=1,\dotsc,n-1, so Tt→|π⟩⟨1n|Tt→|π⟩⟨1n|

T^t \to \ket{\pi}\bra{1_n}. Whatever you put into this operator, you get ππ

\pi out.

Spectral gap

The example of the random walk on the cycle shows that there is generally a unique stationary distribution and suggests that the speed of convergence is determined by how close the other eigenvalues are to 11

1. Specifically, suppose for simplicity that the eigenvalues of TT

T are 1=λ0≥λ1≥⋯≥01=λ0≥λ1≥⋯≥0

1 = \lambda_0 \ge \lambda_1\ge\cdots \ge 0 (real and positive). Then, the convergence time is on the order of ∼1/(1−λ1)∼1/(1−λ1)

\sim 1/(1-\lambda_1).

Typically, the distance of the eigenvalues from 11

1 reflects the size of the physical system. Even from the simple example, we can get some physical intuition from this. If kk

k is small, then the spectral gap is cos(2πk/n)=1−O(k2/n2)cos⁡(2πk/n)=1−O(k2/n2)

\cos(2\pi k/n) = 1-O(k^2/n^2). Thus, the convergence time is ∼1/(1−λ1)∼n2∼1/(1−λ1)∼n2

\sim 1/(1-\lambda_1) \sim n^2, which is indeed the convergence time for a random walk on a cycle.

References

S. Gharibian, Y. Huang, Z. Landau, and S. W. Shin, “Quantum Hamiltonian complexity,” Found. Trends Theor. Comput. Sci., vol. 10, no. 3, pp. front matter, 159–282, 2014.

R. W. Keener, Theoretical statistics. Springer, New York, 2010, p. xviii+538.

E. Crosson, D. Bacon, and K. R. Brown, “Making Classical Ground State Spin Computing Fault-Tolerant,” Physical Review E, vol. 82, no. 3, Sep. 2010.

C. Moore and S. Mertens, The nature of computation. Oxford University Press, Oxford, 2011, p. xviii+985.

F. Martinelli, “Lectures on Glauber dynamics for discrete spin models,” in Lectures on probability theory and statistics (Saint-Flour, 1997), vol. 1717, Springer, Berlin, 1999, pp. 93–191.

F. Martinelli and E. Olivieri, “Finite volume mixing conditions for lattice spin systems and exponential approach to equilibrium of Glauber dynamics,” in Cellular automata and cooperative systems (Les Houches, 1992), vol. 396, Kluwer Acad. Publ., Dordrecht, 1993, pp. 473–490.

M. S. Leifer and D. Poulin, “Quantum graphical models and belief propagation,” Ann. Physics, vol. 323, no. 8, pp. 1899–1946, 2008.

M. J. Kastoryano and F. G. S. L. Brandão, “Quantum Gibbs samplers: the commuting case,” Comm. Math. Phys., vol. 344, no. 3, pp. 915–957, 2016.

F. G. S. L. Brandão and M. J. Kastoryano, “Finite correlation length implies efficient preparation of quantum thermal states,” ArXiv e-prints, Sep. 2016.

O. Fawzi and R. Renner, “Quantum conditional mutual information and approximate Markov chains,” Comm. Math. Phys., vol. 340, no. 2, pp. 575–611, 2015.

T. M. Cover and J. A. Thomas, Elements of information theory, Second. Wiley-Interscience [John Wiley & Sons], Hoboken, NJ, 2006, p. xxiv+748.

R. Durrett, Essentials of stochastic processes. Springer, Cham, 2016, p. ix+275.

R. Durrett, Probability: theory and examples, Fourth., vol. 31. Cambridge University Press, Cambridge, 2010, p. x+428.

M. Mitzenmacher and E. Upfal, Probability and computing, Second. Cambridge University Press, Cambridge, 2017, p. xx+467.

F. Cesi, “Quasi-factorization of the entropy and logarithmic Sobolev inequalities for Gibbs random fields,” Probab. Theory Related Fields, vol. 120, no. 4, pp. 569–584, 2001.

This is the opposite of the probabilists’ convention, i.e., the transition probability matrix that we define here is the transpose of the one usually found in most probability theory textbooks. ↩

As a side note, it may be a good research question to investigate to what extent quantum algorithms can be used to compute summations whose terms are possibly negative. In quantum Monte Carlo, the quantum Hamiltonian is converted to a classical energy function; this conversion always works, but sometimes you end up with complex energies, which is terrible for estimating the partition function because terms can cancel each other out. ↩

You may recognize this as the total variation norm. ↩

Petz wrote about quantum relative entropy in 1991, way before it was cool. ↩

Windows On Theory published first on Windows On Theory

#Technology #Finance

0 notes