Solutions to Blum, Hopcroft & Kannan’s Foundations of Data Science

(1) \(E(x)=\frac12\), \(E(x^2)=\frac13\), \(E(x-y)=0\), \(E(xy)=\frac14\), \(E(x-y)^2=\frac16\). Indeed \(E(x^2)=\int_0^1 t^2\,dt=\frac13\), independence gives \(E(xy)=E(x)E(y)=\frac14\), and

\begin{equation*} E(x-y)^2=E(x^2)-2E(x)E(y)+E(y^2)=\tfrac13-\tfrac12+\tfrac13=\tfrac16 . \end{equation*}

(2) \(E(x)=0\), \(E(x^2)=\frac1{12}\), \(E(x-y)=0\), \(E(xy)=0\), \(E(x-y)^2=\frac16\):

\begin{equation*} E(x^2)=\int_{-1/2}^{1/2}t^2\,dt=\tfrac1{12},\qquad E(x-y)^2=\tfrac1{12}-0+\tfrac1{12}=\tfrac16 . \end{equation*}

(3) \(d/6\). The coordinates of the two points \(\mathbf x,\mathbf y\) are independent uniforms on \([0,1]\), so by linearity and part 1,

\begin{equation*} E|\mathbf x-\mathbf y|^2=\sum_{i=1}^d E(x_i-y_i)^2=\frac d6 . \end{equation*}

All \(\binom{30}{2}=435\) distances cluster tightly around \(\sqrt{100/6}\approx 4.08\) and all angles cluster around \(90^\circ\). The program (seed 0):

import numpy as np, matplotlib.pyplot as plt
from itertools import combinations
rng = np.random.default_rng(0)
X = rng.uniform(-0.5, 0.5, size=(30, 100))
dist, ang = [], []
for i, j in combinations(range(30), 2):
    dist.append(np.linalg.norm(X[i] - X[j]))
    c = X[i] @ X[j] / (np.linalg.norm(X[i]) * np.linalg.norm(X[j]))
    ang.append(np.degrees(np.arccos(c)))
fig, ax = plt.subplots(1, 2, figsize=(9, 3.5))
ax[0].hist(dist, bins=30); ax[0].set_xlabel('pairwise distance')
ax[1].hist(ang, bins=30); ax[1].set_xlabel('pairwise angle (degrees)')
plt.tight_layout()
for ext in ('svg', 'pdf'): plt.savefig('fds-ch02-cube-pairs.' + ext)

The distances have mean \(4.080\), standard deviation \(0.245\) and range \([3.26,\,4.73]\); the angles have mean \(89.93^\circ\), standard deviation \(5.91^\circ\) and range \([69.3^\circ,\,107.0^\circ]\); the mean norm is \(2.889\). This matches Exercise 2.1 (part 3, shifted cube) with \(d=100\):

\begin{equation*} \begin{aligned} E|\mathbf x-\mathbf y|^2&=\tfrac{d}{6},\quad &\sqrt{100/6}&\approx 4.082,\\ E|\mathbf x|^2&=\tfrac{d}{12},\quad &\sqrt{100/12}&\approx 2.887, \end{aligned} \end{equation*}

and \(|\mathbf x-\mathbf y|^2\) is a sum of \(d\) independent terms, so by the law of large numbers (Theorem 2.4) its relative fluctuation is \(O(1/\sqrt d)\): since \(\operatorname{Var}(x_i-y_i)^2=\frac1{15}-\frac1{36}=\frac7{180}\), the distance has standard deviation \(\approx\sqrt{100\cdot 7/180}\,/\,(2\cdot 4.082)\approx 0.24\). Likewise \(\cos\theta=\mathbf x\cdot\mathbf y/(|\mathbf x||\mathbf y|)\) has mean \(0\) and standard deviation \(\approx\sqrt{d}\,\frac1{12}/\frac{d}{12}=1/\sqrt d=0.1\), i.e. about \(5.7^\circ\) around \(90^\circ\).

Put mass \(\frac1a\) at \(a\) and the rest at \(0\).

(1)

(2) For \(a\ge1\) let \(p(0)=1-\frac1a\) and \(p(a)=\frac1a\); this is a probability distribution because \(0<\frac1a\le 1\), and \(x\ge0\). Then

\begin{equation*} E(x)=a\cdot\tfrac1a=1,\qquad \mathrm{Prob}(x\ge a)=\tfrac1a=\frac{E(x)}{a}, \end{equation*}

so equality holds in Markov’s inequality (Theorem 2.1).

Take \(p( c)=p(-c)=\frac1{2c^2}\) and \(p(0)=1-\frac1{c^2}\), a probability distribution because \(c\ge1\) gives \(0<\frac1{c^2}\le1\). By symmetry \(E(x)=0\), and

\begin{equation*} \mathrm{Var}(x)=E(x^2)=c^2\cdot\tfrac1{c^2}=1,\qquad \mathrm{Prob}(|x-E(x)|\ge c)=\mathrm{Prob}(x=\pm c)=\tfrac1{c^2}=\frac{\mathrm{Var}(x)}{c^2}, \end{equation*}

so Chebyshev’s inequality (Theorem 2.3) holds with equality.

The true value is \(\mathrm{Prob}(x\ge3)=\frac14\), and \(x\ge0\) so \(|x|=x\). All three parts use \(E(x^r)=\frac14\int_0^4 t^r\,dt=\frac{4^r}{r+1}\) and Markov’s inequality (Theorem 2.1) applied to the nonnegative variable \(x^r\).

(1) \(\mathrm{Prob}(x\ge3)\le \frac{E(x)}{3}=\frac23\).

(2) \(\mathrm{Prob}(x\ge3)=\mathrm{Prob}(x^2\ge9)\le\frac{E(x^2)}{9}=\frac{16/3}{9}=\frac{16}{27}\approx0.593\).

(3) For \(r>0\),

\begin{equation*} \mathrm{Prob}(x\ge3)=\mathrm{Prob}(x^r\ge3^r)\le\frac{E(x^r)}{3^r}=\frac{(4/3)^r}{r+1}. \end{equation*}

Setting the derivative of \(r\ln\frac43-\ln(r+1)\) to zero gives \(r+1=1/\ln\frac43\), i.e. \(r\approx2.476\), where the bound attains its minimum

\begin{equation*} \frac{(4/3)^r}{r+1}=\tfrac34\,e\,\ln\tfrac43\approx0.587 . \end{equation*}

Among integers, \(r=2\) and \(r=3\) tie at \(\frac{16}{27}\), while \(r=4\) gives \(\frac{256}{405}\approx0.632\) and the bound tends to \(\infty\) as \(r\to\infty\); so no moment bound gets below \(0.587\).

All three bounds equal the true probability \(\frac1a\) for every \(a\ge1\). Indeed \(E(x^k)=a^k\cdot\frac1a=a^{k-1}\), so Markov’s inequality (Theorem 2.1) applied to \(x^k\) gives

\begin{equation*} \mathrm{Prob}(x\ge a)=\mathrm{Prob}(x^k\ge a^k)\le\frac{E(x^k)}{a^k}=\frac{a^{k-1}}{a^k}=\frac1a, \qquad k=1,2,4, \end{equation*}

and \(\mathrm{Prob}(x\ge a)=p(x=a)=\frac1a\). The plot (for \(1\le a\le10\); numerically the curves differ by at most \(10^{-16}\)) therefore shows four coincident curves:

(1) \(c=3\), since \(\int_1^\infty c\,x^{-4}\,dx=\frac c3\) must equal \(1\).

(2) With seed 0 the sample mean was \(1.4779\), within \(0.022\) of \(E(x)=\frac32\). The relevant moments are

\begin{equation*} E(x)=\int_1^\infty \frac{3}{x^3}\,dx=\frac32,\qquad E(x^2)=\int_1^\infty \frac{3}{x^2}\,dx=3,\qquad \mathrm{Var}(x)=3-\tfrac94=\tfrac34 . \end{equation*}

Samples come from the inverse CDF: \(F(x)=1-x^{-3}\), so \(x=u^{-1/3}\) with \(u\) uniform on \((0,1)\).

import numpy as np
rng = np.random.default_rng(0)
x = rng.uniform(size=100) ** (-1/3)     # inverse CDF of 3/x^4 on [1, inf)
print(x.mean(), x.mean() - 1.5, x.std(ddof=1))

The output was mean \(1.4779\), error \(-0.0221\), sample standard deviation \(0.825\) (largest sample \(7.15\)). The error is of the expected size: the variance is finite, so by the law of large numbers (Theorem 2.4) the average of \(n=100\) samples has standard deviation

\begin{equation*} \sqrt{\mathrm{Var}(x)/n}=\sqrt{0.75/100}\approx0.087, \end{equation*}

and \(|{-0.022}|\) is about a quarter of that.

\(\mathbb{E}|\mathbf{x}|^2 = d/2\). For \(\mathbf{x}=(x_1,\dots,x_d)\sim G\) each \(x_i \sim N(0,\tfrac12)\), so by linearity

\begin{equation*} \begin{aligned} \mathbb{E}|\mathbf{x}|^2&=\sum_{i=1}^d \mathbb{E}\,x_i^2 =\sum_{i=1}^d \Big(\operatorname{Var}(x_i)+(\mathbb{E}x_i)^2\Big)\\ &=d\cdot\tfrac12=\frac d2 . \end{aligned} \end{equation*}

\(\operatorname{Var}(x_1)=1/d\). The uniform distribution on the sphere is invariant under the reflection \(x_1\mapsto -x_1\), so \(\mathbb{E}x_1=0\); it is also invariant under permutations of the coordinates, so \(\mathbb{E}x_1^2=\cdots=\mathbb{E}x_d^2\). Since \(|\mathbf{x}|=1\) always,

\begin{equation*} d\,\mathbb{E}x_1^2=\sum_{i=1}^d\mathbb{E}x_i^2=\mathbb{E}|\mathbf{x}|^2=1 \quad\Longrightarrow\quad \operatorname{Var}(x_1)=\mathbb{E}x_1^2-(\mathbb{E}x_1)^2=\frac1d . \end{equation*}

\(\varepsilon \ge 1-0.01^{1/1000}\approx 0.00459\). The shell is \(S\setminus(1-\varepsilon)S\), and by the scaling identity of Section 2.3, \(\operatorname{vol}((1-\varepsilon)S)=(1-\varepsilon)^{d}\operatorname{vol}(S)\). With \(d=1000\),

\begin{equation*} \begin{aligned} 1-(1-\varepsilon)^{1000}\ge 0.99 &\iff (1-\varepsilon)^{1000}\le 0.01\\ &\iff \varepsilon\ge 1-e^{-\ln 100/1000}=0.0045946\ldots \end{aligned} \end{equation*}

Let \(g(x)=e^x-1-x\). Then \(g^{\prime}(x)=e^x-1\) is negative for \(x<0\) and positive for \(x>0\), so \(g\) decreases on \((-\infty,0]\), increases on \([0,\infty)\), and

\begin{equation*} g(x)\ge g(0)=0\quad\text{for all }x,\qquad\text{i.e.}\quad 1+x\le e^x . \end{equation*}

Since \(g\ge 0\), the error is \(|e^x-(1+x)|=g(x)\), and by the monotonicity just shown \(\{x: g(x)\le 0.01\}\) is the interval between the two roots of \(e^x-1-x=0.01\). Solving numerically (bisection),

\begin{equation*} -0.1448\le x\le 0.1382 . \end{equation*}

\(d=5\), with \(V(5)=8\pi^2/15\approx 5.264\). By Lemma 2.6, \(V(d)=2\pi^{d/2}/(d\,\Gamma(\frac d2))=\pi^{d/2}/\Gamma(\frac d2+1)\), and \(\Gamma(x+1)=x\Gamma(x)\) gives

\begin{equation*} \frac{V(d)}{V(d-1)}=\frac{\sqrt\pi\,\Gamma(\frac{d+1}{2})}{\Gamma(\frac d2+1)},\qquad \frac{V(d)}{V(d-2)}=\frac{\pi\,\Gamma(\frac d2)}{\Gamma(\frac d2+1)}=\frac{2\pi}{d}. \end{equation*}

From \(V(1)=2\), \(V(2)=\pi\), \(V(3)=\frac43\pi\), \(V(4)=\frac{\pi^2}{2}\), \(V(5)=\frac{8\pi^2}{15}\), \(V(6)=\frac{\pi^3}{6}\) the successive ratios are

\begin{equation*} \frac{V(d)}{V(d-1)}=\frac\pi2,\ \frac43,\ \frac{3\pi}{8},\ \frac{16}{15},\ \frac{5\pi}{16} \qquad(d=2,\dots,6), \end{equation*}

so \(V\) increases up to \(d=5\) and \(V(6)<V(5)\). For \(d\ge 7\) we have \(2\pi/d<1\), so \(V(d)<V(d-2)\); iterating within each parity class,

\begin{equation*} V(d)<V(5)\ (d\ge 7\text{ odd}),\qquad V(d)<V(6)<V(5)\ (d\ge 8\text{ even}). \end{equation*}

Hence \(\max_d V(d)=V(5)\).

Take the cube \([0,1]^d\); a \(k\)-face is obtained by choosing \(k\) free coordinates and fixing each of the other \(d-k\) at \(0\) or \(1\).

(1) \(\binom dk 2^{d-k}\) (choose the free coordinates, then one of \(2^{d-k}\) fixings).

(2) By the binomial theorem,

\begin{equation*} \sum_{k=0}^d\binom dk 2^{d-k}1^k=(2+1)^d=3^d . \end{equation*}

(3) The surface consists of the \(\binom{d}{d-1}2^1=2d\) facets, each a unit \((d-1)\)-cube of area \(1\), so the surface area is \(2d\).

(4) Each of the \(2d\) facets is a \((d-1)\)-cube of side 2, of area \(2^{d-1}\), so the surface area is \(2d\cdot 2^{d-1}=d\,2^d\).

(5) For \(\mathbf{x}\in[0,1]^d\) the distance to the boundary is \(\min_i\min(x_i,1-x_i)\), so the points at distance at least \(\varepsilon\) from the surface form \([\varepsilon,1-\varepsilon]^d\). Using \(1+x\le e^x\) (Exercise 2.11),

\begin{equation*} \operatorname{vol}\big([\varepsilon,1-\varepsilon]^d\big)=(1-2\varepsilon)^d\le e^{-2\varepsilon d}. \end{equation*}

Hence at least a \(1-e^{-2\varepsilon d}\) fraction of the unit volume lies within distance \(\varepsilon\) of the surface; taking \(\varepsilon=c/d\), all but \(e^{-2c}\) of the volume lies within \(c/d\) of the surface.

\(dA=\sin\theta\,d\theta\,d\phi\), the area is \(2\pi(1-\cos\alpha)\), and for \(\alpha=36^\circ\) it equals \(\pi(3-\sqrt5)/2\approx 1.200\). Put the cone’s axis along the \(z\)-axis and use spherical coordinates with polar angle \(\theta\) (from the axis) and azimuth \(\phi\); on the unit sphere a \(d\theta\) step has length \(d\theta\) and a \(d\phi\) step has length \(\sin\theta\,d\phi\), giving the area element above. The cone of half-angle \(\alpha\) is \(\{\theta\le\alpha\}\), so

\begin{equation*} \begin{aligned} A(\alpha)&=\int_0^{2\pi}\!\!\int_0^{\alpha}\sin\theta\,d\theta\,d\phi=2\pi(1-\cos\alpha),\\ A(36^\circ)&=2\pi\Big(1-\tfrac{1+\sqrt5}{4}\Big)=\frac{\pi(3-\sqrt5)}{2}\approx 1.19998, \end{aligned} \end{equation*}

using \(\cos 36^\circ=\frac{1+\sqrt5}{4}\).

Flat top: \(V(3) = \tfrac43\pi\); hemispherical top: \(\tfrac12 A(4) = \pi^2\).

(i) The flat top is the slice \(\{x_1 = 1,\ x_2^2+x_3^2+x_4^2 \le 1\}\), a solid unit 3-ball, so its 3-dimensional measure is \(V(3) = \tfrac43\pi\).

(ii) The rounded top is the boundary of the half 4-ball \(\{x_1 \ge 1,\ (x_1-1)^2+x_2^2+x_3^2+x_4^2 \le 1\}\) minus its flat base, i.e. half of the unit 3-sphere in \(\mathbb{R}^4\); by Lemma 2.6

\begin{equation*} \tfrac12 A(4) = \frac12\cdot\frac{2\pi^{2}}{\Gamma(2)} = \pi^2 . \end{equation*}

(1)

\begin{equation*} \text{area} = h\, r^{d-2} A(d-1) + 2\, r^{d-1} V(d-1). \end{equation*}

(2)

\begin{equation*} \text{volume} = h\, r^{d-1} V(d-1). \end{equation*}

For the cylinder \(\{x : 0 \le x_1 \le h,\ x_2^2+\cdots+x_d^2 \le r^2\}\) the cross-section is a \((d-1)\)-ball of radius \(r\), with volume \(r^{d-1}V(d-1)\) and surface \(r^{d-2}A(d-1)\); the lateral surface is \(h\) times the latter, the two end caps are cross-sections, and the volume is \(h\) times the cross-section.

Every constant radius gives volume tending to \(0\); constant volume needs \(r \approx \sqrt{d/(2\pi e)}\).

A ball of radius \(r\) has volume \(r^d V(d)\), and by Lemma 2.6, \(V(d) = \pi^{d/2}/\Gamma(\tfrac d2+1)\). Stirling with \(n = d/2\) gives

\begin{equation*} r^d V(d) \approx r^d \pi^{d/2}\left(\frac{2e}{d}\right)^{d/2} = \left(\frac{2\pi e\, r^2}{d}\right)^{d/2}. \end{equation*}

• Radius \(2\): the volume is \(\approx (8\pi e/d)^{d/2}\). Exactly, \(r^dV(d)/\big(r^{d-2}V(d-2)\big) = 2\pi r^2/d = 8\pi/d\), which exceeds \(1\) exactly when \(d < 8\pi \approx 25.1\); the volume therefore rises, peaks at \(d = 24\) (volume \(\approx 3.24\times 10^4\)), falls below \(1\) at \(d = 63\), and then tends to \(0\) faster than any exponential (e.g. \(\approx 3.0\times10^{-10}\) at \(d=100\)).
• Constant radius \(r > 2\): the same picture with the peak moved to \(d \approx 2\pi r^2\); the base \(2\pi e r^2/d\) drops below \(1\) once \(d > 2\pi e r^2\), so the volume still tends to \(0\).
• Constant volume: the base must stay at \(1\), i.e.

\begin{equation*} r = \sqrt{\frac{d}{2\pi e}} \approx 0.242\sqrt d . \end{equation*}

There is no contradiction: \(V(d)\) and \(V(3)\) are measured in different units, so they cannot be compared.

\(V(3)\) is a 3-dimensional measure (in units of length cubed) and \(V(d)\) is a \(d\)-dimensional measure (in units of length to the \(d\)). The containment does give an inequality, but between \(d\)-dimensional measures: the copy \(\{x : x_4=\cdots=x_d=0,\ x_1^2+x_2^2+x_3^2\le 1\}\) of the 3-ball lies in a 3-dimensional subspace, so its \(d\)-dimensional volume is \(0 \le V(d)\). The numbers \(V(d) < V(3)\) only mean that the \(d\)-ball holds fewer unit \(d\)-cubes than the 3-ball holds unit 3-cubes. The two counts are unrelated: the unit \(d\)-cube \([-\tfrac12,\tfrac12]^d\) has volume \(1\) but its corners are at distance \(\sqrt d/2\), so for large \(d\) it sticks far outside the unit ball.

(1)

\begin{equation*} V(d) = V(d-1)\int_{-1}^{1} (1-t^2)^{\frac{d-1}{2}}\, dt . \end{equation*}

The slice \(\{x_1 = t\}\) of the unit ball is a \((d-1)\)-ball of radius \(\sqrt{1-t^2}\), whose volume is \((\sqrt{1-t^2})^{d-1}V(d-1)\); integrating over \(t \in [-1,1]\) gives the display.

(2) With \(V(1) = 2\) (the interval \([-1,1]\)):

\begin{equation*} \begin{aligned} V(2) &= 2\int_{-1}^{1}\sqrt{1-t^2}\,dt = 2\cdot\frac{\pi}{2} = \pi,\\ V(3) &= \pi\int_{-1}^{1}(1-t^2)\,dt = \pi\left(2-\frac23\right) = \frac43\pi, \end{aligned} \end{equation*}

where \(\int_{-1}^1\sqrt{1-t^2}\,dt = \int_{-\pi/2}^{\pi/2}\cos^2\theta\,d\theta = \pi/2\) by \(t = \sin\theta\).

Take \(c \ge 2\) with \(\frac{8}{c}e^{-c^2/8} < \epsilon\); then \(s \ge c/\sqrt{d-1}\) gives (for \(d = 2\), \(s \ge c \ge 2\) makes \(A\cap B\) null, so take \(d \ge 3\)) \(\mathrm{Prob}(x\in A\cap B) \le \frac8c e^{-c^2/8} < \epsilon\).

Let \(u\) be the unit vector with \(B\) centered at \(su\), and \(H = \{x : u\cdot x \ge s/2\}\). Each component is uniform on a ball of volume \(V(d)\), so

\begin{equation*} \mathrm{Prob}(x \in A\cap B) = \tfrac12\frac{\mathrm{vol}(A\cap B)}{V(d)} + \tfrac12\frac{\mathrm{vol}(A\cap B)}{V(d)} = \frac{\mathrm{vol}(A\cap B)}{V(d)} . \end{equation*}

The reflection \(R(x) = x + (s - 2u\cdot x)u\) in the hyperplane \(u\cdot x = s/2\) swaps \(0\) and \(su\), hence maps \(A\cap B\) onto itself and \((A\cap B)\setminus H\) into \((A\cap B)\cap H\). Therefore

\begin{equation*} \begin{aligned} \mathrm{vol}(A\cap B) &\le 2\,\mathrm{vol}(A \cap H) \le 2\,\mathrm{vol}\big(A\cap\{|u\cdot x| \ge s/2\}\big)\\ &\le 2\cdot\frac{2}{c^{\prime}}e^{-c^{\prime\,2}/2}\,V(d), \qquad c^{\prime} = \frac{s\sqrt{d-1}}{2}, \end{aligned} \end{equation*}

by Theorem 2.7 applied in direction \(u\) (rotational symmetry of the ball), which needs \(c^{\prime} \ge 1\) and \(d \ge 3\). If \(s \ge c/\sqrt{d-1}\) then \(c^{\prime} \ge c/2 \ge 1\), and since \(\frac{4}{c^{\prime}}e^{-c^{\prime\,2}/2}\) is decreasing in \(c^{\prime}\),

\begin{equation*} \mathrm{Prob}(x\in A\cap B) \le \frac{4}{c/2}e^{-c^2/8} = \frac8c e^{-c^2/8} < \epsilon . \end{equation*}

The height is \(h = 1/\sqrt d\) (with radius \(r = \sqrt{(d-1)/d}\)).

As the statement intends, the cylinder’s axis is the north direction \(x_d\); let its base be a \((d-1)\)-ball of radius \(r\) centred at \((c,a)\) with \(a \ge 0\), and let its height be \(h\). Its farthest point from the origin is on the top rim, so it lies in the hemisphere iff \((|c|+r)^2 + (a+h)^2 \le 1\); taking \(c = 0\), \(a = 0\) loosens this, so the optimal cylinder stands on the flat face, centred, with \(r^2 + h^2 = 1\). By Exercise 2.16 its volume is

\begin{equation*} h\,r^{d-1}V(d-1) = V(d-1)\, h\,(1-h^2)^{\frac{d-1}{2}}, \qquad 0 < h < 1 . \end{equation*}

Maximise \(g(h) = \ln h + \frac{d-1}{2}\ln(1-h^2)\), which is strictly concave:

\begin{equation*} g^{\prime}(h) = \frac1h - \frac{(d-1)h}{1-h^2} = 0 \iff 1 - h^2 = (d-1)h^2 \iff h = \frac{1}{\sqrt d}, \end{equation*}

and then \(r = \sqrt{1 - 1/d}\).

The volume is \(\left(2/\sqrt{d}\right)^{d} = (4/d)^{d/2}\).

A hypercube of side \(s\) has two opposite vertices \(u, w\) with \(|u-w| = s\sqrt{d}\) (its main diagonal). If the cube lies in the unit ball then both vertices do, so

\begin{equation*} s\sqrt{d} = |u-w| \le |u| + |w| \le 2, \qquad\text{i.e.}\qquad s \le \frac{2}{\sqrt{d}} . \end{equation*}

Conversely the cube \([-1/\sqrt{d},\,1/\sqrt{d}]^d\) has side \(2/\sqrt{d}\) and every point satisfies \(|x|^2 \le d\cdot \frac1d = 1\), so it fits. Hence the largest cube has side \(2/\sqrt d\) and volume

\begin{equation*} \left(\frac{2}{\sqrt d}\right)^{d} = \left(\frac{4}{d}\right)^{d/2}. \end{equation*}

The ratio is

\begin{equation*} R_d(\epsilon) = \frac{\int_\epsilon^1 (1-t^2)^{\frac{d-1}{2}}\,dt}{\int_0^1 (1-t^2)^{\frac{d-1}{2}}\,dt}, \end{equation*}

with the values in the table below (rows “vol” read “area” as the \(d\)-dimensional volume, as in the proof of Theorem 2.7; rows “surf” read it as surface area, which replaces the exponent by \(\frac{d-3}{2}\)).

As in the proof of Theorem 2.7, the slice of the ball at \(x_1 = t\) is a \((d-1)\)-ball of radius \(\sqrt{1-t^2}\), of volume \((1-t^2)^{\frac{d-1}{2}}V(d-1)\); the factor \(V(d-1)\) cancels in the ratio. Substituting \(u = t^2\) turns the numerator and denominator into incomplete and complete Beta integrals, so \(R_d(\epsilon) = 1 - I_{\epsilon^2}\!\left(\tfrac12, \tfrac{d+1}{2}\right)\) (surface: \(1 - I_{\epsilon^2}(\tfrac12, \tfrac{d-1}{2})\)) with \(I\) the regularized incomplete Beta function, evaluated numerically:

There is no conflict: the equatorial slabs for different poles overlap, and almost all of the volume lies in all of them at once.

Fix poles (unit vectors) \(v_1,\dots,v_k\) and let \(S_j = \{x : |x\cdot v_j| \le c/\sqrt{d-1}\}\) be the slab about the equator of \(v_j\). By rotational symmetry Theorem 2.7 applies to each \(v_j\) (it was proved for \(v = e_1\) only for convenience), so each \(S_j\) misses at most a \(\frac2c e^{-c^2/2}\) fraction of the volume, and by the union bound

\begin{equation*} \begin{aligned} \frac{\mathrm{vol}\big(B \setminus (S_1\cap\cdots\cap S_k)\big)}{\mathrm{vol}(B)} &\le \sum_{j=1}^k \frac{\mathrm{vol}(B\setminus S_j)}{\mathrm{vol}(B)}\\ &\le \frac{2k}{c}\, e^{-c^2/2}. \end{aligned} \end{equation*}

Taking \(c = \sqrt{2\ln(4k)}\) (so \(c \ge 1\)) makes the right side at most \(\frac12\) for every \(d \ge 3\), so more than half the ball lies within \(\sqrt{2\ln (4k)}/\sqrt{d-1}\) of every one of the \(k\) equators simultaneously; for any fixed \(k\) this width tends to \(0\) as \(d\to\infty\). This is exactly the box argument after Theorem 2.8 with the \(d\) coordinate directions as poles.

The low-dimensional picture misleads because there a point near one equator is far from another. In high dimensions a typical point \(x\) of the ball has all its coordinates of size \(O(1/\sqrt d)\) in any orthonormal frame (the Discussion after Theorem 2.8), so \(x\cdot v = O(1/\sqrt d)\) for each fixed direction \(v\): a typical point is nearly orthogonal to every fixed direction, and so it sits near every equator.

The volume lives in the intersection of the slab and the annulus, a thin belt around the equator just under the surface, and that intersection is large because the slab is thin in only one of the \(d\) directions.

Let \(S = \{|x_1| \le c/\sqrt{d-1}\}\) and \(A = \{1-\epsilon \le |x| \le 1\}\). By Section 2.3 the complement of \(A\) holds at most an \(e^{-\epsilon d}\) fraction of the volume, and by Theorem 2.7 (with \(c\ge1\), \(d\ge3\)) the complement of \(S\) holds at most \(\frac2c e^{-c^2/2}\), so

\begin{equation*} \frac{\operatorname{vol}(S\cap A)}{\operatorname{vol}(B)} \;\ge\; 1 - e^{-\epsilon d} - \frac{2}{c}\,e^{-c^2/2}, \end{equation*}

which tends to \(1\) for, e.g., \(\epsilon = \frac{\ln d}{d}\) and \(c = \sqrt{2\ln d}\): both the slab width \(\sqrt{2\ln d}/\sqrt{d-1}\) and the annulus width \(\ln d/d\) go to \(0\).

The two constraints barely interact. A point of \(S\) has \(x_1^2 \le c^2/(d-1) = O\!\left(\tfrac{\ln d}{d}\right)\) in

\begin{equation*} |x|^2 = x_1^2 + (x_2^2 + \cdots + x_d^2), \end{equation*}

so being near the surface is a condition on the other \(d-1\) coordinates, which are free: the slab contains the whole \((d-1)\)-ball \(\{x_1 = 0\}\), whose volume is itself concentrated near its boundary sphere by Section 2.3. In the picture of the Discussion after Theorem 2.8, a typical point looks like \((\pm\frac1{\sqrt d},\dots,\pm\frac1{\sqrt d})\): norm about \(1\) and first coordinate \(O(1/\sqrt d)\) at once.

With seed 0 the five bands contain \(326, 336, 334, 348, 334\) of the 500 points, \(53\) points lie in all five, and the bands must have width \(1.002 \approx 7.08/\sqrt{50}\) to contain every point.

Points are generated as in Section 2.5 (normalized Gaussian vectors); point \(x\) is in the band of pole \(p\) when \(|x\cdot p| \le 1/\sqrt{50}\):

import numpy as np
rng = np.random.default_rng(0)
d, n = 50, 500
def sphere(m):
    g = rng.standard_normal((m, d))
    return g / np.linalg.norm(g, axis=1, keepdims=True)
X, P = sphere(n), sphere(5)
D = np.abs(X @ P.T)                 # |x . pole|, shape (500, 5)
inband = D <= 1 / np.sqrt(d)        # band of width 2/sqrt(d)
print("per band:", inband.sum(axis=0))
print("in all five:", inband.all(axis=1).sum())
print("width for all:", 2 * D.max())

For a uniform point \(x\) and fixed unit \(p\), \((x\cdot p)^2 \sim \mathrm{Beta}\big(\tfrac12,\tfrac{d-1}{2}\big)\) (the surface version of the slice computation in Theorem 2.7), so

\begin{equation*} \Pr\big(|x\cdot p| \le \tfrac{1}{\sqrt{50}}\big) = I_{1/50}\big(\tfrac12, \tfrac{49}{2}\big) = 0.678, \end{equation*}

predicting \(339\) points per band; given \(x\) the five independent poles act independently, so \(500\cdot 0.678^5 = 71.5\) points are expected in all five (the run’s \(53\) is a low draw). To capture all \(2500\) pairs the half-width must exceed \(\max|x\cdot p|\), which by the union-bound reasoning of Theorem 2.8 is about \(\sqrt{2\ln 2500}/\sqrt{d} \approx 0.56\); the run needed width \(2\max|x\cdot p| = 1.002\), about seven times \(2/\sqrt{50}\).

Each hyperplane can be brought to distance \(t \approx 3.89/\sqrt d\) from the origin (gap \(\approx 7.8/\sqrt d\)); Theorem 2.7 certifies the slightly larger \(t = 4.13/\sqrt{d-1}\).

Let \(v\) be the unit normal and the hyperplanes \(x\cdot v = \pm t\). By rotational symmetry \(x_i\cdot v\) has the law of the first coordinate \(x_1\) of a uniform point of the ball, whatever \(v\) is, so with \(q(t) = \Pr(|x_1| > t)\) the union bound gives

\begin{equation*} \Pr\big(\exists i : |x_i\cdot v| > t\big) \le 100\, q(t). \end{equation*}

Guaranteed bound. Theorem 2.7 (with \(c\ge 1\), \(d\ge3\)) gives \(q\big(c/\sqrt{d-1}\big) \le \frac2c e^{-c^2/2}\), so it suffices that

\begin{equation*} 100\cdot\frac{2}{c}\,e^{-c^2/2} \le 0.01 \iff \frac{2}{c}\,e^{-c^2/2} \le 10^{-4} \Longleftarrow c \ge 4.13 , \end{equation*}

i.e. \(t = 4.13/\sqrt{d-1}\) works for every \(d \ge 3\).

Sharp value. Given \(v\), the 100 points are independent, so the success probability is exactly \((1-q(t))^{100}\), and

\begin{equation*} (1-q(t))^{100} \ge 0.99 \iff q(t) \le 1 - 0.99^{1/100} = 1.005\times 10^{-4}. \end{equation*}

For large \(d\), \(\sqrt d\,x_1 \to N(0,1)\) (Section 2.5: \(x_1 = \rho\, g_1/|g|\) with \(|g| \approx \sqrt d\) by Theorem 2.9 and \(\rho \to 1\)), so \(q(t) \approx 2\big(1-\Phi(t\sqrt d)\big)\) and the condition becomes \(t\sqrt d \ge \Phi^{-1}(1 - 5.02\times10^{-5}) = 3.89\). The exact Beta law \(x_1^2 \sim \mathrm{Beta}\big(\tfrac12,\tfrac{d+1}{2}\big)\) gives \(t\sqrt d = 3.74,\ 3.87,\ 3.89\) for \(d = 10^2, 10^3, 10^4\).

The dot product is only \(O(\sqrt d)\) while \(|\mathbf x||\mathbf y| \approx d\), so \(\cos\theta = O(1/\sqrt d)\) with high probability (taking \(\mathbf x, \mathbf y\) independent, as intended).

Since the \(2d\) coordinates are independent \(N(0,1)\), the cross terms vanish and

\begin{equation*} \begin{aligned} E(\mathbf x\cdot\mathbf y) &= \sum_{i=1}^d E(x_i)E(y_i) = 0,\\ \operatorname{Var}(\mathbf x\cdot\mathbf y) &= \sum_{i,j} E(x_ix_j)E(y_iy_j) = \sum_{i=1}^d E(x_i^2)E(y_i^2) = d . \end{aligned} \end{equation*}

By Chebyshev’s inequality (Theorem 2.3), for every \(k>0\),

\begin{equation*} \Pr\big(|\mathbf x\cdot\mathbf y| \ge k\sqrt d\big) \le \frac{1}{k^2}. \end{equation*}

By the Gaussian Annulus Theorem (Theorem 2.9) with \(\beta = \sqrt d/2 \le \sqrt d\), each of \(|\mathbf x|, |\mathbf y|\) is at least \(\sqrt d/2\) except with probability \(3e^{-cd/4}\). On the intersection of the three good events, which has probability at least \(1 - \frac1{k^2} - 6e^{-cd/4}\),

\begin{equation*} |\cos\theta| = \frac{|\mathbf x\cdot\mathbf y|}{|\mathbf x|\,|\mathbf y|} \le \frac{k\sqrt d}{d/4} = \frac{4k}{\sqrt d}. \end{equation*}

Taking \(k = d^{1/4}\) gives \(|\cos\theta| \le 4d^{-1/4}\) with probability at least \(1 - d^{-1/2} - 6e^{-cd/4}\), so as \(d\to\infty\) the bound tends to \(0\) and the probability to \(1\): \(\mathbf x\) and \(\mathbf y\) are almost orthogonal, with \(\theta = \frac\pi2 \pm O(k/\sqrt d)\).

With probability \(1-O(1/d)\) the angle \(\theta\) satisfies \(\cos\theta\le \sqrt{6\ln d/(d-1)}\,/(1-2\ln d/d)^2\), which is below \(1/\sqrt2\) for every \(d\ge 85\), so \(\theta\ge 45^\circ\).

Take the two vectors uniform in the unit ball (the angle depends only on their directions, which are uniform on the sphere also for Gaussian or sphere-uniform vectors). Regard them as \(\mathbf x_1,\mathbf x_2\) among \(n=d\) independent uniform points \(\mathbf x_1,\dots,\mathbf x_d\) of the ball. By Theorem 2.8 with \(n=d\), with probability \(1-O(1/d)\) we have \(|\mathbf x_i|\ge 1-\frac{2\ln d}{d}\) for \(i=1,2\) and \(|\mathbf x_1\cdot\mathbf x_2|\le \frac{\sqrt{6\ln d}}{\sqrt{d-1}}\). On this event

\begin{equation*} \begin{aligned} \cos\theta\le\frac{|\mathbf x_1\cdot\mathbf x_2|}{|\mathbf x_1|\,|\mathbf x_2|} &\le \frac{\sqrt{6\ln d/(d-1)}}{\left(1-2\ln d/d\right)^2} \xrightarrow[d\to\infty]{}0 . \end{aligned} \end{equation*}

The right side is at most \(1/\sqrt2=\cos 45^\circ\) for all \(d\ge 85\) (Check!), hence

\begin{equation*} \operatorname{Prob}\left(\theta\ge 45^\circ\right)\ge 1-O(1/d)\longrightarrow 1 . \end{equation*}

The projected volume has density proportional to \((1-t^2/d)^{(d-1)/2}\), which tends to \(e^{-t^2/2}\): the projection is asymptotically the standard Gaussian \(N(0,1)\).

Take the line to be the \(x_1\)-axis. As in the proof of Theorem 2.7, the slice \(x_1=t\) of the ball of radius \(\sqrt d\) is a \((d-1)\)-ball of radius \(\sqrt{d-t^2}\), so the projected volume at \(t\in[-\sqrt d,\sqrt d]\) is

\begin{equation*} \begin{aligned} V(d-1)\,(d-t^2)^{\frac{d-1}{2}} &= V(d-1)\,d^{\frac{d-1}{2}}\left(1-\frac{t^2}{d}\right)^{\frac{d-1}{2}},\\ \left(1-\frac{t^2}{d}\right)^{\frac{d-1}{2}} &=\exp\left(\frac{d-1}{2}\ln\left(1-\frac{t^2}{d}\right)\right) \longrightarrow e^{-t^2/2}. \end{aligned} \end{equation*}

Since \(\ln(1-u)\le -u\), the integrand is dominated by \(e^{-t^2(d-1)/(2d)}\le e^{-t^2/4}\) for \(d\ge2\), so its integral tends to \(\int e^{-t^2/2}\,dt=\sqrt{2\pi}\) and the normalized projected density tends to \(\frac{1}{\sqrt{2\pi}}e^{-t^2/2}\).

Intuition: almost all the volume of the ball of radius \(\sqrt d\) lies in a thin shell at radius about \(\sqrt d\) (Section 2.3), which is exactly where a unit-variance spherical Gaussian puts almost all its mass (Theorem 2.9). Both distributions are rotation invariant, so uniform-in-the-ball and \(N(0,I_d)\) are nearly the same distribution, and the projection of \(N(0,I_d)\) on any unit direction is exactly \(N(0,1)\).

(1) Normalize spherical Gaussian vectors, as in Section 2.5: \(\mathbf x=\mathbf g/|\mathbf g|\) with \(\mathbf g\sim N(0,I_d)\) is uniform on the unit sphere because the Gaussian density is rotation invariant.

import numpy as np, matplotlib.pyplot as plt
rng = np.random.default_rng(2031)
def sphere_points(n, d):
    g = rng.standard_normal((n, d))
    return g / np.linalg.norm(g, axis=1, keepdims=True)
X = sphere_points(200, 50)                      # part 2
fig, axes = plt.subplots(1, 4, figsize=(11, 2.8), sharey=True)
for k, ax in enumerate(axes):                   # part 3
    u = sphere_points(1, 50)[0]                 # random line through 0
    p = X @ u                                   # projections onto the line
    print(k, p.mean(), p.std(), abs(p).max())
    ax.hist(p, bins=np.linspace(-1, 1, 41))
    ax.set_xlim(-1, 1); ax.set_title(f'line {k+1}')
plt.show()

(2) The call sphere_points(200, 50) produces the \(200\times 50\) matrix X of points. (3) For four random lines the projections had standard deviations \(0.149,\ 0.133,\ 0.144,\ 0.135\) (against \(1/\sqrt{50}=0.141\)), means within \(0.016\) of \(0\), and largest absolute value between \(0.40\) and \(0.48\). Between \(92.5\%\) and \(98\%\) of the points fell within \(2/\sqrt{50}=0.283\) of the origin. All four histograms are bell-shaped around \(0\) and leave the ends \(\pm1\) of the line empty.

(4) Relative to every line, the surface area is concentrated near the equator perpendicular to it: the projection onto a unit direction is approximately \(N(0,1/d)\), so almost all of the area lies within \(O(1/\sqrt d)\) of the hyperplane through the origin orthogonal to the line, and the polar caps near \(\pm1\) carry almost none (Theorem 2.7 and Section 2.4; the surface analogue holds by the same calculation).

(1) The projection is exactly \(N(0,1)\), whatever the line. For a unit vector \(\mathbf u\) along the line and \(\mathbf x\) with i.i.d. \(N(0,1)\) coordinates,

\begin{equation*} \mathbf u\cdot\mathbf x=\sum_{i=1}^d u_ix_i \sim N\left(0,\sum_{i=1}^d u_i^2\right)=N(0,1), \end{equation*}

a linear combination of independent Gaussians (this is the computation in the proof of Theorem 2.10). If \(\mathbf u\) is itself random and independent of \(\mathbf x\), conditioning on \(\mathbf u\) gives \(N(0,1)\) for every value of \(\mathbf u\), so the projection is still \(N(0,1)\). (2) Near the sphere of radius \(2\sqrt d\). Writing \(\mathbf x=2\mathbf z\) with \(\mathbf z\sim N(0,I_d)\), we get \(E|\mathbf x|^2=4d\), and by Theorem 2.9 applied to \(\mathbf z\), for \(\beta\le\sqrt d\),

\begin{equation*} \operatorname{Prob}\left(\bigl||\mathbf x|-2\sqrt d\bigr|\le 2\beta\right) \ge 1-3e^{-c\beta^2}, \end{equation*}

so the points lie in an annulus of width \(O(1)\) around radius \(2\sqrt d\).

In 3 dimensions the \(4950\) pairwise distances spread over all of \([0,2]\) with density rising linearly; in 100 dimensions they cluster tightly at \(\sqrt2\).

import numpy as np, matplotlib.pyplot as plt
from scipy.spatial.distance import pdist
rng = np.random.default_rng(2033)
fig, axes = plt.subplots(1, 2, figsize=(9, 3), sharex=True)
for ax, d in zip(axes, (3, 100)):
    g = rng.standard_normal((100, d))
    X = g / np.linalg.norm(g, axis=1, keepdims=True)   # uniform on sphere
    D = pdist(X)                                        # all 4950 distances
    print(d, D.mean(), D.std(), D.min(), D.max())
    ax.hist(D, bins=np.linspace(0, 2, 51), density=True)
    ax.set_title(f'd = {d}')
plt.show()

The run gave:

Both match the exact picture. For unit vectors \(|\mathbf x-\mathbf y|^2=2-2\,\mathbf x\cdot\mathbf y\). In \(d=3\), \(\mathbf x\cdot\mathbf y\) is uniform on \([-1,1]\) (Archimedes), so \(r=|\mathbf x-\mathbf y|\) has density \(r/2\) on \([0,2]\), with mean \(4/3\) and standard deviation \(\sqrt{2-16/9}=0.471\) (the dashed line in the left panel). In \(d=100\), \(\mathbf x\cdot\mathbf y\) is the projection of \(\mathbf y\) on the direction \(\mathbf x\), approximately \(N(0,1/d)\) (Exercises 2.31 and 2.32), so

\begin{equation*} r=\sqrt{2-2\,\mathbf x\cdot\mathbf y} \approx\sqrt2-\frac{\mathbf x\cdot\mathbf y}{\sqrt2}, \qquad \operatorname{sd}( r)\approx\frac{1}{\sqrt{2d}}=0.0707 . \end{equation*}

Yes: a new random vertex is close to perpendicular to essentially all ten poles, on average \(9.6\) of the ten within \(10^\circ\) of \(90^\circ\) and \(9.97\) of the ten within \(15^\circ\).

For a fixed pole \(\mathbf p\in\{\pm1\}^d\) and an independent uniform \(\mathbf q\in\{\pm1\}^d\), the products \(p_iq_i\) are independent uniform \(\pm1\) variables, whatever \(\mathbf p\) is, so

\begin{equation*} \cos\theta=\frac{\mathbf p\cdot\mathbf q}{d}=\frac1d\sum_{i=1}^d p_iq_i, \qquad E\cos\theta=0,\qquad \operatorname{Var}\cos\theta=\frac1d , \end{equation*}

and by Chebyshev (Theorem 2.3) \(\operatorname{Prob}(|\cos\theta|>t/\sqrt d)\le 1/t^2\); exactly, since \((\mathbf p\cdot\mathbf q+d)/2\sim\mathrm{Bin}(d,1/2)\), for \(d=128\) the probabilities of \(|\cos\theta|\le 2/\sqrt d=0.177\) and \(|\cos\theta|\le 3/\sqrt d=0.265\) are \(0.958\) and \(0.997\). As for the ball, the claim holds for every choice of pole.

Experiment (seed 2034): ten random poles, \(1000\) new random \(\pm1\) vectors, \(d=128\). The \(10\,000\) cosines had standard deviation \(0.0879\) (theory \(1/\sqrt{128}=0.0884\)) and largest absolute value \(0.359\), so every angle lay between \(68.9^\circ\) and \(111.1^\circ\). Counting, for each new vector, the poles it is close to perpendicular to:

The averages agree with \(10\times0.958=9.58\) and \(10\times0.997=9.97\).

Yes to all three: in each case all but a \(O(1/t^2)\) fraction lies within distance \(t\) of the equator, while the cube reaches distance \(\sqrt d/2\) from it (at the vertices \(\mathbf 0\) and \(\mathbf 1\)).

The cube is \([0,1]^d\) and the distance from \(\mathbf x\) to the equator is \(|S-d/2|/\sqrt d\) with \(S=\sum_i x_i\). Pick a point uniformly from the relevant set; by Chebyshev (Theorem 2.3),

\begin{equation*} \operatorname{Prob}\left(\frac{|S-d/2|}{\sqrt d}\ge t\right) \le\frac{\operatorname{Var}S}{t^2d} \qquad\text{when } ES=\frac d2 . \end{equation*}

(1) A uniform vertex has independent coordinates \(x_i\in\{0,1\}\) with \(\operatorname{Var}x_i=\frac14\), so \(ES=\frac d2\), \(\operatorname{Var}S=\frac d4\), and all but a \(\frac{1}{4t^2}\) fraction of the \(2^d\) vertices lie within distance \(t\) of the equator. (2) A uniform point of the cube has independent \(x_i\sim U[0,1]\) with \(\operatorname{Var}x_i=\frac1{12}\), so \(\operatorname{Var}S=\frac d{12}\) and all but a \(\frac{1}{12t^2}\) fraction of the volume lies within distance \(t\) of the equator. (3) The surface is the union of the \(2d\) facets \(\{x_i=b\}\), \(b\in\{0,1\}\), each of area \(1\), so it suffices to treat one facet. On it \(S=b+\sum_{j\ne i}x_j\) with the \(d-1\) free coordinates uniform, so \(ES=\frac d2\pm\frac12\), \(\operatorname{Var}S=\frac{d-1}{12}\), and for \(t>\frac{1}{2\sqrt d}\)

\begin{equation*} \begin{aligned} \operatorname{Prob}\left(\frac{|S-d/2|}{\sqrt d}\ge t\right) &\le\operatorname{Prob}\left(|S-ES|\ge t\sqrt d-\tfrac12\right)\\ &\le\frac{(d-1)/12}{\left(t\sqrt d-\frac12\right)^2} \le\frac{1}{12\left(t-\frac1{2\sqrt d}\right)^2}. \end{aligned} \end{equation*}

Averaging over the facets, all but this fraction of the surface area lies within distance \(t\) of the equator.

(1) If \(\mathbf{v}=\sum_i a_i\mathbf{e_i}=\sum_i b_i\mathbf{e_i}\) then \(\sum_i (a_i-b_i)\mathbf{e_i}=\mathbf{0}\), so linear independence forces \(a_i=b_i\) for every \(i\).

(2) \(|\mathbf{z}|^2=\frac12\):

\begin{equation*} \begin{aligned} \mathbf{z}&=\tfrac{\sqrt2}{2}(1,0)+\left(-\tfrac{\sqrt2}{2},\tfrac{\sqrt2}{2}\right) =\left(0,\tfrac{\sqrt2}{2}\right),\\ |\mathbf{z}|^2&=\tfrac12 . \end{aligned} \end{equation*}

(3) No: the cross terms \(2a_ia_j\,\mathbf{e_i}\cdot\mathbf{e_j}\) in \(|\sum_i a_i\mathbf{e_i}|^2\) can be negative. In the basis of part 2 take

\begin{equation*} \begin{aligned} \mathbf{y}&=\tfrac{\sqrt2}{10}\mathbf{e_1}+\mathbf{e_2} =\left(-\tfrac{2\sqrt2}{5},\tfrac{\sqrt2}{2}\right), &|\mathbf{y}|^2&=\tfrac{41}{50},\\ \mathbf{z}&=\tfrac{3\sqrt2}{5}\mathbf{e_1}+\tfrac65\,\mathbf{e_2} =\left(0,\tfrac{3\sqrt2}{5}\right), &|\mathbf{z}|^2&=\tfrac{36}{50}. \end{aligned} \end{equation*}

Here \(0<\frac{\sqrt2}{10}<\frac{3\sqrt2}{5}\) and \(0<1<\frac65\), yet \(|\mathbf{z}|<|\mathbf{y}|\).

(4) Solving \((x,y)=a_1\mathbf{e_1}+a_2\mathbf{e_2}=\left(a_1-\frac{\sqrt2}{2}a_2,\ \frac{\sqrt2}{2}a_2\right)\) gives

\begin{equation*} a_2=\sqrt2\,y,\qquad a_1=x+y . \end{equation*}

• (a) \((0,1)=1\cdot\mathbf{e_1}+\sqrt2\,\mathbf{e_2}\), i.e. \((1,\sqrt2)_e\).
• (b) \(\left(\frac{\sqrt2}{2},\frac{\sqrt2}{2}\right)=\sqrt2\,\mathbf{e_1}+1\cdot\mathbf{e_2}\), i.e. \((\sqrt2,1)_e\).
• (c) \((1,2)=3\,\mathbf{e_1}+2\sqrt2\,\mathbf{e_2}\), i.e. \((3,2\sqrt2)_e\).

With the Gaussian projection \(f(\mathbf v)=(\mathbf u_1\cdot\mathbf v,\dots,\mathbf u_k\cdot\mathbf v)\) of Section 2.7 (seed 0), the maximum over the \(\binom{20}{2}=190\) pairs of the relative error \(\big||f(\mathbf v_i)-f(\mathbf v_j)|-\sqrt k\,|\mathbf v_i-\mathbf v_j|\big|/(\sqrt k\,|\mathbf v_i-\mathbf v_j|)\) is:

import numpy as np
from scipy.spatial.distance import pdist
rng = np.random.default_rng(0)
d, n = 900, 20
X = rng.standard_normal((n, d))
X = 30 * X / np.linalg.norm(X, axis=1, keepdims=True)   # uniform on sphere of radius 30
D = pdist(X)                                             # original distances
for k in [100, 50, 10, 5, 4, 3, 2, 1]:
    U = rng.standard_normal((k, d))                      # rows u_1, ..., u_k
    Dk = pdist(X @ U.T)                                  # |f(v_i) - f(v_j)|
    rel = np.abs(Dk - np.sqrt(k) * D) / (np.sqrt(k) * D)
    print(k, 100 * rel.max(), 100 * np.median(rel))

The original distances lie in \([40.8,44.2]\), near \(30\sqrt2\) (Theorem 2.8). For fixed \(\mathbf w=\mathbf v_i-\mathbf v_j\) the vector \(f(\mathbf w)/|\mathbf w|\) is standard Gaussian in \(\mathbb R^k\) (proof of Theorem 2.10), so its length is \(\sqrt k\pm O(1)\) by Theorem 2.9 and the relative error has standard deviation about \(1/\sqrt{2k}\) (\(7\%\) at \(k=100\)): the error decays like \(1/\sqrt k\), as the Johnson-Lindenstrauss Lemma (Theorem 2.11) predicts, and is useless for \(k\le5\).

Both methods give 1000 vectors pairwise more than \(60^\circ\) apart; the projection method is slightly better on every measure. Over the \(\binom{1000}{2}=499{,}500\) pairs (seed 0, angles measured after normalising):

import numpy as np
rng = np.random.default_rng(0)
def stats(V):
    V = V / np.linalg.norm(V, axis=1, keepdims=True)
    C = np.abs(V @ V.T)[np.triu_indices(len(V), 1)]
    ang = np.degrees(np.arccos(C))
    return C.max(), ang.min(), np.sqrt((C**2).mean()), (ang < 80).mean()
Q, _ = np.linalg.qr(rng.standard_normal((1000, 100)))  # orthonormal basis of a random subspace
A = np.eye(1000) @ Q            # method 1: row i = coordinates of the projection of e_i
B = rng.standard_normal((1000, 100))                    # method 2
print(stats(A)); print(stats(B))

For method 2, \(E\cos^2\theta=1/100\) for independent Gaussian directions (Section 2.4.2), so the RMS is \(0.1\). For method 1 the rows \(\mathbf q_i\) of \(Q\) form a tight frame, \(\sum_i\mathbf q_i\mathbf q_i^T=Q^TQ=I_{100}\), hence

\begin{equation*} \sum_{i,j}(\mathbf q_i\cdot\mathbf q_j)^2=\operatorname{tr}\big((QQ^T)^2\big)=\operatorname{tr}(QQ^T)=100, \end{equation*}

the minimum possible for a rank-100 Gram matrix of trace 100; with \(|\mathbf q_i|^2\approx\frac1{10}\) the average off-diagonal \(\cos^2\theta\) is about \(\frac{9}{999}\), RMS \(0.0949\), matching the observed \(0.0950\).

Fit the line \(\mathbf x(t)=\mathbf a+\mathbf b\,t\) to all readings by least squares and report \(\hat{\mathbf a}+\hat{\mathbf b}\,n\). With readings \(\mathbf y_t=\mathbf a+\mathbf b t+\boldsymbol\eta_t\) at minutes \(t=1,\dots,n\) and independent spherical noise \(\boldsymbol\eta_t\sim N(\mathbf 0,\sigma^2I)\), the likelihood is

\begin{equation*} c\,\exp\Big(-\frac{1}{2\sigma^2}\sum_{t=1}^n|\mathbf y_t-\mathbf a-\mathbf b t|^2\Big), \end{equation*}

so, exactly as in Section 2.9, the maximum likelihood \((\mathbf a,\mathbf b)\) minimises the sum of squared distances. Setting the gradient to zero (as in Lemma 2.12), with \(\bar t=\frac{n+1}{2}\) and \(\bar{\mathbf y}=\frac1n\sum_t\mathbf y_t\),

\begin{equation*} \begin{aligned} \hat{\mathbf b}&=\frac{\sum_t (t-\bar t)(\mathbf y_t-\bar{\mathbf y})}{\sum_t (t-\bar t)^2},\qquad \hat{\mathbf a}=\bar{\mathbf y}-\hat{\mathbf b}\,\bar t,\\ \hat{\mathbf x}(n)&=\bar{\mathbf y}+\tfrac{n-1}{2}\,\hat{\mathbf b}. \end{aligned} \end{equation*}

This is unbiased, and since \(\sum_t(t-\bar t)^2=\frac{n(n^2-1)}{12}\) each coordinate has variance

\begin{equation*} \begin{aligned} \sigma^2\Big(\frac1n+\frac{(n-\bar t)^2}{\sum_t(t-\bar t)^2}\Big) &=\sigma^2\Big(\frac1n+\frac{3(n-1)}{n(n+1)}\Big)\\ &=\frac{2(2n-1)}{n(n+1)}\,\sigma^2\approx\frac{4\sigma^2}{n}, \end{aligned} \end{equation*}

far below the variance \(\sigma^2\) of the latest reading alone.

Write \(\varphi(x)=\frac{1}{\sqrt{2\pi}}e^{-x^2/2}\) and \(\Phi\) for its distribution function.

(1) The rectangle \([-1,1]\times[0,\varphi(1)]\), of area \(2\varphi(1)=\sqrt{2/(\pi e)}\approx0.4839\). A rectangle \([l,l+w]\times[0,h]\) lies under the curve iff \(h\le\min\big(\varphi(l),\varphi(l+w)\big)\) (\(\varphi\) is unimodal), and one endpoint is at distance at least \(w/2\) from \(0\), so

\begin{equation*} \text{area}\le w\,\varphi(w/2)=2a\varphi(a)\quad(a=w/2),\qquad \frac{d}{da}\,2a\varphi(a)=2(1-a^2)\varphi(a), \end{equation*}

which is positive for \(a<1\) and negative for \(a>1\); the maximum is at \(a=1\).

(2) The rectangle \([-a,a]\times[0,\frac1{2a}]\) with \(a\approx1.486\) (width \(2.973\), height \(0.3364\)), whose symmetric difference with the Gaussian is about \(0.370\). Both regions have area \(1\), so

\begin{equation*} |\,\text{Gaussian}\ \triangle\ \text{rectangle}\,|=\int|\varphi-R|=2-2\int\min(\varphi,R)\,dx, \end{equation*}

and we maximise the overlap. For fixed width \(w\) and height \(h=1/w\) the overlap \(\int_l^{l+w}\min(\varphi,h)\) integrates the symmetric unimodal function \(\min(\varphi,h)\) over a window of length \(w\), which is largest for the centred window \(l=-a\). Since part 1 rules out \(h\le\varphi(a)\), two cases remain:

• (i) \(h\ge\varphi(0)\), i.e. \(a\le\sqrt{\pi/2}\): overlap \(=2\Phi(a)-1\), increasing in \(a\), so not optimal;
• (ii) \(h<\varphi(0)\): with \(\varphi( c)=h=\frac1{2a}\), \(0<c<a\), overlap

\begin{equation*} G(a)=2ch+2\big(\Phi(a)-\Phi( c)\big). \end{equation*}

The \(c\)-derivative \(2h-2\varphi( c)\) vanishes, so

\begin{equation*} \begin{aligned} G^{\prime}(a)=2\varphi(a)-\frac{c}{a^2}=0&\iff 2a^2\varphi(a)=c,\\ c&=\sqrt{2\ln\frac{2a}{\sqrt{2\pi}}}. \end{aligned} \end{equation*}

Solving numerically gives \(a=1.4864\), \(c=0.5841\), \(h=0.3364\), \(G=0.8149\), so the symmetric difference is \(2-2G=0.3701\) (confirmed by direct numerical integration, which also increases when the rectangle is shifted or rescaled).

By independence the variances add, so

\begin{equation*} \operatorname{var}(m_s)=\frac1{n^2}\sum_{i=1}^n\operatorname{var}(x_i)=\frac{\sigma^2}{n}, \qquad E(m_s)=\mu . \end{equation*}

Using \(x_i-m_s=(x_i-\mu)-(m_s-\mu)\) and \(\sum_i(x_i-\mu)=n(m_s-\mu)\),

\begin{equation*} \begin{aligned} \sum_{i=1}^n(x_i-m_s)^2&=\sum_i(x_i-\mu)^2-2(m_s-\mu)\sum_i(x_i-\mu)+n(m_s-\mu)^2\\ &=\sum_i(x_i-\mu)^2-n(m_s-\mu)^2 , \end{aligned} \end{equation*}

hence

\begin{equation*} E(\sigma_s^2)=\frac1n\Big(n\sigma^2-n\operatorname{var}(m_s)\Big) =\frac{n-1}{n}\,\sigma^2 . \end{equation*}

Therefore \(\frac{1}{n-1}\sum_i(x_i-m_s)^2=\frac{n}{n-1}\sigma_s^2\) has expectation exactly \(\sigma^2\), i.e. dividing by \(n-1\) gives an unbiased estimator.

With seed 0 the ten values are

\begin{equation*} 0.126,\ -0.132,\ 0.640,\ 0.105,\ -0.536,\ 0.362,\ 1.304,\ 0.947,\ -0.704,\ -1.265, \end{equation*}

the averaged center is \(m=0.0847\), and the three variance estimates are

\begin{equation*} \begin{aligned} \tfrac1{10}\textstyle\sum_i x_i^2&=0.5566 &&\text{(real center 0)},\\ \tfrac1{10}\textstyle\sum_i (x_i-m)^2&=0.5494 &&\text{(estimated center, }n=10),\\ \tfrac1{9}\textstyle\sum_i (x_i-m)^2&=0.6105 &&\text{(estimated center, }n=9). \end{aligned} \end{equation*}

import numpy as np
rng = np.random.default_rng(0)
x = rng.standard_normal(10)
m = x.mean()
print(m, (x**2).sum() / 10, ((x - m)**2).sum() / 10, ((x - m)**2).sum() / 9)
X = rng.standard_normal((100000, 10)); M = X.mean(1, keepdims=True)   # repeat 100000 times
for v in (X**2).mean(1), ((X - M)**2).sum(1) / 10, ((X - M)**2).sum(1) / 9:
    print(v.mean(), ((v - 1)**2).mean())

The estimated-center sum is always the smaller one, by Lemma 2.12 (the centroid minimises the sum of squared distances); precisely \(\sum_i(x_i-m)^2=\sum_i x_i^2-10m^2=5.566-0.072=5.494\). One sample of ten is too noisy to rank the estimators (all three are well below the true value \(1\) here), so the run was repeated \(100{,}000\) times:

Dividing by \(10\) after estimating the center is biased low by the factor \(\frac{n-1}{n}=0.9\) of Exercise 2.41; dividing by \(9\) removes the bias, at the price of slightly larger variance. (Theory column: \(\sum_i x_i^2\sim\chi^2_{10}\), \(\sum_i(x_i-m)^2\sim\chi^2_9\), and \(\operatorname{var}\chi^2_k=2k\).)

\(n \ge 2\ln(200d)/\varepsilon^2 = O(\log d/\varepsilon^2)\) samples suffice for the \(\ell_\infty\) guarantee, and \(n \ge 100\,d/\varepsilon^2 = O(d/\varepsilon^2)\) suffice for the \(\ell_2\) guarantee, where \(\mathbf{m}_s=\frac1n\sum_{j=1}^n \mathbf{x}_j\) is the sample mean.

Since the \(\mathbf{x}_j\) are independent \(N(\boldsymbol{\mu},I_d)\), the vector \(\mathbf{z}=\sqrt n\,(\mathbf{m}_s-\boldsymbol{\mu})\) is a standard Gaussian in \(\mathbf{R}^d\), so each coordinate \(z_i\sim N(0,1)\). For \(Z\sim N(0,1)\) and \(\lambda>0\), \(E(e^{\lambda Z})=e^{\lambda^2/2}\), so Markov’s inequality (Theorem 2.1, applied to the nonnegative \(e^{\lambda Z}\)) with \(\lambda=t\) gives \(\text{Prob}(Z\ge t)\le e^{-t^2/2}\), hence \(\text{Prob}(|Z|\ge t)\le 2e^{-t^2/2}\). A union bound over the \(d\) coordinates with \(t=\varepsilon\sqrt n\) gives

\begin{equation*} \begin{aligned} \text{Prob}\big(\|\boldsymbol{\mu}-\mathbf{m}_s\|_\infty>\varepsilon\big) &\le \sum_{i=1}^d \text{Prob}\big(|z_i|>\varepsilon\sqrt n\big) \le 2d\,e^{-n\varepsilon^2/2}\\ &\le 0.01 \iff n\ge \frac{2\ln(200d)}{\varepsilon^2}. \end{aligned} \end{equation*}

For the Euclidean norm, \(E\|\mathbf{m}_s-\boldsymbol{\mu}\|_2^2=\frac1n E|\mathbf{z}|^2=\frac dn\), so Markov’s inequality applied to the nonnegative variable \(\|\mathbf{m}_s-\boldsymbol{\mu}\|_2^2\) gives

\begin{equation*} \text{Prob}\big(\|\boldsymbol{\mu}-\mathbf{m}_s\|_2>\varepsilon\big) \le \frac{d}{n\varepsilon^2}\le 0.01 \iff n\ge \frac{100\,d}{\varepsilon^2}. \end{equation*}

Method (2): For \(\beta=\sqrt{\ln(300)/c}\) (with \(c\) the constant of Theorem 2.9, and \(d\ge\beta^2\) so that \(\beta\le\sqrt d\)), the Gaussian Annulus Theorem gives \(|\mathbf{z}|\le\sqrt d+\beta\) with probability at least \(1-3e^{-c\beta^2}=0.99\). So \(n\ge(\sqrt d+\beta)^2/\varepsilon^2=(1+o(1))\,d/\varepsilon^2\) suffices.

With seed 2044 the estimates are \(m=4.887\), then \(5.828\), \(5.816\) and \(6.462\) for \(\sigma^2\). The density is \(N(5,9)\), so we draw ten points with \(\sigma=3\):

import numpy as np
rng = np.random.default_rng(2044)
x = rng.normal(loc=5, scale=3, size=10)
m = x.mean()
print(m, abs(m - 5))                 # (a)
print(np.sum((x - 5)**2) / 10)       # (b)
print(np.sum((x - m)**2) / 10)       # (c)
print(np.sum((x - m)**2) / 9)        # (d)

The sample, rounded, is \(2.651, 5.091, 3.713, 3.927, 4.670, 6.621, 11.417, 3.355, 3.757, 3.666\).

(a) \(m=4.887\), so \(|m-5|=0.113\). The standard deviation of \(m\) is \(3/\sqrt{10}=0.949\) (Theorem 2.4), so this draw came out unusually close.

(b) \(\frac1{10}\sum(x_i-5)^2=5.828\), an error of \(-3.17\). The estimator is unbiased, but its standard deviation is \(\sqrt{2\cdot 81/10}=4.02\), so this error is typical for ten points.

(c) \(\frac1{10}\sum(x_i-m)^2=5.816\), an error of \(-3.18\). By Lemma 2.12, \(m\) minimizes \(\sum(x_i-\mu)^2\), and in fact

\begin{equation*} \sum_{i=1}^{10}(x_i-m)^2=\sum_{i=1}^{10}(x_i-5)^2-10(m-5)^2 = 58.283-0.128=58.155, \end{equation*}

so (c) is never larger than (b). By Exercise 2.41 its mean is \(\frac{9}{10}\cdot 9=8.1\).

(d) \(\frac19\sum(x_i-m)^2=6.462\), an error of \(-2.54\). This is the unbiased estimator (Exercise 2.41), and it is the closest of the three here.

(1) Volume lives at the surface. Shrinking a body by a factor \(1-\varepsilon\) multiplies its volume by \((1-\varepsilon)^d\le e^{-\varepsilon d}\) (Section 2.3). So almost all of the unit ball lies in a shell of width \(O(1/d)\), and the volume \(V(d)=\pi^{d/2}/\Gamma(\frac d2+1)\) of the unit ball goes to \(0\) (Section 2.4.1).

(2) Volume also lives at every equator. By Theorem 2.7, all but a \(\frac2c e^{-c^2/2}\) fraction of the ball has \(|x_1|\le c/\sqrt{d-1}\). Since any direction can serve as north, two random points are nearly orthogonal with high probability (Theorem 2.8).

(3) A Gaussian is a thin annulus, not a bump. By the Gaussian Annulus Theorem 2.9, a standard \(d\)-dimensional Gaussian puts nearly all its mass within \(O(1)\) of the sphere of radius \(\sqrt d\), even though the density is largest at the origin. A spherical Gaussian therefore looks like a uniform distribution on a sphere.

(4) Distances survive drastic compression. By the Random Projection Theorem 2.10 and the Johnson-Lindenstrauss Lemma (Theorem 2.11), projecting onto \(k=O(\varepsilon^{-2}\ln n)\) random Gaussian directions keeps all \(\binom n2\) pairwise distances within a factor \(1\pm\varepsilon\) (after scaling by \(\sqrt k\)). The target dimension depends on \(n\), not on \(d\), which makes nearest-neighbor search in huge dimensions practical.

(5) Concentration turns randomness into near-certainty. The Law of Large Numbers (Theorem 2.4) and the Master Tail Bounds Theorem 2.5 make sums of many independent coordinates concentrate sharply. Two unit-variance Gaussians whose centers are \(\Omega(d^{1/4})\) apart can therefore be separated by distance alone (Section 2.8). Fitting a spherical Gaussian reduces to the sample mean and sample spread (Lemmas 2.12 and 2.13), with the \(n-1\) correction of Exercise 2.41.

Your intuition was built in three dimensions, and in a thousand dimensions it will lie to you. That is exactly why the subject is worth learning.

Start with an orange. In three dimensions the peel is a small part of the fruit. In a hundred dimensions, shrinking an orange by just \(5\%\) throws away more than \(99\%\) of its volume, because volume scales like \((0.95)^{100}\approx 0.006\). A high-dimensional orange is almost all peel. Stranger still, almost all of it also sits near the equator, whichever direction you call north. Pick two random directions and they will be nearly perpendicular.

Who cares about a thousand dimensions? You do, every day. A photo is a point whose coordinates are its pixel values. A song, a document described by its word counts, and your movie ratings are points too. Search engines, recommendation systems and face recognition all work by measuring distances between such points in spaces with thousands or millions of coordinates.

High dimensions are not only a curse; they also contain hidden gifts. A bell curve in high dimensions is not a bump but a thin soap bubble of radius \(\sqrt d\). Because of this, two overlapping clouds of data can still be told apart. And the Johnson-Lindenstrauss Lemma says you can squash a million-dimensional data set into a few hundred random directions while keeping every pairwise distance almost exact. That one fact makes many large-scale computations possible.

The tools are things you already know or will soon learn: calculus, linear algebra and basic probability. Combined, they let you reason confidently about places you can never draw. Few subjects give you surprises this deep for so little machinery. Come and find out where your intuition breaks.

With \(\bar x=\frac1n\sum_i x_i\) and \(\bar y=\frac1n\sum_i y_i\), the minimizer is

\begin{equation*} m=\frac{\sum_i (x_i-\bar x)(y_i-\bar y)}{\sum_i (x_i-\bar x)^2} =\frac{n\sum_i x_iy_i-\sum_i x_i\sum_i y_i}{n\sum_i x_i^2-\left(\sum_i x_i\right)^2}, \qquad b=\bar y-m\bar x, \end{equation*}

valid whenever the \(x_i\) are not all equal (if all \(x_i\) equal \(\bar x\), every pair with \(b=\bar y-m\bar x\) is optimal).

The error \(E(m,b)=\frac1n\sum_i (y_i-mx_i-b)^2\) is a convex quadratic in \((m,b)\), so the minimum is where both partial derivatives vanish:

\begin{equation*} \begin{aligned} \frac{\partial E}{\partial b}=0 &\iff \sum_i (y_i-mx_i-b)=0 \iff b=\bar y-m\bar x,\\ \frac{\partial E}{\partial m}=0 &\iff \sum_i x_i(y_i-mx_i-b)=0 \iff \sum_i x_i\bigl((y_i-\bar y)-m(x_i-\bar x)\bigr)=0, \end{aligned} \end{equation*}

and since \(\sum_i \bar x\bigl((y_i-\bar y)-m(x_i-\bar x)\bigr)=0\) we may replace \(x_i\) by \(x_i-\bar x\) in the last sum, giving \(\sum_i (x_i-\bar x)(y_i-\bar y)=m\sum_i (x_i-\bar x)^2\), i.e. the formula for \(m\). The Hessian is \(\frac2n\begin{pmatrix}\sum x_i^2 & \sum x_i\\ \sum x_i & n\end{pmatrix}\), positive definite exactly when the \(x_i\) are not all equal (its determinant is \(\frac4{n^2}\, n\sum_i(x_i-\bar x)^2>0\)), so this critical point is the unique minimizer.

Center the data and take the last right-singular vector: \((a_1,\dots,a_5)=v_5\) and \(a_6=v_5\cdot c\). Here the \(n\) people are the rows \(p_1,\dots,p_n\in\mathbb{R}^5\), \(c=\frac1n\sum_i p_i\) is the centroid, and \(v_1,\dots,v_5\) are the right-singular vectors of the centered \(n\times 5\) matrix \(A\) with rows \(p_i-c\). By Lemma 3.14 the best 4-dimensional affine subspace passes through \(c\), so it is \(c+V_4\) for the best 4-dimensional subspace \(V_4\) of the centered data. By Theorem 3.1 that is \(V_4=\operatorname{span}(v_1,\dots,v_4)\). Its unit normal is \(v_5\), which gives

\begin{equation*} c+V_4=\{x\in\mathbb{R}^5 : v_5\cdot(x-c)=0\}=\{x : v_5\cdot x=v_5\cdot c\}, \end{equation*}

and the sum of squared perpendicular distances is \(\sum_i (v_5\cdot(p_i-c))^2=\sigma_5^2\). A good 4-dimensional fit means \(\sigma_5\) is small compared with \(\sigma_1,\dots,\sigma_4\). In practice each variable is first rescaled (say to unit variance), since the five variables are in incomparable units.

Why perpendicular distance: (i) Symmetry. Vertical distance along, say, the blood-pressure axis means regressing blood pressure on the other four. That treats blood pressure as the only noisy variable and the other four as exact. The fitted sheet then depends on which variable is singled out (five choices give five different answers). Perpendicular distance treats all five measurements, each of them noisy, in the same way, and the answer does not depend on how the axes are labelled (it is invariant under rotations). (ii) Degeneracy. If the true relation hardly involves the chosen variable (\(a_5\approx 0\)), the sheet is nearly parallel to that axis. Vertical distances to it are then huge or undefined (when \(a_5=0\) the sheet cannot even be written as blood pressure \(=\) a function of the rest), while perpendicular distances stay small and the SVD finds the sheet anyway. (iii) Optimality. The perpendicular criterion is exactly the one for which the SVD is optimal (Theorem 3.1, Lemma 3.14). It finds the sheet carrying the most variance of the data, which is the “4-dimensional sheet” picture in the question.

The best fit lines are (1) \(x+y=8\); (2) every line through \((5,3)\); (3) \(y=-2x+11\). Only in case 2 does a best fit line, namely \(y=\tfrac35x\), pass through the origin. In each case let \(A\) be the matrix of centered points. By Lemma 3.13 the best fit line passes through the centroid \(c\). Its direction is the top right-singular vector of \(A\), i.e. the top eigenvector of \(A^TA\) (Section 3.8). The best fit line of the centered points is the line through \(0\) in that direction.

(1) Here \(c=(5,3)\) and the centered points are \(\pm(-1,1)\), so \(A^TA=\begin{pmatrix}2&-2\\-2&2\end{pmatrix}\), with top eigenvector \(\frac1{\sqrt2}(1,-1)\) (eigenvalue \(4\)). The centered fit is \(y=-x\), which contains both centered points. The original fit is \((5,3)+\lambda(1,-1)\), i.e. \(x+y=8\). It does not pass through the origin.

(2) Here \(c=(5,3)\) and the centered points are \((\pm1,\pm1)\), so \(A^TA=4I\). Every unit \(v\) gives \(|Av|^2=4\), so every line through \(0\) is a best fit for the centered data. Every line through \((5,3)\) is a best fit for the original data, each with sum of squared distances \(|A|_F^2-4=4\). The fit is not unique, and exactly one best fit line, \(y=\tfrac35 x\), passes through the origin.

(3) Here \(c=(4,3)\) and the centered points are \((-1,-\tfrac12),\ (-1,2),\ (1,-2),\ (1,\tfrac12)\), so

\begin{equation*} A^TA=\begin{pmatrix}4&-3\\-3&8.5\end{pmatrix},\qquad \lambda^2-12.5\lambda+25=0\ \Rightarrow\ \lambda\in\{10,\ 2.5\}. \end{equation*}

For \(\lambda=10\), \((4-10)x-3y=0\) gives \(v_1=\frac1{\sqrt5}(1,-2)\). The centered fit is \(y=-2x\) (sum of squared distances \(2.5\)). The original fit is \(y-3=-2(x-4)\), i.e.

\begin{equation*} y=-2x+11, \end{equation*}

which does not pass through the origin.

In (1) every line through the origin is a best fit, so the fit is not unique. In (2) the unique best fit is the \(x\)-axis. By Section 3.2 the best fit line through \(0\) maximizes \(|Av|^2=v^TA^TAv\) over unit vectors \(v\), where \(A\) has the points as rows.

(1) Here \(A=I\), so \(|Av|^2=v_1^2+v_2^2=1\) for every unit \(v\). Every line through the origin has the same sum of squared distances \(2-1=1\), so the best fit line is not unique.

(2) Here \(A=\begin{pmatrix}0&1\\2&0\end{pmatrix}\), so for \(|v|=1\)

\begin{equation*} |Av|^2=4v_1^2+v_2^2=1+3v_1^2\le 4, \end{equation*}

with equality iff \(v=\pm(1,0)\). The best fit line is the \(x\)-axis \(\{y=0\}\), with sum of squared distances \(5-4=1\). It is unique because the maximizer \(v\) is unique up to sign (\(\sigma_1=2>\sigma_2=1\)).

In both cases \(v_1=\frac1{\sqrt2}(1,1)\) and \(v_2=\frac1{\sqrt2}(1,-1)\). The \(v_i\) are the eigenvectors of \(M^TM\) and \(\sigma_i^2\) its eigenvalues (Section 3.8). The \(u_i=Mv_i/\sigma_i\) come from Lemma 3.10, and \(M=\sum_i\sigma_iu_iv_i^T\) by Theorem 3.4.

(a) Here

\begin{equation*} M^TM=\begin{pmatrix}10&1\\1&10\end{pmatrix},\qquad \sigma_1^2=11,\ \sigma_2^2=9, \end{equation*}

so \(\sigma_1=\sqrt{11}\), \(\sigma_2=3\), and

\begin{equation*} u_1=\frac{M v_1}{\sqrt{11}}=\frac{1}{\sqrt{22}}\begin{pmatrix}2\\3\\3\end{pmatrix},\qquad u_2=\frac{M v_2}{3}=\frac{1}{\sqrt2}\begin{pmatrix}0\\-1\\1\end{pmatrix}. \end{equation*}

Hence

\begin{equation*} M=\sqrt{11}\,u_1v_1^T+3\,u_2v_2^T =\begin{pmatrix}\frac{2}{\sqrt{22}}&0\\ \frac{3}{\sqrt{22}}&-\frac1{\sqrt2}\\ \frac{3}{\sqrt{22}}&\frac1{\sqrt2}\end{pmatrix} \begin{pmatrix}\sqrt{11}&0\\0&3\end{pmatrix} \begin{pmatrix}\frac1{\sqrt2}&\frac1{\sqrt2}\\ \frac1{\sqrt2}&-\frac1{\sqrt2}\end{pmatrix}. \end{equation*}

(b) Here

\begin{equation*} M^TM=\begin{pmatrix}14&6\\6&14\end{pmatrix},\qquad \sigma_1^2=20,\ \sigma_2^2=8, \end{equation*}

so \(\sigma_1=2\sqrt5\), \(\sigma_2=2\sqrt2\), and

\begin{equation*} u_1=\frac{M v_1}{2\sqrt5}=\frac{1}{\sqrt{10}}\begin{pmatrix}1\\1\\2\\2\end{pmatrix},\qquad u_2=\frac{M v_2}{2\sqrt2}=\frac{1}{2}\begin{pmatrix}-1\\1\\-1\\1\end{pmatrix}. \end{equation*}

Hence

\begin{equation*} M=2\sqrt5\,u_1v_1^T+2\sqrt2\,u_2v_2^T =\begin{pmatrix}\frac1{\sqrt{10}}&-\frac12\\ \frac1{\sqrt{10}}&\frac12\\ \frac2{\sqrt{10}}&-\frac12\\ \frac2{\sqrt{10}}&\frac12\end{pmatrix} \begin{pmatrix}2\sqrt5&0\\0&2\sqrt2\end{pmatrix} \begin{pmatrix}\frac1{\sqrt2}&\frac1{\sqrt2}\\ \frac1{\sqrt2}&-\frac1{\sqrt2}\end{pmatrix}. \end{equation*}

(Check!) Multiplying out recovers \(M\) in both cases.

(1) The estimate is \(v_1\approx\frac{1}{\sqrt{4097}}(1,64)\approx(0.0156,\ 0.9999)\). The columns of \(A\) are orthogonal (\(1\cdot2-1\cdot2-1\cdot2+1\cdot2=0\)), so with the power method of Section 3.7,

\begin{equation*} B=A^TA=\begin{pmatrix}4&0\\0&16\end{pmatrix},\qquad B^3x=\begin{pmatrix}64\\4096\end{pmatrix},\qquad \frac{B^3x}{|B^3x|}=\frac{1}{\sqrt{4097}}\begin{pmatrix}1\\64\end{pmatrix}. \end{equation*}

(2) \(B\) is diagonal, so its eigenvectors are the coordinate vectors and, by Lemma 3.10,

\begin{equation*} \begin{aligned} v_1&=(0,1),& \sigma_1&=\sqrt{16}=4,& u_1&=\tfrac{Av_1}{4}=\tfrac12(1,1,-1,-1)^T,\\ v_2&=(1,0),& \sigma_2&=\sqrt4=2,& u_2&=\tfrac{Av_2}{2}=\tfrac12(1,-1,1,-1)^T, \end{aligned} \end{equation*}

so \(A=4u_1v_1^T+2u_2v_2^T\).

(3) \(v_1=(0,1)\) is the dominant direction in restaurant space, the “taste axis” that explains most of the variation in ratings. Here it is simply restaurant 2, the one people disagree about most strongly. \(u_1=\frac12(1,1,-1,-1)\) gives each person’s score on that axis (up to the factor \(\sigma_1\)): persons 1 and 2 love restaurant 2 and persons 3 and 4 dislike it, so \(u_1\) splits the diners into two taste groups. The gap \(\sigma_1-\sigma_2=2\) measures how dominant this single factor is. A large gap means the ratings are close to the rank-one matrix \(\sigma_1u_1v_1^T\) (Theorem 3.6: the error is \(\sigma_2^2\) out of \(\sigma_1^2+\sigma_2^2\) in squared Frobenius norm), so one “liking” score per person explains most of the data. A small gap means a second taste factor (here restaurant 1, which splits persons \(\{1,3\}\) from \(\{2,4\}\)) matters almost as much.

Orthonormal rows mean \(AA^T=I\), and for a square matrix this forces \(A^TA=I\), i.e. orthonormal columns:

\begin{equation*} \begin{aligned} AA^T=I &\Rightarrow \det(A)\det(A^T)=\det(A)^2=1\ \Rightarrow\ A \text{ is invertible},\\ &\Rightarrow A^T=A^{-1}(AA^T)=A^{-1}\\ &\Rightarrow A^TA=A^{-1}A=I. \end{aligned} \end{equation*}

The \((i,j)\) entry of \(A^TA\) is the dot product of columns \(i\) and \(j\), so the columns are orthonormal.

The singular vectors are \(\mathbf{v}_i=\mathbf{w}_i\) (up to sign) with singular values \(\sigma_i=a_in/k\), where \(\mathbf{w}_i\) is \((k/n)^{1/2}\) on block \(B_i\) and \(0\) elsewhere; when all \(a_i=a\), the singular vectors are an arbitrary orthonormal basis of \(W=\operatorname{span}(\mathbf{w}_1,\ldots,\mathbf{w}_k)\).

Each block has size \(m=n/k\), and the \(i^{\text{th}}\) block is \(a_i\mathbf{1}_{B_i}\mathbf{1}_{B_i}^T=a_im\,\mathbf{w}_i\mathbf{w}_i^T\). The \(\mathbf{w}_i\) are orthonormal (disjoint supports, \(m\cdot\frac kn=1\)), so with \(c_j=\mathbf{w}_j^T\mathbf{v}\) and \(\sigma_j=a_jn/k\),

\begin{equation*} A=\sum_{j=1}^k\sigma_j\mathbf{w}_j\mathbf{w}_j^T,\qquad |A\mathbf{v}|^2=\sum_{j=1}^k\sigma_j^2c_j^2,\qquad \sum_{j=1}^kc_j^2\le|\mathbf{v}|^2 . \end{equation*}

Now run the greedy definition of Section 3.3 by induction on \(i\). If \(\mathbf{v}\) is a unit vector orthogonal to \(\mathbf{w}_1,\ldots,\mathbf{w}_{i-1}\), then \(c_j=0\) for \(j<i\), and

\begin{equation*} \begin{aligned} \sigma_i^2-|A\mathbf{v}|^2 &=\sigma_i^2\Bigl(1-\sum_{j\ge i}c_j^2\Bigr) +\sum_{j>i}(\sigma_i^2-\sigma_j^2)c_j^2\ \ge 0 . \end{aligned} \end{equation*}

Since \(\sigma_i>\sigma_j\) for \(j>i\), equality forces \(c_j=0\) for \(j>i\) and \(c_i^2=1\), i.e. \(\mathbf{v}=\pm\mathbf{w}_i\). So the maximum is \(\sigma_i\), attained only at \(\mathbf{v}_i=\pm\mathbf{w}_i\). After \(k\) steps any \(\mathbf{v}\perp\mathbf{w}_1,\ldots,\mathbf{w}_k\) has \(A\mathbf{v}=0\), so the process stops: there are exactly \(k\) singular vectors, and \(\mathbf{u}_i=A\mathbf{v}_i/\sigma_i=\mathbf{v}_i\).

If \(a_1=\cdots=a_k=a\), then \(A=\sigma P_W\), where \(\sigma=an/k\) and \(P_W=\sum_j\mathbf{w}_j\mathbf{w}_j^T\) is the orthogonal projection onto \(W\). So \(|A\mathbf{v}|=\sigma|P_W\mathbf{v}|\), which equals \(\sigma\) for every unit \(\mathbf{v}\in W\) and is smaller otherwise. At step \(i\) the greedy algorithm can therefore pick any unit vector of \(W\) orthogonal to \(\mathbf{v}_1,\ldots,\mathbf{v}_{i-1}\). The possible singular vector lists \((\mathbf{v}_1,\ldots,\mathbf{v}_k)\) are exactly the orthonormal bases of \(W\), i.e. \(\mathbf{v}_i=\sum_jQ_{ji}\mathbf{w}_j\) with \(Q\) any \(k\times k\) orthogonal matrix. All singular values equal \(an/k\), and \(\mathbf{u}_i=\mathbf{v}_i\). Each such vector is constant on every block, but it can mix blocks, so the singular vectors no longer single out individual blocks; only their span \(W\) is determined.

\(\mathbf{v}_1\) is the dominant topic, a nonnegative weighting of the terms. \(\mathbf{u}_1\) scores each document by how strongly it expresses that topic, and this is the intrinsic-relevance ranking of Section 3.9.4.

As in Section 3.5, the \(n\times d\) matrix \(A\) has a row \(\mathbf{a}_i\) for each document, and \(a_{ij}\) is the number of times term \(j\) occurs in document \(i\). Then

\begin{equation*} \mathbf{v}_1=\arg\max_{|\mathbf{v}|=1}\sum_{i=1}^n(\mathbf{a}_i\cdot\mathbf{v})^2,\qquad \mathbf{u}_1=\frac{A\mathbf{v}_1}{\sigma_1},\qquad A^T\mathbf{u}_1=\sigma_1\mathbf{v}_1 , \end{equation*}

where the last identity is Lemma 3.10. (i) Right singular vector. \(\mathbf{v}_1\in\mathbb{R}^d\) spans the best-fit line through the documents (Theorem 3.1 with \(k=1\)). It is the single synthetic query that agrees best, in the sum-of-squares sense, with the whole collection. The entries of \(A\) are nonnegative, so \(|(A|\mathbf{v}|)_i|\ge|(A\mathbf{v})_i|\) for every \(i\), and \(|\mathbf{v}_1|\) (taken coordinatewise) is also a maximizer. Hence we may take \(\mathbf{v}_1\ge0\). Its large coordinates are the terms that co-occur heavily across the collection, i.e. the vocabulary of the main topic. (ii) Left singular vector. \(\mathbf{u}_1\ge0\) and \((\mathbf{u}_1)_i=\mathbf{a}_i\cdot\mathbf{v}_1/\sigma_1\), so \((\mathbf{u}_1)_i\) is proportional to the projection of document \(i\) onto the topic direction. This is the score by which Section 3.9.4 ranks documents. By Exercise 3.11, \(\mathbf{u}_1\) also maximizes \(|\mathbf{u}^TA|\), so it gives the weighting of the documents whose combined term vector \(\mathbf{u}^TA\) is longest. The two vectors reinforce each other: \((\mathbf{v}_1)_j\propto\sum_i a_{ij}(\mathbf{u}_1)_i\), so a term is important when it occurs often in relevant documents, and \((\mathbf{u}_1)_i\propto\sum_ja_{ij}(\mathbf{v}_1)_j\), so a document is relevant when it uses important terms often. This is the same circular hub/authority definition as in Section 3.9.4.

Compare entries. Let \(X\) be \(n\times r\) and \(Y\) be \(d\times r\), so \(X_{ai}=(\mathbf{x}_i)_a\) and \(Y_{bi}=(\mathbf{y}_i)_b\), and \(C_{ij}=c_i\delta_{ij}\). For every \(a,b\),

\begin{equation*} (XCY^T)_{ab}=\sum_{i=1}^r\sum_{j=1}^rX_{ai}C_{ij}Y_{bj} =\sum_{i=1}^rc_i(\mathbf{x}_i)_a(\mathbf{y}_i)_b =\Bigl(\sum_{i=1}^rc_i\mathbf{x}_i\mathbf{y}_i^T\Bigr)_{ab}. \end{equation*}

Hence \(XCY^T=\sum_{i=1}^rc_i\mathbf{x}_i\mathbf{y}_i^T\). \(\blacksquare\)

\(\mathbf{u}_1^TA=\sigma_1\mathbf{v}_1^T\). The \(\mathbf{u}_i\) are orthonormal (Theorem 3.7) and so are the \(\mathbf{v}_i\). For any unit \(\mathbf{u}\),

\begin{equation*} \begin{aligned} \mathbf{u}^TA&=\sum_{i=1}^r\sigma_i(\mathbf{u}^T\mathbf{u}_i)\mathbf{v}_i^T \quad\Rightarrow\quad \mathbf{u}_1^TA=\sigma_1\mathbf{v}_1^T,\ \ |\mathbf{u}_1^TA|=\sigma_1,\\ |\mathbf{u}^TA|^2&=\sum_{i=1}^r\sigma_i^2(\mathbf{u}^T\mathbf{u}_i)^2 \le\sigma_1^2\sum_{i=1}^r(\mathbf{u}^T\mathbf{u}_i)^2\le\sigma_1^2|\mathbf{u}|^2=\sigma_1^2 . \end{aligned} \end{equation*}

The last inequality is Bessel’s inequality for the orthonormal set \(\{\mathbf{u}_i\}\). Hence \(\max_{|\mathbf{u}|=1}|\mathbf{u}^TA|\le\sigma_1\), and \(\mathbf{u}=\mathbf{u}_1\) attains this bound. \(\blacksquare\)

(1) Write \(A=\sum_i\sigma_i\mathbf{u}_i\mathbf{v}_i^T\) (Theorem 3.4). The \(\mathbf{u}_i\) are orthonormal (Theorem 3.7), so

\begin{equation*} \begin{aligned} A^TA&=\sum_{i,j}\sigma_i\sigma_j\mathbf{v}_i(\mathbf{u}_i^T\mathbf{u}_j)\mathbf{v}_j^T =\sum_{i=1}^r\sigma_i^2\mathbf{v}_i\mathbf{v}_i^T . \end{aligned} \end{equation*}

(2) The \(\mathbf{v}_i\) are orthonormal, so \(A^TA\mathbf{v}_j=\sum_i\sigma_i^2\mathbf{v}_i(\mathbf{v}_i^T\mathbf{v}_j)=\sigma_j^2\mathbf{v}_j\). Thus \(\mathbf{v}_j\) is an eigenvector of \(A^TA\) with eigenvalue \(\sigma_j^2\). (3) By part 1, completing \(\mathbf{v}_1,\ldots,\mathbf{v}_r\) to an orthonormal basis gives the spectral decomposition of \(A^TA\). So for every SVD of \(A\), the list \(\sigma_1^2\ge\cdots\ge\sigma_r^2\) is the list of nonzero eigenvalues of \(A^TA\) in decreasing order, with multiplicity. The singular values are therefore determined by \(A\). Under the hypothesis, each eigenspace \(E_j=\{\mathbf{x}:A^TA\mathbf{x}=\sigma_j^2\mathbf{x}\}\) is a line \(\mathbb{R}\mathbf{e}_j\). Then:

• (i) any right singular vector \(\mathbf{v}_j\) is a unit vector in \(E_j\) by part 2, so \(\mathbf{v}_j=\pm\mathbf{e}_j/|\mathbf{e}_j|\);
• (ii) the corresponding left singular vector is \(\mathbf{u}_j=A\mathbf{v}_j/\sigma_j\), so it is determined up to the same sign.

So the singular vectors are unique up to sign. \(\blacksquare\)

\begin{equation*} \|A_k\|_F^2=\sum_{i=1}^k\sigma_i^2,\qquad \|A_k\|_2^2=\sigma_1^2,\qquad \|A-A_k\|_F^2=\sum_{i=k+1}^r\sigma_i^2,\qquad \|A-A_k\|_2^2=\sigma_{k+1}^2 . \end{equation*}

Let \(S\subseteq\{1,\ldots,r\}\) and \(B=\sum_{i\in S}\sigma_i\mathbf{u}_i\mathbf{v}_i^T\). The \(\mathbf{u}_i\) are orthonormal (Theorem 3.7) and so are the \(\mathbf{v}_i\), so with \(\mathbf{e}_j\) the standard basis,

\begin{equation*} \begin{aligned} |B\mathbf{v}|^2&=\sum_{i\in S}\sigma_i^2(\mathbf{v}_i^T\mathbf{v})^2 \le\max_{i\in S}\sigma_i^2\,|\mathbf{v}|^2 \quad(\text{equality at }\mathbf{v}=\mathbf{v}_{\min S}),\\ \|B\|_F^2&=\sum_j|B\mathbf{e}_j|^2 =\sum_{i\in S}\sigma_i^2\sum_j(\mathbf{v}_i^T\mathbf{e}_j)^2 =\sum_{i\in S}\sigma_i^2 . \end{aligned} \end{equation*}

So \(\|B\|_2^2=\max_{i\in S}\sigma_i^2\) and \(\|B\|_F^2=\sum_{i\in S}\sigma_i^2\). Take \(S=\{1,\ldots,k\}\) for \(A_k\) and \(S=\{k+1,\ldots,r\}\) for \(A-A_k\); the \(\sigma_i\) are nonincreasing, so the maxima are \(\sigma_1^2\) and \(\sigma_{k+1}^2\) (the last also being Lemma 3.8). \(\blacksquare\)

\(\mathbf{u}_i=\varepsilon_i\mathbf{v}_i\) with \(\varepsilon_i=\pm1\), and \(A=VDV^T\) with \(D=\operatorname{diag}(\varepsilon_1\sigma_1,\ldots,\varepsilon_r\sigma_r)\); the book’s claim holds only up to these signs (for \(A=\operatorname{diag}(1,-2)\), \(\mathbf{v}_1=\mathbf{e}_2\) but \(\mathbf{u}_1=A\mathbf{v}_1/2=-\mathbf{e}_2\)), with \(\mathbf{u}_i=\mathbf{v}_i\) for all \(i\) exactly when \(A\) is positive semidefinite.

By Lemma 3.10 and \(A^T=A\), we have \(A\mathbf{v}_i=\sigma_i\mathbf{u}_i\) and \(A\mathbf{u}_i=\sigma_i\mathbf{v}_i\). Hence

\begin{equation*} A(\mathbf{v}_i+\mathbf{u}_i)=\sigma_i(\mathbf{v}_i+\mathbf{u}_i),\qquad A(\mathbf{v}_i-\mathbf{u}_i)=-\sigma_i(\mathbf{v}_i-\mathbf{u}_i). \end{equation*}

Suppose both vectors were nonzero. They would be eigenvectors of the symmetric \(A\) for the distinct eigenvalues \(\pm\sigma_i\), hence orthogonal. Both would lie in \(E_i=\{\mathbf{x}:A^TA\mathbf{x}=A^2\mathbf{x}=\sigma_i^2\mathbf{x}\}\), so \(\dim E_i\ge2\). But by Exercise 3.12, \(A^TA=\sum_j\sigma_j^2\mathbf{v}_j\mathbf{v}_j^T\) with the \(\sigma_j\) distinct and positive, so \(E_i=\mathbb{R}\mathbf{v}_i\) is one-dimensional. Therefore \(\mathbf{u}_i=\mathbf{v}_i\) or \(\mathbf{u}_i=-\mathbf{v}_i\); write \(\mathbf{u}_i=\varepsilon_i\mathbf{v}_i\). Theorem 3.4 then gives, with \(V=[\mathbf{v}_1,\ldots,\mathbf{v}_r]\),

\begin{equation*} \begin{aligned} A&=\sum_{i=1}^r\sigma_i\mathbf{u}_i\mathbf{v}_i^T =\sum_{i=1}^r\varepsilon_i\sigma_i\mathbf{v}_i\mathbf{v}_i^T=VDV^T, \qquad D=\operatorname{diag}(\varepsilon_i\sigma_i). \end{aligned} \end{equation*}

Here \(\varepsilon_i\sigma_i=\mathbf{v}_i^TA\mathbf{v}_i\). If \(A\) is positive semidefinite this is nonnegative, so every \(\varepsilon_i=+1\) and the left and right singular vectors coincide. \(\blacksquare\)

Use the power method of Section 3.7: draw a random Gaussian \(\mathbf{x}\), normalize, and iterate \(\mathbf{x}\leftarrow A^TA\mathbf{x}/|A^TA\mathbf{x}|\), computing each step as two matrix-vector products (\(A\mathbf{x}\), then \(A^T(A\mathbf{x})\)), never forming \(A^TA\) (Section 3.7.1).

Write \(\mathbf{x}=\sum_{i=1}^d c_i\mathbf{v}_i\) in the right singular basis (completed to an orthonormal basis of \(\mathbb{R}^d\)). Since \(A^TA=\sum_i\sigma_i^2\mathbf{v}_i\mathbf{v}_i^T\),

\begin{equation*} \begin{aligned} (A^TA)^k\mathbf{x}&=\sum_{i}\sigma_i^{2k}c_i\mathbf{v}_i\\ &=\sigma_1^{2k}\Bigl(c_1\mathbf{v}_1 +\sum_{i\ge 2}\bigl(\tfrac{\sigma_i}{\sigma_1}\bigr)^{2k}c_i\mathbf{v}_i\Bigr), \end{aligned} \end{equation*}

so if \(\sigma_1>\sigma_2\) and \(c_1\ne0\) the normalized iterate tends to \(\pm\mathbf{v}_1\) geometrically at rate \((\sigma_2/\sigma_1)^2\). Without a gap, Theorem 3.11 still applies: if \(|\mathbf{x}^T\mathbf{v}_1|\ge\delta\), then after \(k=\ln(1/(\varepsilon\delta))/(2\varepsilon)\) steps the unit iterate \(\mathbf{w}\) has component at most \(\varepsilon\) orthogonal to the span \(V\) of singular vectors with \(\sigma_i>(1-\varepsilon)\sigma_1\), hence \(|A\mathbf{w}|^2\ge(1-\varepsilon^2)(1-\varepsilon)^2\sigma_1^2\); and by Lemma 3.12 a normalized Gaussian start has \(\delta\ge 1/(20\sqrt d)\) with probability at least \(9/10-3e^{-d/96}\). Then \(\sigma_1\approx|A\mathbf{w}|\) and \(\mathbf{u}_1=A\mathbf{v}_1/\sigma_1\).

For the next singular vectors, any of the following works.

(i) Deflation. By Theorem 3.4, \(A-\sigma_1\mathbf{u}_1\mathbf{v}_1^T=\sum_{i\ge2}\sigma_i\mathbf{u}_i\mathbf{v}_i^T\), whose first singular vector is \(\mathbf{v}_2\); run the power method on it (as \(\mathbf{x}\mapsto A\mathbf{x}-\sigma_1\mathbf{u}_1(\mathbf{v}_1^T\mathbf{x})\), keeping sparsity), and repeat.

(ii) Orthogonal restriction. By the greedy definition of Section 3.3, \(\mathbf{v}_2=\arg\max_{\mathbf{v}\perp\mathbf{v}_1,|\mathbf{v}|=1}|A\mathbf{v}|\); run the power method but project each iterate onto \(\{\mathbf{v}_1\}^\perp\). This is the same iteration as (i), since \(A^TA\) maps \(\mathbf{v}_1^\perp\) into itself.

(iii) Subspace (block) iteration, as in Section 3.7.1: start with \(k\) random vectors, and repeatedly multiply the current orthonormal basis by \(A^TA\) and re-orthonormalize (Gram-Schmidt/QR); the orthonormalization prevents the block collapsing onto \(\mathbf{v}_1\), and the span converges to \(\mathrm{span}(\mathbf{v}_1,\dots,\mathbf{v}_k)\) (when \(\sigma_k>\sigma_{k+1}\)), which by Theorem 3.1 is the best-fit \(k\)-dimensional subspace; a final SVD of the small \(n\times k\) matrix \(AQ\) (\(Q\) the current basis) gives the individual vectors.

\(\sigma_1=\sqrt{15+\sqrt{221}}\approx5.4650\), \(\sigma_2=\sqrt{15-\sqrt{221}}\approx0.3660\), and

\begin{equation*} \begin{aligned} A&=\sigma_1\begin{pmatrix}0.4046\\0.9145\end{pmatrix} \begin{pmatrix}0.5760&0.8174\end{pmatrix}\\ &\quad+\sigma_2\begin{pmatrix}0.9145\\-0.4046\end{pmatrix} \begin{pmatrix}-0.8174&0.5760\end{pmatrix}. \end{aligned} \end{equation*}

Iterate with \(B=A^TA=\begin{pmatrix}10&14\\14&20\end{pmatrix}\) from \(\mathbf{x}_0=(1,0)^T\):

\begin{equation*} \begin{aligned} B\mathbf{x}_0&=(10,\,14), & B^2\mathbf{x}_0&=(296,\,420),\\ B^3\mathbf{x}_0&=(8840,\,12544), & B^4\mathbf{x}_0&=(264016,\,374640), \end{aligned} \end{equation*}

with second-to-first ratios \(1.4,\ 1.41892,\ 1.419005,\ 1.4190049\); normalizing \(B^4\mathbf{x}_0\) gives \(\mathbf{v}_1=(0.57605,\,0.81742)^T\), already correct to 7 digits because the error contracts by \((\sigma_2/\sigma_1)^2\approx0.0045\) per step (Section 3.7). The limit is exact: \(\mathbf{v}_1\) is the eigenvector of \(B\) for its top eigenvalue \(\lambda_1=15+\sqrt{221}\) (roots of \(\lambda^2-30\lambda+4=0\)), namely \(\mathbf{v}_1\propto(14,\ 5+\sqrt{221})\), ratio \((5+\sqrt{221})/14=1.4190049\ldots\). Then

\begin{equation*} \sigma_1=|A\mathbf{v}_1|=\sqrt{\lambda_1}\approx5.46499,\qquad \mathbf{u}_1=\frac{A\mathbf{v}_1}{\sigma_1}=(0.40455,\,0.91451)^T . \end{equation*}

In \(\mathbb{R}^2\) the second right singular vector is forced to be the unit vector orthogonal to \(\mathbf{v}_1\) (Section 3.3), \(\mathbf{v}_2=(-0.81742,\,0.57605)^T\) (equivalently, power iteration on the deflated \(B-\sigma_1^2\mathbf{v}_1\mathbf{v}_1^T\) returns it in one step), so

\begin{equation*} \sigma_2=|A\mathbf{v}_2|=\sqrt{15-\sqrt{221}}\approx0.36597,\qquad \mathbf{u}_2=\frac{A\mathbf{v}_2}{\sigma_2}=(0.91451,\,-0.40455)^T, \end{equation*}

and \(A=\sigma_1\mathbf{u}_1\mathbf{v}_1^T+\sigma_2\mathbf{u}_2\mathbf{v}_2^T\) by Theorem 3.4.

\(\sigma_1=43.4304\) with

\begin{equation*} \mathbf{v}_1=(0.3198,\,0.3696,\,0.3981,\,0.4039,\,0.3873,\,0.3500,\,0.2951,\,0.2272,\,0.1514,\,0.0736)^T, \end{equation*}

and the next three singular values are \(23.9832,\ 14.1109,\ 10.4970\). Program (seed 0; each step applies \(A^TA\) as \(A^T(A\,\cdot)\), and np.linalg.qr plays the role of orth):

import numpy as np
rng = np.random.default_rng(0)
n = 10
A = np.array([[i+j+1 if i+j+1 <= n else 0 for j in range(n)] for i in range(n)], float)

def power(A, iters=100):
    x = rng.standard_normal(A.shape[1])
    for _ in range(iters):
        x = A.T @ (A @ x); x /= np.linalg.norm(x)
    return x, np.linalg.norm(A @ x)

def block_power(A, k=4, iters=100):
    Q, _ = np.linalg.qr(rng.standard_normal((A.shape[1], k)))
    for _ in range(iters):
        Q, _ = np.linalg.qr(A.T @ (A @ Q))
    return Q, np.linalg.norm(A @ Q, axis=0)

v1, s1 = power(A)
Q, s = block_power(A)

(1) power returns the \(\mathbf{v}_1\) and \(\sigma_1=|A\mathbf{v}_1|=43.43043\) above, and \(\mathbf{u}_1=A\mathbf{v}_1/\sigma_1=\mathbf{v}_1\). It agrees with a library SVD to \(4\times10^{-16}\) (up to sign), as expected from Section 3.7: the error contracts by \((\sigma_2/\sigma_1)^2=0.305\) per step, so 100 steps are far more than enough.

(2) block_power returns (columns, up to sign)

with \(|A\mathbf{v}_i|=43.4304,\ 23.9832,\ 14.1109,\ 10.4970\); each column matches the library singular vector to within \(2\times10^{-15}\). Re-orthonormalizing after each multiplication keeps the four vectors from collapsing onto \(\mathbf{v}_1\) (Section 3.7.1): the span converges to \(\mathrm{span}(\mathbf{v}_1,\dots,\mathbf{v}_4)\) at rate \((\sigma_5/\sigma_4)^2=(8.392/10.497)^2=0.639\) per step, and because the QR basis is nested (column \(j\) spans the same space as the first \(j\) columns of the iterated block), column \(j\) converges to \(\mathbf{v}_j\) individually. Since \(A\) is symmetric, \(\mathbf{u}_i=A\mathbf{v}_i/\sigma_i=\pm\mathbf{v}_i\); here \(\mathbf{u}_i=\mathbf{v}_i\) for \(i=1,3\) and \(\mathbf{u}_i=-\mathbf{v}_i\) for \(i=2,4\), because \(A\) has eigenvalues \(43.43,\,-23.98,\,14.11,\,-10.50\). Multiplying by \(A\) as the exercise says, instead of \(A^TA\), gives the same subspaces, since two applications of \(A\) are one application of \(A^TA=A^2\).

(1) For every \(\mathbf{x}\),

\begin{equation*} \mathbf{x}^TAA^T\mathbf{x}=(A^T\mathbf{x})^T(A^T\mathbf{x})=|A^T\mathbf{x}|^2\ge0 . \end{equation*}

(2) Take the graph simple and undirected (so \(a_{ij}=a_{ji}\in\{0,1\}\), \(a_{ii}=0\), and \(d_{ii}=\deg i\)). Each row of \(B\) has exactly two nonzero entries, \(\pm1\), in the columns of the two endpoints of \(e\). Hence

\begin{equation*} (B^TB)_{ij}=\sum_{e}b_{ei}b_{ej}= \begin{cases} \sum_e b_{ei}^2=\#\{e: i\in e\}=\deg i=d_{ii}, & i=j,\\ -1 \ \text{ if } \{i,j\}\in E,\ \ 0 \text{ otherwise}, & i\ne j, \end{cases} \end{equation*}

since for \(i\ne j\) the only edge with \(b_{ei}b_{ej}\ne0\) is \(e=\{i,j\}\), where one endpoint carries \(-1\) and the other \(+1\). So \((B^TB)_{ij}=d_{ij}-a_{ij}\) for all \(i,j\), i.e. \(L=B^TB\), and by part 1 (applied to the real matrix \(B^T\))

\begin{equation*} \mathbf{x}^TL\mathbf{x}=|B\mathbf{x}|^2=\sum_{\{i,j\}\in E}(x_i-x_j)^2\ge0 . \end{equation*}

Let \(A=A^T\) be real, \(A\mathbf{v}=\lambda\mathbf{v}\) with \(\mathbf{v}\in\mathbb{C}^n\setminus\{0\}\), \(\lambda\in\mathbb{C}\), and put \(q=\bar{\mathbf{v}}^TA\mathbf{v}\). Then \(q=\lambda\,\bar{\mathbf{v}}^T\mathbf{v}=\lambda|\mathbf{v}|^2\), while, using \(\bar A=A\) and that a \(1\times1\) matrix equals its transpose,

\begin{equation*} \begin{aligned} \bar q&=\mathbf{v}^TA\bar{\mathbf{v}}=(\mathbf{v}^TA\bar{\mathbf{v}})^T\\ &=\bar{\mathbf{v}}^TA^T\mathbf{v}=\bar{\mathbf{v}}^TA\mathbf{v}=q . \end{aligned} \end{equation*}

So \(q\in\mathbb{R}\), and since \(|\mathbf{v}|^2=\sum_i|v_i|^2>0\), \(\lambda=q/|\mathbf{v}|^2\in\mathbb{R}\).

Let \(C=\sum_i\frac1{\sigma_i}\mathbf{v}_i\mathbf{u}_i^T\) (well defined, as all singular values are positive). An invertible \(n\times n\) matrix has rank \(n\), so there are \(r=n\) terms and \(\mathbf{v}_1,\dots,\mathbf{v}_n\) is an orthonormal basis of \(\mathbb{R}^n\); the \(\mathbf{u}_i\) are orthonormal by Theorem 3.7. Hence

\begin{equation*} \begin{aligned} CA&=\sum_{i,j}\frac{\sigma_j}{\sigma_i} \mathbf{v}_i(\mathbf{u}_i^T\mathbf{u}_j)\mathbf{v}_j^T\\ &=\sum_{i=1}^n\mathbf{v}_i\mathbf{v}_i^T=VV^T=I, \end{aligned} \end{equation*}

where \(V=(\mathbf{v}_1\cdots\mathbf{v}_n)\) is orthogonal, so \(V^TV=I\) forces \(VV^T=I\). A left inverse of an invertible matrix is its inverse: \(C=CAA^{-1}=A^{-1}\).

\(BA=\sum_{i=1}^r\mathbf{v}_i\mathbf{v}_i^T\), the orthogonal projection onto \(\mathrm{span}(\mathbf{v}_1,\dots,\mathbf{v}_r)\). Indeed, by orthonormality of the \(\mathbf{u}_i\) (Theorem 3.7),

\begin{equation*} \begin{aligned} BA&=\sum_{i,j=1}^r\frac{\sigma_j}{\sigma_i} \mathbf{v}_i(\mathbf{u}_i^T\mathbf{u}_j)\mathbf{v}_j^T\\ &=\sum_{i=1}^r\mathbf{v}_i\mathbf{v}_i^T, \end{aligned} \end{equation*}

so for \(\mathbf{x}=\sum_{j=1}^r c_j\mathbf{v}_j\), by orthonormality of the \(\mathbf{v}_i\),

\begin{equation*} BA\mathbf{x}=\sum_{i=1}^r\mathbf{v}_i\sum_{j=1}^r c_j\,\mathbf{v}_i^T\mathbf{v}_j =\sum_{i=1}^r c_i\mathbf{v}_i=\mathbf{x}. \end{equation*}

(1) Since \(\sigma_1\ge\sigma_2\ge\cdots\), Lemma 3.2 gives

\begin{equation*} k\sigma_k^2 \le \sigma_1^2+\cdots+\sigma_k^2 \le \sum_{i}\sigma_i^2 = \|A\|_F^2 , \end{equation*}

so \(\sigma_k\le \|A\|_F/\sqrt{k}\) (with \(\sigma_k=0\) if \(k>\operatorname{rank}A\)).

(2) Take \(B=A_k\) (or \(B=A\) if \(\operatorname{rank}A\le k\)). By Lemma 3.8 and part 1 applied to index \(k+1\),

\begin{equation*} \|A-A_k\|_2=\sigma_{k+1}\le \frac{\|A\|_F}{\sqrt{k+1}}\le \frac{\|A\|_F}{\sqrt{k}} . \end{equation*}

(3) No, except for \(k=1\). For \(k=1\), \(\|A-A_1\|_F^2=\|A\|_F^2-\sigma_1^2\le\|A\|_F^2\). For \(k\ge 2\) take \(A=I_{2k+1}\), all of whose \(2k+1\) singular values equal \(1\). By Theorem 3.6, every \(B\) of rank at most \(k\) has

\begin{equation*} \|A-B\|_F^2\ge\|A-A_k\|_F^2=\sum_{i>k}\sigma_i^2=k+1>2+\frac1k=\frac{\|A\|_F^2}{k}, \end{equation*}

because \(k^2-k-1>0\) for \(k\ge2\).

Preprocess by computing the SVD and storing \(A_k=\sum_{i=1}^k\sigma_i\mathbf{u}_i\mathbf{v}_i^T\) in factored form, with \(k=\min(\lceil 1/\varepsilon^2\rceil,\operatorname{rank}A)\): keep the \(d\times k\) matrix \(V_k=[\mathbf{v}_1,\ldots,\mathbf{v}_k]\) and the \(n\times k\) matrix \(W_k=[\sigma_1\mathbf{u}_1,\ldots,\sigma_k\mathbf{u}_k]\). For each query output

\begin{equation*} \mathbf{y}_j=W_k\,(V_k^T\mathbf{x}_j)=A_k\mathbf{x}_j . \end{equation*}

Computing \(\mathbf{z}=V_k^T\mathbf{x}_j\in\mathbb{R}^k\) takes \(O(dk)\) operations and \(W_k\mathbf{z}\) takes \(O(nk)\), so with \(k\le 1/\varepsilon^2+1\) the cost is \(O\big((d+n)/\varepsilon^2\big)\) per query (for \(\varepsilon\le1\); for \(\varepsilon>1\) the output \(\mathbf{y}_j=0\) already works, since \(|A\mathbf{x}_j|\le\|A\|_2|\mathbf{x}_j|\le\|A\|_F|\mathbf{x}_j|\)).

If \(k=\operatorname{rank}A\) then \(A_k=A\) and \(\mathbf{y}_j=A\mathbf{x}_j\) exactly. Otherwise \(k=\lceil1/\varepsilon^2\rceil\), and by Lemma 3.8 and Exercise 3.22(1) applied to index \(k+1\),

\begin{equation*} \begin{aligned} |\mathbf{y}_j-A\mathbf{x}_j|=|(A-A_k)\mathbf{x}_j| &\le\|A-A_k\|_2\,|\mathbf{x}_j|=\sigma_{k+1}|\mathbf{x}_j|\\ &\le\frac{\|A\|_F}{\sqrt{k+1}}\,|\mathbf{x}_j|\le\varepsilon\|A\|_F|\mathbf{x}_j| , \end{aligned} \end{equation*}

since \(k+1>1/\varepsilon^2\).

Take \(c_1=\pm1\) and \(c_i=0\) for \(i\ge2\); the maximum is \(\sigma_1^2\). Indeed, for any feasible \(c\),

\begin{equation*} \sum_{i=1}^r c_i^2\sigma_i^2\le\sigma_1^2\sum_{i=1}^r c_i^2=\sigma_1^2 , \end{equation*}

with equality if and only if \(c_i=0\) for every \(i\) with \(\sigma_i<\sigma_1\). Thus the maximizers are exactly the unit vectors \(c\) supported on \(\{i:\sigma_i=\sigma_1\}\); when \(\sigma_1>\sigma_2\) these are only \(c=\pm\mathbf{e}_1\).

(1) It is the first right-singular vector \(\mathbf{v}_1\) of \(A\) (synthetic documents are taken of unit length, else the objective is unbounded). For a unit \(\mathbf{x}\) with rows \(\mathbf{a}_1,\ldots,\mathbf{a}_m\) of \(A\),

\begin{equation*} \sum_{i=1}^m(\mathbf{a}_i\cdot\mathbf{x})^2=|A\mathbf{x}|^2 , \end{equation*}

which by definition (Section 3.3) is maximized by \(\mathbf{x}=\mathbf{v}_1\), with value \(\sigma_1^2\). Find it by the power method (Section 3.7, Theorem 3.11) or any SVD routine.

(2) The center of gravity \(\mathbf{c}=\frac1m\sum_i\mathbf{a}_i\) maximizes, in direction, the sum of the (signed) similarities, since \(\sum_i\mathbf{a}_i\cdot\mathbf{x}=m\,\mathbf{c}\cdot\mathbf{x}\) is maximized over unit \(\mathbf{x}\) by \(\mathbf{c}/|\mathbf{c}|\); \(\mathbf{v}_1\) maximizes the sum of squared similarities instead, and is the best-fit line through the origin rather than a mean (it is the best-fit line through \(\mathbf{c}\) only after centering, Section 3.9.1). They differ in general: for the three documents \(\mathbf{e}_1,\mathbf{e}_1,\mathbf{e}_2\) we get \(\mathbf{v}_1=\mathbf{e}_1\) (since \(A^TA=\operatorname{diag}(2,1)\)), with value \(2\), whereas \(\mathbf{c}=(\tfrac23,\tfrac13)\) and \(\mathbf{c}/|\mathbf{c}|=(2,1)/\sqrt5\) scores only \(\tfrac45+\tfrac45+\tfrac15=\tfrac95\).

(3) Take \(\mathbf{v}_1,\ldots,\mathbf{v}_k\), the top \(k\) right-singular vectors. For orthonormal \(\mathbf{x}_1,\ldots,\mathbf{x}_k\) spanning a subspace \(V\),

\begin{equation*} \begin{aligned} \sum_{j=1}^k\sum_{i=1}^m(\mathbf{a}_i\cdot\mathbf{x}_j)^2 &=\sum_{i=1}^m|\operatorname{proj}_V\mathbf{a}_i|^2 , \end{aligned} \end{equation*}

and by Theorem 3.1 this is maximized by \(V=V_k=\operatorname{span}(\mathbf{v}_1,\ldots,\mathbf{v}_k)\), with value \(\sigma_1^2+\cdots+\sigma_k^2\). Any orthonormal basis of \(V_k\) is optimal; the singular vectors are the canonical one.

(4) Permuting rows and columns only relabels coordinates, so assume \(A=\operatorname{diag}(A_1,\ldots,A_k)\) with \(A_i\) the \(i\)-th block. Then

\begin{equation*} A^TA=\operatorname{diag}(A_1^TA_1,\ldots,A_k^TA_k). \end{equation*}

A right-singular vector \(\mathbf{v}\) with singular value \(\sigma>0\) satisfies \(A^TA\mathbf{v}=\sigma^2\mathbf{v}\) (Section 3.8). Writing \(\mathbf{v}=(\mathbf{w}_1,\ldots,\mathbf{w}_k)\) in blocks, the block-diagonal form splits this into

\begin{equation*} A_i^TA_i\,\mathbf{w}_i=\sigma^2\mathbf{w}_i\qquad(i=1,\ldots,k). \end{equation*}

So every nonzero block \(\mathbf{w}_i\) is an eigenvector of \(A_i^TA_i\) with eigenvalue \(\sigma^2\), which means \(\mathbf{w}_i/|\mathbf{w}_i|\) is a right-singular vector of \(A_i\) with singular value \(\sigma\) (the right-singular vectors of \(A_i\) for \(\sigma\) are exactly the unit vectors of this eigenspace). Blocks that are zero are the degenerate case, and \(\mathbf{u}=A\mathbf{v}/\sigma\) splits the same way.

(5) Read “distinct” as all nonzero singular values of \(A\) being distinct, and “similar” as: within a cluster any two documents are joined by a chain of documents, consecutive ones sharing a term. Then:

(i) Each singular pair \((\mathbf{u}_j,\mathbf{v}_j)\) of \(A\) lives in a single block. If \(\mathbf{v}_j\) had nonzero blocks \(\mathbf{w}_i,\mathbf{w}_l\) with \(i\ne l\), then by (4) the vectors \((0,\ldots,\mathbf{w}_i,\ldots,0)\) and \((0,\ldots,\mathbf{w}_l,\ldots,0)\) would be independent eigenvectors of \(A^TA\) for \(\sigma_j^2\), so \(\sigma_j\) would be a repeated singular value. Hence \(\mathbf{v}_j\) is supported on one term cluster, and \(\mathbf{u}_j=A\mathbf{v}_j/\sigma_j\) on the matching document cluster.

(ii) Each cluster is covered by one support. \(A_iA_i^T\ge0\) entrywise, and it is irreducible by the chain assumption, so by Perron-Frobenius its top eigenvalue \(\sigma_1(A_i)^2\) is simple, with an eigenvector \(\mathbf{p}_i\) that is strictly positive on all documents of cluster \(i\). Because \(\sigma_1(A_i)\) is a simple singular value of \(A\), the corresponding \(\mathbf{u}_j\) is \(\pm(0,\ldots,\mathbf{p}_i,\ldots,0)\), so its support is exactly cluster \(i\).

Algorithm: compute the SVD \(A=\sum_{j=1}^r\sigma_j\mathbf{u}_j\mathbf{v}_j^T\) and merge documents that lie in a common support \(\operatorname{supp}(\mathbf{u}_j)\) (union-find). By (i) no merge crosses clusters, and by (ii) each cluster gets fully merged. So the classes are exactly the \(k\) document clusters. Each term goes to the cluster of any document containing it (equivalently, by the supports of the \(\mathbf{v}_j\)). When the \(k\) largest singular values are the block tops \(\sigma_1(A_i)\), which is the typical case of strong within-cluster similarity, the \(k\) supports of \(\mathbf{u}_1,\ldots,\mathbf{u}_k\) already are the clusters.

With \(S=\{i:x_i=1\}\), \(s_1=\sum_{i\in S}u_i=\mathbf{u}^T\mathbf{x}\), \(s_2=\sum_{i\notin S}u_i=\mathbf{u}^T(\mathbf{1}-\mathbf{x})\) and \(T=s_1+s_2=\sum_iu_i\) (fixed),

\begin{equation*} \begin{aligned} \mathbf{x}^T\mathbf{u}\mathbf{u}^T(\mathbf{1}-\mathbf{x}) &=(\mathbf{u}^T\mathbf{x})\big(\mathbf{u}^T(\mathbf{1}-\mathbf{x})\big)\\ &=s_1s_2=\frac{T^2-(s_1-s_2)^2}{4}. \end{aligned} \end{equation*}

Since \(x\leftrightarrow S\) is a bijection between \(\{0,1\}^n\) and partitions \((S,\bar S)\) of the coordinates, maximizing the objective is the same as minimizing \(|s_1-s_2|\), i.e. making the two subset sums as equal as possible.

At 5% of the singular values the photo is already recognizable but blurred and streaky. It looks essentially perfect from 25% on, and in every case over 97.5% of \(\|A\|_F^2\) is captured. I used the \(600\times512\) colour photo grace_hopper.jpg bundled with matplotlib. Each RGB channel \(A^{( c)}\) was converted to double and replaced by its best rank-\(k\) approximation \(A^{( c)}_k\) (Theorem 3.6), with \(k\) equal to 5%, 10%, 25%, 50% of the \(r=512\) singular values.

import numpy as np, matplotlib; import matplotlib.pyplot as plt
A = plt.imread(matplotlib.cbook.get_sample_data('grace_hopper.jpg')).astype(float)
r = min(A.shape[:2])
svd = [np.linalg.svd(A[:, :, c], full_matrices=False) for c in range(3)]
tot = sum((s**2).sum() for _, s, _ in svd)
for p in [5, 10, 25, 50]:
    k = round(p * r / 100)
    R = np.stack([U[:, :k] * s[:k] @ Vt[:k] for U, s, Vt in svd], axis=2)
    en = sum((s[:k]**2).sum() for _, s, _ in svd) / tot
    print(p, k, en, np.sqrt(en), np.linalg.norm(A - R) / np.linalg.norm(A))
    plt.imsave('recon%d.jpg' % p, np.clip(R, 0, 255).astype(np.uint8))

(1) At \(k=26\) (5%) the face, cap and flag are clearly recognizable, but fine texture is blurred and there are visible horizontal/vertical streaks. These are typical rank-\(k\) artifacts, since each term \(\sigma_i\mathbf{u}_i\mathbf{v}_i^T\) is an outer product. At \(k=51\) (10%) the streaks are faint and the image is good. At \(k=128\) and \(k=256\) the reconstructions cannot be told apart from the original by eye.

(2) Pooling the three channels, the captured fraction \(\sum_{i\le k}\sigma_i^2/\sum_i\sigma_i^2=\|A_k\|_F^2/\|A\|_F^2\) (Lemma 3.2) and the relative error are as follows; the ratio of norms \(\|A_k\|_F/\|A\|_F\) is the square root of the first column.

The fraction is so high even at 5% because a nonnegative photo has a huge \(\sigma_1\) (the mean brightness pattern). Visual quality is tracked better by the relative error \(\sqrt{\sum_{i>k}\sigma_i^2}/\|A\|_F\), which only drops below about 5% at \(k=128\).

Correlated: the top 10 singular vectors capture 99.2% of \(\|A\|_F^2\) and 3 already reach 95%. Independent: the top 10 capture 82.5% and 95% needs 40. Here “fraction captured by the top \(k\)” is \(\sum_{i\le k}\sigma_i^2/\|A\|_F^2=\|A_k\|_F^2/\|A\|_F^2\) (Lemma 3.2). Seed 0 was used, and the correlated matrix was made by replacing each entry of an i.i.d. uniform matrix by the average of its \(9\times9\) neighbourhood and then rescaling to \([0,1]\).

import numpy as np
from scipy.ndimage import uniform_filter
rng = np.random.default_rng(0)
X = rng.random((100, 100))                       # independent U(0,1)
C = uniform_filter(X, size=9, mode='reflect')    # 9x9 local averages
C = (C - C.min()) / (C.max() - C.min())          # rescale to [0,1]
for name, A in [('correlated', C), ('independent', X)]:
    s = np.linalg.svd(A, compute_uv=False)
    e = np.cumsum(s**2) / (s**2).sum()
    adj = np.corrcoef(A[:, :-1].ravel(), A[:, 1:].ravel())[0, 1]
    print(name, adj, e[0], e[9], np.searchsorted(e, 0.95) + 1)

If the fraction is measured instead by \(\|A_k\|_F/\|A\|_F\), the top 10 give 99.6% and 90.9%, and 95% needs 1 and 25 vectors.

(1) Smoothing removes the high-frequency content, so \(A\) is close to a low-rank matrix: one singular vector (essentially the mean level) carries 92% of the energy and the next two bring it past 95%. By Theorem 3.6, \(A_{10}\) is then a near-perfect approximation, with \(\|A-A_{10}\|_F/\|A\|_F=\sqrt{1-0.992}\approx0.09\).

(2) Write \(A=\tfrac12 J+E\), where \(J\) is the all-ones matrix and \(E\) has i.i.d. mean-zero entries of variance \(\tfrac1{12}\). The rank-one part gives \(\sigma_1^2\approx\|\tfrac12J\|_F^2=n^2/4\), out of \(\mathbb{E}\|A\|_F^2=n^2/3\), i.e. 75%, matching the observed 75.4%. The remaining \(n^2/12\) is noise spread over all 100 directions with no preferred low-dimensional subspace. Its singular values decay slowly, so reaching 95% requires capturing 80% of the noise energy, which takes 40 singular vectors.

Compute \(A_k\) in \(O(n^3)\) time, then replace the dynamic program by \((k+1)^k\le e\,k^k\) linear programs of polynomial size, one for each guess of the vector \((\mathbf u_1^T\mathbf x,\dots,\mathbf u_k^T\mathbf x)\) on a grid of mesh \(2\sqrt n/k\). (The book’s list of \(\tilde w(S)\) can hold \((2n^{3/2}k^2)^{2k}=n^{\Theta(k)}\) vectors, so the second phase has to change.) We may assume \(k\le n\) (otherwise enumerate all \(2^n\le k^k\) cuts) and \(a_{ii}=0\) (self-loops are never cut). Write

\begin{equation*} \begin{aligned} f(\mathbf x)&=\mathbf x^TA(\mathbf 1-\mathbf x)=\sum_{i\ne j}a_{ij}x_i(1-x_j),\\ g(\mathbf x)&=\mathbf x^TA_k(\mathbf 1-\mathbf x)=\sum_{i=1}^k\sigma_i\,\alpha_i(\mathbf x)\beta_i(\mathbf x), \end{aligned} \end{equation*}

where \(\alpha_i(\mathbf x)=\mathbf u_i^T\mathbf x\) and \(\beta_i(\mathbf x)=\mathbf v_i^T(\mathbf 1-\mathbf x)\). For \(\mathbf x\in[0,1]^n\) we have \(|\mathbf x|,|\mathbf 1-\mathbf x|\le\sqrt n\), so \(|\alpha_i|,|\beta_i|\le\sqrt n\), and as in the proof of Item 1 (via Lemma 3.8), \(|f(\mathbf x)-g(\mathbf x)|\le n\sigma_{k+1}\le n^2/\sqrt{k+1}\). Also \(\sum_i\sigma_i\le\sqrt k\,\|A\|_F\le\sqrt k\,n\) by Cauchy-Schwarz.

Algorithm. Let \(T=\{-\sqrt n+2j\sqrt n/k: j=0,\dots,k\}\), so every \(t\in[-\sqrt n,\sqrt n]\) is within \(\sqrt n/k\) of \(T\). For each \(\mathbf a\in T^k\): (i) solve the LP: maximize \(h_{\mathbf a}(\mathbf x)=\sum_i\sigma_ia_i\beta_i(\mathbf x)\) over \(\mathbf x\in[0,1]^n\) with \(|\alpha_i(\mathbf x)-a_i|\le\sqrt n/k\) for all \(i\) (skip if infeasible), getting \(\hat{\mathbf x}\); (ii) round \(\hat{\mathbf x}\) to \(\mathbf x^{\prime}\in\{0,1\}^n\) one coordinate at a time, each time choosing the endpoint that does not decrease \(f\). Since \(a_{ii}=0\), \(f\) is affine in each coordinate, so \(f(\mathbf x^{\prime})\ge f(\hat{\mathbf x})\). Output the best \(\mathbf x^{\prime}\) found.

Accuracy. Let \(\mathbf x^*\) be a maximum cut and pick \(\mathbf a\in T^k\) with \(|\alpha_i(\mathbf x^*)-a_i|\le\sqrt n/k\); then \(\mathbf x^*\) is LP-feasible. For every feasible \(\mathbf x\),

\begin{equation*} |g(\mathbf x)-h_{\mathbf a}(\mathbf x)|\le\sum_i\sigma_i|\alpha_i-a_i||\beta_i|\le\frac{\sqrt n}{k}\sqrt n\cdot\sqrt k\,n=\frac{n^2}{\sqrt k}. \end{equation*}

Hence

\begin{equation*} \begin{aligned} f(\mathbf x^{\prime})\ge f(\hat{\mathbf x})&\ge g(\hat{\mathbf x})-\tfrac{n^2}{\sqrt{k+1}}\ge h_{\mathbf a}(\hat{\mathbf x})-\tfrac{2n^2}{\sqrt k}\\ &\ge h_{\mathbf a}(\mathbf x^*)-\tfrac{2n^2}{\sqrt k}\ge g(\mathbf x^*)-\tfrac{3n^2}{\sqrt k}\ge f(\mathbf x^*)-\tfrac{4n^2}{\sqrt k}, \end{aligned} \end{equation*}

which is the guarantee of Theorem 3.19.

Time. The top \(k\) singular triples cost \(O(n^3)\). Each of the \((k+1)^k=k^k(1+1/k)^k\le e\,k^k\) guesses costs one LP with \(n\) variables and \(2n+2k\le 4n\) constraints, solvable in \(\mathrm{poly}(n)\) time (interior point), plus \(O(n^2)\) for the rounding (keep the \(n\) partial derivatives of \(f\) and update them in \(O(n)\) per fixed coordinate). Total: \(O(n^3+\mathrm{poly}(n)\,k^k)\). \(\blacksquare\)

(1) Put \(s=\sum_k|\mathbf x_k|^2\). Expanding \(d_{ij}^2=|\mathbf x_i|^2+|\mathbf x_j|^2-2\mathbf x_i\mathbf x_j^T\) and using \(\sum_k\mathbf x_k=0\) to kill the cross terms in the sums,

\begin{equation*} \begin{aligned} \frac1n\sum_k d_{ik}^2&=|\mathbf x_i|^2+\frac sn,\qquad \frac1n\sum_k d_{kj}^2=|\mathbf x_j|^2+\frac sn,\qquad \frac1{n^2}\sum_{k,l}d_{kl}^2=\frac{2s}n, \end{aligned} \end{equation*}

so the bracket equals \(|\mathbf x_i|^2+|\mathbf x_j|^2-2\mathbf x_i\mathbf x_j^T-|\mathbf x_i|^2-|\mathbf x_j|^2-\frac{2s}n+\frac{2s}n=-2\mathbf x_i\mathbf x_j^T\), as claimed. In matrix form, with \(J=I-\frac1n\mathbf 1\mathbf 1^T\) and \(D^{(2)}=(d_{ij}^2)\),

\begin{equation*} B:=XX^T=-\tfrac12\,J D^{(2)} J . \end{equation*}

(2) (Classical multidimensional scaling.) (i) Form \(B=-\frac12 JD^{(2)}J\) from part 1 (\(O(n^2)\) work: subtract row and column means of \(D^{(2)}\), add back the grand mean, multiply by \(-\frac12\)). (ii) \(B=XX^T\) is symmetric positive semidefinite of rank at most \(d\), so compute its spectral decomposition (for a PSD matrix this is its SVD, Theorem 3.4) \(B=\sum_{t=1}^{d}\lambda_t\mathbf w_t\mathbf w_t^T\) with \(\lambda_1\ge\dots\ge\lambda_d\ge0\). (iii) Output \(\hat X=[\sqrt{\lambda_1}\mathbf w_1,\dots,\sqrt{\lambda_d}\mathbf w_d]\), an \(n\times d\) matrix with \(\hat X\hat X^T=B\). Then \(\hat X=XQ\) for an orthogonal \(d\times d\) matrix \(Q\): from \(\hat X\hat X^T=XX^T\) we get \(|\hat X^T\mathbf y|=|X^T\mathbf y|\) for all \(\mathbf y\), so \(X^T\mathbf y\mapsto\hat X^T\mathbf y\) is a well-defined isometry between the row spaces, which extends to an orthogonal \(Q^T\) on \(\mathbf R^d\). Thus \(\hat X\) recovers the points up to exactly the rotation/reflection ambiguity in the statement. If \(D\) is only approximately Euclidean, \(B\) may have small negative eigenvalues; use the top \(d\) eigenvalues, clipped at \(0\). \(\blacksquare\)

(1) Both maps come out recognisable: the US one is a \(180^\circ\) rotation of the true map and the Chinese one a mirror image. The program implements Exercise 3.30, part 2, with \(d=2\):

import numpy as np
def cmdscale(D, d=2):
    n = len(D)
    J = np.eye(n) - np.ones((n, n)) / n
    B = -0.5 * J @ (D**2) @ J          # B = X X^T, Exercise 3.30(1)
    w, V = np.linalg.eigh(B)
    idx = np.argsort(w)[::-1]
    w, V = w[idx], V[:, idx]
    return V[:, :d] * np.sqrt(w[:d]), w
# D_us: 20x20 table above (diagonal 0); D_cn: 12x12 table
X_us, w_us = cmdscale(D_us)
X_cn, w_cn = cmdscale(D_cn)

In the figure, the US coordinates are plotted as \((-x_1,-x_2)\) and the Chinese ones as \((-x_1,x_2)\), so that north is up and west is left.

(i) US cities. The eigenvalues of \(B\) are \(\lambda_1=1.460\times10^7\), \(\lambda_2=2.797\times10^6\), then \(\lambda_3=2.14\times10^4\), \(\lambda_4=1.85\times10^4\), down to \(\lambda_{20}=-4.93\times10^4\). The top two carry \(99.7\%\) of the positive eigenvalue mass, and the negative eigenvalues total only \(0.46\%\) of it. So the data are almost exactly a planar Euclidean configuration (Exercise 3.30). The map has Seattle, San Francisco, Los Angeles and Phoenix on the west, Boston, New York, Philadelphia and Washington in the north-east, and Miami in the south-east. Distances recomputed from the map match the table with median relative error \(0.10\%\) (worst \(9.2\%\)). The leftover non-Euclidean part comes from inconsistent entries (for example, \(d(\mathrm{NY},\mathrm{DC})=230>83+123=d(\mathrm{NY},\mathrm{PHI})+d(\mathrm{PHI},\mathrm{DC})\)) and from the Earth’s curvature. (ii) Chinese cities. The eigenvalues are \(2.49\times10^7\), \(1.06\times10^7\), \(2.99\times10^6\), \(2.21\times10^6,\dots\), down to \(-2.53\times10^6\). The top two carry only \(86.7\%\) of the positive mass, and the negative eigenvalues total \(8.1\%\), so this table is far from Euclidean. The main culprit is the entry \(d(\mathrm{Beijing},\mathrm{Chongqing})=3026\), which violates the triangle inequality by \(947\) km, since \(d(\mathrm{Beijing},\mathrm{Tianjin})+d(\mathrm{Tianjin},\mathrm{Chongqing})=125+1954=2079\). This is almost certainly a misprint, and \(28\) ordered triples violate the triangle inequality in all. Even so, the mirrored map is geographically right. Harbin, Changchun and Shenyang lie in the north-east. Beijing and Tianjin are next to each other, with Hohhot and Yinchuan to their west. Shanghai is on the east coast, Chongqing and Nanning are in the south-west and south, and Lhasa and Urumqi are far to the west. The median relative distance error is \(4.2\%\).

(2) Yes, but not by feeding the airline distances in directly. An airline distance is a great-circle arc \(s_{ij}=R\theta_{ij}\) on the sphere of radius \(R\approx3959\) miles, not a straight-line distance in \(\mathbf R^3\). First convert each arc to the chord

\begin{equation*} c_{ij}=2R\sin\frac{s_{ij}}{2R}=|\mathbf p_i-\mathbf p_j| , \end{equation*}

which is the Euclidean distance between the cities’ positions \(\mathbf p_i\in\mathbf R^3\). Then apply Exercise 3.30 with \(d=3\). The matrix \(B\) has rank \(3\), and its top three eigenvectors give the \(\mathbf p_i\) (shifted to their centroid) up to a rotation or reflection. The cities then lie on a sphere of radius \(R\), which is the globe. Feeding in the raw arcs instead gives a non-Euclidean \(D\) (arcs exceed chords), so \(B\) generally has rank above \(3\) and negative eigenvalues and the model is distorted.

Wherever the ten nearest neighbours of a point are its neighbours along the curve, there is one large singular value and the test correctly reports dimension \(1\). The book’s curve, however, folds back on itself, and where neighbours from the other branch are close the test reports dimension \(2\).

import numpy as np
i = np.arange(1, 101)
x = np.sin(np.pi / 100 * i)
data = np.c_[x, np.sqrt(1 - x**2), 0.003 * i]   # book's create_sheath, n=100

def local_svals(data, k=10):
    out = []
    for p in range(len(data)):
        d = np.linalg.norm(data - data[p], axis=1)
        nb = np.argsort(d)[1:k + 1]             # k nearest neighbours
        M = data[nb] - data[p]                  # k x 3 matrix of difference vectors
        out.append(np.linalg.svd(M, compute_uv=False))
    return np.array(out)

S = local_svals(data)
print(np.median(S[5:23], 0), np.median(S[23:46], 0), S[49])

(i) The book’s own matrix (point \(6\) against points \(1\) to \(5\) and \(7\) to \(11\)) has singular values \(0.330,\ 0.0218,\ 3.9\times10^{-5}\). That is one large value, \(\sigma_2/\sigma_1=0.066\) coming from the curvature of the arc, and \(\sigma_3\) essentially zero. (ii) Since \(y=\sqrt{1-\sin^2t}=|\cos t|\) with \(t=\pi i/100\), the \((x,y)\) projection runs along the quarter circle up to \(i=50\) and then retraces it backwards, while \(z\) keeps increasing. So points \(i\) and \(100-i\) differ only by \(0.003|100-2i|\) in \(z\). Neighbours along the curve are \(\approx0.031\) apart, so for \(24\le i\le46\) and \(54\le i\le76\) two to five of the ten nearest neighbours come from the other branch. (iii) The table gives the medians over the groups (the curve’s two ends are left out).

In the middle group \(\sigma_2/\sigma_1\) reaches \(0.80\). The test cannot tell a curve that doubles back within the neighbourhood radius from a genuine surface. (iv) On the unfolded curve \((\sin t,\cos t,0.003i)\), every point has the profile \((0.330,\,0.0218,\,4\times10^{-5})\). The curvature term \(\sigma_2/\sigma_1\) grows with the neighbourhood size: it is \(0.029\), \(0.066\) and \(0.127\) for \(4\), \(10\) and \(20\) neighbours. (v) Adding Gaussian noise of standard deviation \(0.001\) and \(0.01\) to the unfolded data (numpy default_rng(0)) gives median singular values \((0.330,\,0.022,\,0.0030)\) and \((0.334,\,0.051,\,0.024)\): noise lifts the small singular values to roughly the noise level times \(\sqrt{10}\). \(\blacksquare\)

No. Both examples are connected, so Theorem 4.2 gives \(\mathbf{a}(t)\to\pi\), but \(\mathbf{p}(t)\) cycles.

(i) \(P=\begin{pmatrix}0&1\\1&0\end{pmatrix}\) and \(\mathbf{p}(0)=(1,0)\) give

\begin{equation*} \mathbf{p}(t)=\begin{cases}(1,0), & t \text{ even},\\ (0,1), & t \text{ odd},\end{cases} \qquad \mathbf{a}(t)\to\pi=\left(\tfrac12,\tfrac12\right), \end{equation*}

so \(\mathbf{p}(t)\) has no limit.

(ii) Here \(\mathbf{p}(0)=(1,0,0)\) gives \(\mathbf{p}(1)=(0,1,0)\), \(\mathbf{p}(2)=(0,0,1)\), \(\mathbf{p}(3)=(1,0,0)\), and so on with period 3, while \(\mathbf{a}(t)\to\pi=(\tfrac13,\tfrac13,\tfrac13)\).

Generalization: the directed \(n\)-cycle \(p_{i,i+1}=1\) (indices mod \(n\)) is connected. Starting at state \(1\) it has \(\mathbf{p}(t)=\mathbf{e}_{1+(t \bmod n)}\), which has period \(n\), while \(\pi\) is uniform. Another family is the random walk on any connected bipartite graph with sides \(L,R\), for example an even cycle \(C_{2m}\) or a grid. Starting in \(L\), \(\mathbf{p}(t)\) is supported on \(L\) for even \(t\) and on \(R\) for odd \(t\). So the mass on \(L\) alternates between \(1\) and \(0\), and \(\mathbf{p}(t)\) has no limit.

No. Let \(d_1,d_2,\dots\) be the digit sequence: blocks \(B_0,B_1,B_2,\dots\) of lengths \(2^0,2^1,2^2,\dots\), where block \(B_k\) is all \(1\)’s for even \(k\) and all \(0\)’s for odd \(k\). Let \(a(t)=\frac1t\sum_{s\le t}d_s\). Then

\begin{equation*} |a(t)-a(t+1)|=\frac{|a(t)-d_{t+1}|}{t+1}\le\frac{1}{t+1}\to0 . \end{equation*}

Block \(B_k\) ends at \(t_k=2^{k+1}-1\). The number of \(1\)’s up to \(t_k\) is \(\sum_{j\le k,\ j\text{ even}}2^j\), which gives

\begin{equation*} \begin{aligned} k\text{ even}:&\quad a(t_k)=\frac{(2^{k+2}-1)/3}{2^{k+1}-1}\ \longrightarrow\ \frac23,\\ k\text{ odd}:&\quad a(t_k)=\frac{(2^{k+1}-1)/3}{2^{k+1}-1}=\frac13 . \end{aligned} \end{equation*}

Hence \(\limsup a(t)\ge\frac23>\frac13\ge\liminf a(t)\), and \(a(t)\) does not converge.

(a) \(\pi=\frac{1}{211}(16,24,36,54,81)\). The chain is connected, so by Lemma 4.3 it is enough to solve detailed balance \(\pi_i\,p_{i,i+1}=\pi_{i+1}\,p_{i+1,i}\):

\begin{equation*} \pi_{i+1}=\frac{0.6}{0.4}\,\pi_i=\frac32\,\pi_i \ \Rightarrow\ \pi\propto\left(1,\tfrac32,\tfrac94,\tfrac{27}{8},\tfrac{81}{16}\right)\propto(16,24,36,54,81), \end{equation*}

and \(16+24+36+54+81=211\).

(b) \(\pi=(0,0,0,0,1)\). This chain is not strongly connected, so Theorem 4.2 does not apply. Solving \(\pi P=\pi\) directly instead, column by column,

\begin{equation*} \begin{aligned} \pi_1&=0.5\,\pi_1 &&\Rightarrow\ \pi_1=0,\\ \pi_{i}&=0.5\,\pi_{i-1}+0.5\,\pi_i &&\Rightarrow\ \pi_i=\pi_{i-1}=0\quad(i=2,3,4),\\ \pi_5&=0.5\,\pi_4+\pi_5, && \end{aligned} \end{equation*}

so \(\pi_5=1\) is forced and the stationary distribution is unique.

The uniform distribution \(\pi_i=\frac1n\), where \(n\) is the number of states. For all \(i,j\),

\begin{equation*} \pi_i\,p_{ij}=\frac{p_{ij}}{n}=\frac{p_{ji}}{n}=\pi_j\,p_{ji}, \qquad \sum_i\pi_i=1 . \end{equation*}

The chain is connected, and symmetry makes it strongly connected. So Lemma 4.3 shows that \(\pi\) is the stationary distribution, and it is unique by Theorem 4.2.

Let \(S=\{i:p_i>q_i\}\). Since \(\sum_i p_i=\sum_i q_i=1\), we have \(\sum_{i\in S}(p_i-q_i)=\sum_{i\notin S}(q_i-p_i)\). Therefore

\begin{equation*} \begin{aligned} |\mathbf{p}-\mathbf{q}|_1 &=\sum_{i\in S}(p_i-q_i)+\sum_{i\notin S}(q_i-p_i) =2\sum_{i\in S}(p_i-q_i)\\ &=2\sum_i(p_i-q_i)^+ , \end{aligned} \end{equation*}

because \((p_i-q_i)^+\) equals \(p_i-q_i\) on \(S\) and \(0\) off \(S\).

Use MCMC (Section 4.2): run a chain with stationary distribution \(p\) and report the fraction of visited states with \(x_1=1\). The exact sum has \(2^{99}\approx6.3\times10^{29}\) terms, so it is out of reach.

The chain is Gibbs sampling (Section 4.2.2) on \(\{0,1\}^{100}\). At each step, pick a coordinate \(i\) and resample \(x_i\) from its conditional distribution given the other coordinates:

\begin{equation*} \Pr(x_i=1\mid \mathbf{x}_{-i})=\frac{p(\mathbf{x}^{i\leftarrow1})}{p(\mathbf{x}^{i\leftarrow0})+p(\mathbf{x}^{i\leftarrow1})} . \end{equation*}

Only a ratio of two values of \(p\) is needed, so \(p\) may be known only up to its normalizing constant.

Take \(f(\mathbf{x})=\mathbf{1}[x_1=1]\), so that \(E(f)=p(x_1=1)\) and \(f_{\max}=1\). Let \(\gamma\) be the average of \(f\) over the states of a \(t\)-step walk. By (4.2),

\begin{equation*} \left|\,p(x_1=1)-E(\gamma)\,\right|\le\|p-\mathbf{a}(t)\|_1 , \end{equation*}

and this tends to \(0\) by Theorem 4.2. The bias is small once \(t\) exceeds the \(\varepsilon\)-mixing time (Definition 4.1), which Theorem 4.5 bounds through the conductance. Then \(\hat p(x_1=1)=\gamma\) and \(\hat p(x_1=0)=1-\gamma\). To control the variance as well, run \(N\) independent walks of length beyond the mixing time and average their final values of \(f\). By a Chernoff bound, \(N=O(\varepsilon^{-2}\log(1/\delta))\) walks give error at most \(\varepsilon\) with probability at least \(1-\delta\).

The maximum degree of the lattice is \(r=4\) (at \(11\)). The Metropolis-Hastings rule of Section 4.2.1, \(p_{ij}=\frac14\min(1,p_j/p_i)\) for lattice neighbours \(i,j\), gives the following chain:

For example, \(p_{00,01}=p_{00,10}=\frac14\) and \(p_{01,00}=p_{01,02}=\frac18\). Every row sums to \(1\). Detailed balance holds on the two kinds of edges:

\begin{equation*} \begin{aligned} \text{corner–edge}:&\quad \tfrac1{16}\cdot\tfrac14=\tfrac1{64}=\tfrac18\cdot\tfrac18,\\ \text{edge–centre}:&\quad \tfrac18\cdot\tfrac14=\tfrac1{32}=\tfrac14\cdot\tfrac18 . \end{aligned} \end{equation*}

The lattice is connected and the probabilities sum to \(1\), so by Lemma 4.3 the table is the stationary distribution of this chain.

With the random-coordinate scheme of Section 4.2.2, \(p_{(x,y),(x^{\prime},y)} = \tfrac12\, p(x^{\prime} \mid y)\) and \(p_{(x,y),(x,y^{\prime})} = \tfrac12\, p(y^{\prime} \mid x)\), the self-loop collecting the two terms with \(x^{\prime}=x\), \(y^{\prime}=y\). Rows and columns \(1,4\) are outer lines \((\tfrac1{16},\tfrac1{32},\tfrac1{32},\tfrac1{16})\) with mass \(\tfrac3{16}\); rows and columns \(2,3\) are inner lines \((\tfrac1{32},\tfrac18,\tfrac18,\tfrac1{32})\) with mass \(\tfrac5{16}\). Hence the conditionals are

\begin{equation*} \text{outer: } \bigl(\tfrac13,\tfrac16,\tfrac16,\tfrac13\bigr), \qquad \text{inner: } \bigl(\tfrac1{10},\tfrac25,\tfrac25,\tfrac1{10}\bigr), \end{equation*}

and the three vertex types (corner \(\tfrac1{16}\), edge \(\tfrac1{32}\), center \(\tfrac18\)) get the following transitions; all others follow by the symmetries of the square.

Each row of the table sums to \(1\) (Check!). Edges joining two vertices of the same type carry equal probability both ways, and the two mixed types satisfy detailed balance:

\begin{equation*} \begin{aligned} \text{corner–edge: } & \tfrac1{16}\cdot\tfrac1{12} = \tfrac1{192} = \tfrac1{32}\cdot\tfrac16, \\ \text{edge–center: } & \tfrac1{32}\cdot\tfrac15 = \tfrac1{160} = \tfrac18\cdot\tfrac1{20}. \end{aligned} \end{equation*}

Every transition along a row or column is positive, so the chain is connected, and by Lemma 4.3 its stationary distribution is the table.

Reduce the integral to a volume times a telescoping product of expectations, each estimated by an MCMC walk on a grid of the region. Let \(f \ge 0\) be the polynomial on the convex body \(R \subset \mathbb{R}^d\); if \(f\) vanishes somewhere on \(R\), integrate \(f+\eta\) for a small \(\eta>0\) and subtract \(\eta\,\mathrm{Vol}( R)\), so assume \(0 < f_{\min} \le f \le f_{\max}\) on \(R\). Put \(Z(\beta) = \int_R f^{\beta}\,dx\) and pick \(0=\beta_0<\beta_1<\cdots<\beta_k=1\) with steps \(\beta_i-\beta_{i-1} \le \delta = 1/\ln(f_{\max}/f_{\min})\) and \(k = \lceil \ln(f_{\max}/f_{\min})\rceil\). Then

\begin{equation*} \begin{aligned} \int_R f\,dx = Z(1) &= Z(0)\prod_{i=1}^{k}\frac{Z(\beta_i)}{Z(\beta_{i-1})} = \mathrm{Vol}( R)\prod_{i=1}^{k} \mathbb{E}_{\pi_{i-1}}\bigl[f^{\beta_i-\beta_{i-1}}\bigr],\\ \pi_{i-1}(x) &\propto f(x)^{\beta_{i-1}} \text{ on } R . \end{aligned} \end{equation*}

(i) \(\mathrm{Vol}( R)\) is estimated exactly as in Section 4.3, with concentric spheres and uniform samples from \(S_i \cap R\).

(ii) For each factor, impose a grid on \(R\) and run Metropolis-Hastings (Section 4.2.1) on the grid graph with target \(\pi_{i-1}\), i.e. accept a move \(x \to y\) with probability \(\min\bigl(1, (f(y)/f(x))^{\beta_{i-1}}\bigr)\). Average \(g = f^{\beta_i-\beta_{i-1}}\) over the states visited. By (4.2) the bias is at most \(g_{\max}\,\|\pi_{i-1}-a(t)\|_1\), which Theorem 4.5 makes small once \(t\) exceeds the mixing time. Since \(g_{\max}/g_{\min} \le (f_{\max}/f_{\min})^{\delta} = e\), the estimated quantity never exceeds \(e\) times its mean. So \(O(1/{\varepsilon^{\prime}}^2)\) samples give each factor to within \(1 \pm \varepsilon^{\prime}\), and taking \(\varepsilon^{\prime} = \varepsilon/k\) gives the whole product to within \(1 \pm O(\varepsilon)\), just as the sphere ratios bounded by \(e\) do in Section 4.3.

The total cost is polynomial in \(d\), \(\ln(f_{\max}/f_{\min})\) and \(1/\varepsilon\), provided each walk mixes rapidly. For the uniform target this is the convex-body result cited in Section 4.3. For \(\pi \propto f^{\beta}\) it holds when \(f^{\beta}\) is log-concave on \(R\) (for example, a product of nonnegative linear forms); for a general polynomial, rapid mixing has to be checked through the normalized conductance \(\Phi\) of Theorem 4.5.

The new chain is \(P^{\prime} = \tfrac12(I+P)\). It has the same stationary distribution \(\pi\), and it converges at most about twice as slowly, with the gain that \(p(t)\) itself now converges. For \(x \ne y\),

\begin{equation*} \pi_x p^{\prime}_{xy} = \tfrac12\,\pi_x p_{xy} = \tfrac12\,\pi_y p_{yx} = \pi_y p^{\prime}_{yx}, \end{equation*}

and \(P^{\prime}\) is connected because \(P\) is, so by Lemma 4.3 \(\pi\) is the stationary distribution of \(P^{\prime}\). (Directly: \(\pi P^{\prime} = \tfrac12(\pi + \pi P) = \pi\).)

Convergence time. Self-loops do not enter the flow in Definition 4.2, so \(\Phi^{\prime}(S) = \tfrac12\Phi(S)\) for every \(S\). Hence \(\Phi^{\prime} = \Phi/2\), and the bound of Theorem 4.5 grows only by the constant factor \(4\). More precisely, reversibility makes \(P\) similar to the symmetric matrix \(D^{1/2}PD^{-1/2}\) with \(D = \mathrm{diag}(\pi)\), so \(P\) has real eigenvalues \(1 = \lambda_1 > \lambda_2 \ge \cdots \ge \lambda_n \ge -1\). Then \(P^{\prime}\) has eigenvalues

\begin{equation*} \lambda_i^{\prime} = \tfrac12(1+\lambda_i) \in [0,1], \qquad 1-\lambda_2^{\prime} = \tfrac12(1-\lambda_2), \end{equation*}

and \(\|p(t)-\pi\|\) decays like \(\max_{i\ge 2}|\lambda^{\prime}_i|^{t} = \bigl(\tfrac{1+\lambda_2}{2}\bigr)^{t}\). The spectral gap halves, so the convergence time at most doubles: after \(t\) lazy steps the walk has made a \(\mathrm{Bin}(t,\tfrac12)\) number of real moves. In exchange there are no negative eigenvalues. When \(\lambda_n\) is at or near \(-1\) (a bipartite or nearly periodic chain such as the two-state chain of Exercise 4.1), \(p(t)\) under \(P\) oscillates and does not converge, while under \(P^{\prime}\) it converges at rate \(\bigl(\tfrac{1+\lambda_2}{2}\bigr)^{t}\).

The stationary vector is \(\mathbf{p}\), but the walk is not efficient: on the hypercube it mixes in time at least \(2^n\). The chain is well defined because \(p_{ii} = 1 - \sum_{j \sim i} p_j \ge 1 - (1 - p_i) \ge 0\). For adjacent \(i \ne j\),

\begin{equation*} p_i\,p_{ij} = p_i\,p_j = p_j\,p_{ji}, \end{equation*}

and the chain is connected (take \(p_j > 0\) for all \(j\)), so by Lemma 4.3 \(\mathbf{p}\) is the stationary distribution.

On the hypercube with \(p_j = 2^{-n}\), the walk leaves its current vertex with probability \(n2^{-n}\) per step, so it waits \(2^n/n\) steps on average before its first move. Quantitatively, let \(S = \{x : x_1 = 0\}\). Then \(\pi(S) = \tfrac12\), and \(2^{n-1}\) edges cross the cut, each with \(\pi_x p_{xy} = 2^{-n}\cdot 2^{-n}\), so

\begin{equation*} \Phi \le \Phi(S) = \frac{2^{n-1}\cdot 2^{-2n}}{1/2} = 2^{-n}. \end{equation*}

By the lower bound of Section 4.4 the mixing time is at least \(1/\Phi(S) = 2^n\), exponential in \(n\).

Without self-loops the forced chain is \(p_{12} = \tfrac23\), \(p_{13} = \tfrac13\), \(p_{21} = p_{31} = 1\), \(p_{23} = p_{32} = 0\), so the edge \(\{2,3\}\) must get probability \(0\). With the self-loop \(p_{11} = \tfrac12\) all three edges can be used: \(p_{12} = \tfrac5{12}\), \(p_{13} = \tfrac1{12}\), \(p_{21} = \tfrac58\), \(p_{23} = \tfrac38\), \(p_{31} = \tfrac14\), \(p_{32} = \tfrac34\).

Work with edge weights as in Section 4.4, normalized so that \(\pi_x = w_x = \sum_y w_{xy}\). The flow of each edge is \(w_{xy} = \pi_x p_{xy}\), and \(p_{xy} = w_{xy}/\pi_x\).

(i) No self-loops. Write \(a = w_{12}\), \(b = w_{13}\), \(c = w_{23}\). Then

\begin{equation*} a+b = \tfrac12,\quad a+c = \tfrac13,\quad b+c = \tfrac16 \;\Longrightarrow\; a = \tfrac13,\ b = \tfrac16,\ c = 0 . \end{equation*}

Dividing by \(\pi_x\) gives the chain above. The edge \(\{2,3\}\) cannot carry positive probability, even for a non-reversible chain. With no self-loops, \(\pi_1 = \pi_2 p_{21} + \pi_3 p_{31} \le \pi_2 + \pi_3 = \tfrac12\), and \(\pi_1 = \tfrac12\) forces \(p_{21} = p_{31} = 1\). So as printed the requirement that every pair be joined is unattainable, and the self-loop of the second part removes the obstruction.

(ii) Self-loop with \(p_{11} = \tfrac12\), i.e. \(w_{11} = \pi_1 p_{11} = \tfrac14\). Then

\begin{equation*} \tfrac14 + a + b = \tfrac12,\quad a+c = \tfrac13,\quad b+c = \tfrac16 \;\Longrightarrow\; a = \tfrac5{24},\ b = \tfrac1{24},\ c = \tfrac18 . \end{equation*}

Dividing by \(\pi_x\) gives the chain stated above, whose rows sum to \(1\). All three edges are positive, so the chain is connected, and \(\pi_x p_{xy} = w_{xy} = \pi_y p_{yx}\) shows by Lemma 4.3 that \(\left(\tfrac12,\tfrac13,\tfrac16\right)\) is stationary. Check:

\begin{equation*} \begin{aligned} \pi_1 &= \tfrac12\cdot\tfrac12 + \tfrac13\cdot\tfrac58 + \tfrac16\cdot\tfrac14 = \tfrac12,\\ \pi_2 &= \tfrac12\cdot\tfrac5{12} + \tfrac16\cdot\tfrac34 = \tfrac13 . \end{aligned} \end{equation*}

\(p_{ij} = \tfrac12\min\left(1, \tfrac{p_j}{p_i}\right)\): by Section 4.2.1 with maximum degree \(r = 2\), each neighbour is proposed with probability \(\tfrac12\) and accepted with probability \(\min(1, p_j/p_i)\).

The self-loops are \(p_{11} = \tfrac12\), \(p_{22} = \tfrac14\), \(p_{33} = 0\). The equal flows in the two directions of each edge are the detailed balance that, by Lemma 4.3, makes \(\left(\tfrac12,\tfrac13,\tfrac16\right)\) stationary. In stationarity a step moves along \(\{1,2\}\) with probability \(\tfrac13\), along \(\{1,3\}\) with probability \(\tfrac16\) and along \(\{2,3\}\) with probability \(\tfrac16\), and it stays put with probability \(\tfrac12\cdot\tfrac12 + \tfrac13\cdot\tfrac14 = \tfrac13\).

The square (4-cycle) \(00 - 01 - 11 - 10 - 00\) is bad, and the path \(01 - 00 - 11 - 10\) is good. Gibbs sampling is stuck in the same way as the bad graph, because its graph is that square.

Bad: on the square every neighbour of \(00\) and of \(11\) has probability \(0\). Every proposal from \(00\) is accepted with probability \(\min(1, 0/\tfrac12) = 0\), so

\begin{equation*} p_{00,00} = p_{11,11} = 1 . \end{equation*}

A walk started in the support never leaves its starting state, and both \(\delta_{00}\) and \(\delta_{11}\) are stationary. The set \(S = \{00\}\) has \(\Phi(S) = 0\), so the mixing time, bounded below by \(1/\Phi(S)\) (Section 4.4), is infinite. Averages such as \(\mathbb{E}(x_1) = \tfrac12\) are estimated as \(0\) or \(1\).

Good: on the path (\(r = 2\)), \(00\) proposes \(11\) with probability \(\tfrac12\) and accepts it, and it rejects \(01\), so

\begin{equation*} p_{00,11} = p_{11,00} = \tfrac12, \qquad p_{00,00} = p_{11,11} = \tfrac12 . \end{equation*}

From either support state the distribution after a single step is exactly \(\left(\tfrac12,\tfrac12\right)\) on \(\{00,11\}\), so the chain mixes in one step. The two zero-probability leaves are left at the first accepted move (ratio \(p_j/p_i\) read as \(1\) when \(p_i = 0\)) and never revisited.

Gibbs: Section 4.2.2 changes one coordinate at a time, so its graph is the square. The conditionals are point masses, \(p(x_1 \mid x_2) = \mathbf{1}[x_1 = x_2]\), and likewise for \(x_2\), so resampling either coordinate of \(00\) returns \(00\), and similarly for \(11\). Gibbs never moves between the two modes, because that requires changing both perfectly correlated coordinates at once.

Gibbs sampling has the right stationary distribution but mixes exponentially slowly: it leaves \(\mathbf{0}\) with probability exactly \(2^{-100}\) per step, so a run stays in whichever of the two halves \(\{\mathbf{0}\}\), \(\{0,1\}^{100}\setminus\{\mathbf{0}\}\) it starts in for about \(2^{100}\) steps.

Write \(N=2^{100}\) and \(q=\frac{1}{2(N-1)}\). With a random coordinate \(i\) updated from its conditional (Section 4.2.2), the three kinds of states behave as follows.

(i) At \(\mathbf{x}=\mathbf{0}\): the new \(x_i\) is \(1\) with probability

\begin{equation*} \frac{q}{\tfrac12+q}=\frac{1/(2(N-1))}{N/(2(N-1))}=\frac1N=2^{-100}, \end{equation*}

whatever \(i\) is, so \(\Pr(\text{leave }\mathbf{0})=2^{-100}\) and the expected stay is \(2^{100}\) steps.

(ii) At \(\mathbf{x}\) of weight \(\ge 2\): both values of \(x_i\) give nonzero states of probability \(q\), so \(x_i\) is set to a fair coin, exactly as in the lazy random walk on the hypercube.

(iii) At \(\mathbf{x}=\mathbf{e}_j\): updating \(i=j\) moves to \(\mathbf{0}\) with probability \(\frac{1/2}{1/2+q}=1-2^{-100}\). Any other \(i\) is a fair coin as in (ii).

By Lemma 4.3 the chain is reversible with stationary distribution \(p\), so it is correct in the limit. Take \(S=\{\mathbf{0}\}\), so \(\pi(S)=\pi(\bar S)=\frac12\). Then Definition 4.2 gives

\begin{equation*} \Phi(S)=\frac{p(\mathbf{0})\cdot 2^{-100}}{1/2}=2^{-100}, \end{equation*}

and by the argument after Definition 4.2 the mixing time is at least \(1/\Phi(S)=2^{100}\). In practice:

(i) A run started at \(\mathbf{0}\) stays at \(\mathbf{0}\) for about \(2^{100}\) steps, which is a sample from the point mass and not from \(p\).

(ii) A run started elsewhere samples uniformly from the nonzero states. Its rate of entering \(\mathbf{0}\) is \(2^{-101}/\tfrac12=2^{-100}\) per step, so it never sees the state that carries half the mass.

Estimates such as \(\Pr(\mathbf{x}=\mathbf{0})\) come out as \(0\) or \(1\) instead of \(\frac12\) unless the run is far longer than \(2^{100}\) steps.

Run Metropolis-Hastings (Section 4.2.1) on edge sets \(F\subseteq E(G)\), with target \(\pi(F)\propto|F|\) and a symmetric toggle/swap proposal.

For \(k\ge2\) a connected subgraph is determined by its edge set \(F\), because its vertices are the endpoints of \(F\). The state space is therefore

\begin{equation*} \Omega=\{F\subseteq E(G): (V(F),F)\text{ is connected and } |V(F)|=k\}, \end{equation*}

and every state has \(|F|\ge k-1\ge1\). (For \(k=1\) every subgraph has \(0\) edges, so the request is empty.) Let \(m=|E(G)|\). Start from any \(F_0\in\Omega\), for example the edges of a BFS tree on the first \(k\) vertices reached from some vertex, and repeat:

pick e uniformly from E(G)
with prob 1/2:  F' = F xor {e}                 (toggle e)
otherwise:      pick f uniformly from E(G);
                if e in F and f not in F:  F' = F - {e} + {f}   (swap)
                else F' = F
if F' is in Omega (k vertices, connected):
    accept F' with probability min(1, |F'|/|F|)
output F after many steps

The proposal is symmetric. A toggle and its reverse both have probability \(\frac1{2m}\), and the swap \(F\to F-e+f\) and its reverse \(F^{\prime}\to F^{\prime}-f+e\) both have probability \(\frac1{2m^2}\). Hence for neighbouring states

\begin{equation*} \begin{aligned} \pi(F)\,p_{FF^{\prime}}&=c\,|F|\cdot q\cdot\min\!\Big(1,\frac{|F^{\prime}|}{|F|}\Big)\\ &=c\,q\,\min(|F|,|F^{\prime}|)=\pi(F^{\prime})\,p_{F^{\prime}F}, \end{aligned} \end{equation*}

so by Lemma 4.3 the stationary distribution is \(\pi(F)\propto|F|\), provided the move graph on \(\Omega\) is connected.

Why the move graph should be connected:

(i) Toggles inside a fixed vertex set \(U\) connect all connected spanning edge sets of \(G[U]\), since one can add edges up to \(E(G[U])\) and then delete down to any target.

(ii) Swaps move the vertex set: remove the edge of a leaf of a spanning tree and attach a new neighbour by one edge.

We do not prove that (ii) connects all connected \(k\)-sets of vertices (the exercise waives proofs). A valid proposal is accepted with probability at least \(\frac{k-1}{k}\), since a deletion goes from \(|F|\ge k\) to \(|F|-1\) edges.

Treat a multigraph as a multiset of \(1500\) edges on the labelled vertex set \([1000]\), where a loop adds \(2\) to the degree. The key fact is that whenever \(G\neq H\), some flip of \(G\) increases the multiset intersection \(|G\cap H|\).

(1) Let \(G\neq H\). Pick an edge \(\{u,v\}\) whose multiplicity in \(H\) exceeds its multiplicity in \(G\). Degrees agree, and \(\deg_G=\deg_{G\cap H}+\deg_{G\setminus H}\) (likewise for \(H\)), so for every vertex \(w\)

\begin{equation*} \deg_{G\setminus H}(w)=\deg_{H\setminus G}(w). \end{equation*}

Two cases:

(i) \(u\neq v\). Both \(u\) and \(v\) have positive degree in \(G\setminus H\), so there are edges \(\{u,x\},\{v,y\}\in G\setminus H\). They are distinct, because the only edge containing both \(u\) and \(v\) is \(\{u,v\}\), which is not in \(G\setminus H\). Flip \((u,x),(v,y)\to(u,v),(x,y)\).

(ii) \(u=v\). Then \(\deg_{G\setminus H}(u)\ge2\), and \(G\setminus H\) has no loop at \(u\) because \(m_G(u,u)<m_H(u,u)\). So there are two distinct edges \(\{u,x\},\{u,y\}\in G\setminus H\). Flip \((u,x),(u,y)\to(u,u),(x,y)\).

Either way two edges of \(G\setminus H\) are removed, which costs nothing in \(|G\cap H|\), and a needed copy of \(\{u,v\}\) is gained. So \(|G\cap H|\) rises by at least \(1\). Induction reaches \(H\), and every flip can be undone by a flip, so the network is connected.

(2) Moves must be counted with care. The naive walk (“pick two edges and a re-pairing uniformly”) is not uniform over multigraphs, because edge multiplicities make the number of choices from \(G\) to \(G^{\prime}\) differ from the number from \(G^{\prime}\) to \(G\). On \(2\) vertices, the triple edge \(A\) and \(B=\{11,22,12\}\) give \(P(A\to B)=\frac12\) but \(P(B\to A)=\frac13\), so \(\pi_A/\pi_B=\frac23\). Instead let \(D=1500\cdot1499\), which bounds the number of (unordered pair, re-pairing) choices and hence every network degree. Use

\begin{equation*} \begin{aligned} p_{GG^{\prime}}&=\frac1D \quad (G^{\prime}\text{ a network neighbour of }G),\\ p_{GG}&=1-\frac{\deg(G)}{D}. \end{aligned} \end{equation*}

This is implementable by proposing a uniform (pair, re-pairing) choice and accepting the resulting \(G^{\prime}\neq G\) with probability \(1/N(G,G^{\prime})\), where \(N(G,G^{\prime})\) is the number of choices that produce \(G^{\prime}\). Since \(p_{GG^{\prime}}=p_{G^{\prime}G}\), the uniform vector satisfies \(\pi_Gp_{GG^{\prime}}=\pi_{G^{\prime}}p_{G^{\prime}G}\). The network is connected by part 1, so Lemma 4.3 (with Theorem 4.2) makes the uniform distribution the unique stationary one.

(3) Each flip in part 1 raises \(|G\cap H|\) by at least \(1\), starting from at least \(0\) and ending at \(1500\). Hence

\begin{equation*} \operatorname{diam}\le 1500=\tfrac32\cdot1000 . \end{equation*}

(4) Use rejection: run the chain of part 2, and if the output is disconnected discard it and draw again. The uniform distribution on all cubic multigraphs, conditioned on connectedness, is the uniform distribution on connected ones. The expected number of draws is \(1/\Pr(\text{connected})\).

The walk estimates \(1.0032\) times the true number of grid points in the ball, which gives \(\mathrm{Vol}\approx 0.02720\). The exact value is \(V_{20}=\pi^{10}/10!=0.025807\), and the \(5\%\) gap is almost entirely the bias of a \(0.1\) grid.

The algorithm is the multi-phase method of Section 4.3 with boxes in place of spheres (footnote 15). In grid units the ball is \(B=\{k\in\mathbb{Z}^{20}:|k|^2\le100\}\), and \(\mathrm{Vol}\approx|B|\cdot 0.1^{20}\). Let

\begin{equation*} A(b)=\{k\in B:|k_i|\le b_i\ \forall i\}. \end{equation*}

Start from \(b=(2,\dots,2)\). Here \(|A(b)|=5^{20}\) exactly, since \(20\cdot2^2\le100\). Then raise one coordinate of \(b\) by \(1\) at a time, \(160\) times, until \(b=(10,\dots,10)\) and \(A(b)=B\). In each phase run the walk of Section 4.3 on \(A(b_{\rm new})\): pick one of the \(40\) grid neighbours, move there if it lies in the set, otherwise stay. Its transition matrix is symmetric, so its stationary distribution is uniform by Lemma 4.3. The phase records the fraction of samples lying in \(A(b_{\rm old})\). Ball slices shrink as \(|k_i|\) grows, so each ratio is at least

\begin{equation*} \frac{|A(b_{\rm old})|}{|A(b_{\rm new})|}\ge\frac{2b_i+1}{2b_i+3}\ge\frac57, \end{equation*}

and rejection sampling stays efficient. The estimate is \(|B|\approx5^{20}\big/\prod(\text{fractions})\). In the program, \(2000\) walkers are warm-started from the previous phase, run \(1500\) steps per phase, and are sampled every \(25\) steps in the second half:

import numpy as np, math
rng = np.random.default_rng(0)
d, R2, W, STEPS = 20, 100, 2000, 1500        # grid units: spacing 0.1, radius 10
b = np.full(d, 2)                            # start box |k_i| <= 2 (inside ball)
x = np.zeros((W, d), dtype=int); log_ratio = 0.0
for i in np.repeat(np.arange(d), 8):         # enlarge coordinate i by one unit
    old = b.copy(); b[i] += 1; hits = tot = 0
    for t in range(STEPS):
        j = rng.integers(d, size=W); s = rng.choice([-1, 1], size=W)
        y = x.copy(); y[np.arange(W), j] += s
        ok = (np.abs(y) <= b).all(1) & ((y**2).sum(1) <= R2)
        x[ok] = y[ok]
        if t >= STEPS // 2 and t % 25 == 0:
            hits += (np.abs(x) <= old).all(1).sum(); tot += W
    log_ratio += math.log(hits / tot)        # |A_old| / |A_new|
count = 5.0**d / math.exp(log_ratio)          # grid points in the ball
print(count, count * 0.1**d, math.pi**10 / math.factorial(10))

With seed \(0\) the \(4.8\times10^{8}\) proposals gave

\begin{equation*} \begin{aligned} |B|_{\rm est}&=2.7202\times10^{18}, &\qquad |B|_{\rm exact}&=2.7114\times10^{18},\\ \mathrm{Vol}_{\rm est}&=|B|_{\rm est}\,10^{-20}=0.027202, & V_{20}&=0.025807, \end{aligned} \end{equation*}

where the exact count comes from expanding \(\big(\sum_{|j|\le10}q^{j^2}\big)^{20}\) up to \(q^{100}\).

(1) \(\Theta(n^2)\) for fixed \(\varepsilon\): the bridge is crossed only once every \(\approx n^2\) steps. (2) \(O(1/\varepsilon)\), independent of \(n\). Take both cliques in (1) to have \(n\) vertices, since the exercise does not fix the sizes. Every edge carries flow \(\pi_xp_{xy}=\frac1{2m}\) (Section 4.4.1, connected undirected graph).

(1) Here \(m=n(n-1)+1\). Let \(S\) be one clique. By symmetry \(\pi(S)=\frac12\), and only the bridge leaves \(S\), so

\begin{equation*} \Phi(S)=\frac{1/(2m)}{1/2}=\frac1m . \end{equation*}

Section 4.4.1 gives \(\Phi\ge\frac1m\) for every connected graph, so \(\Phi=\frac1m\approx n^{-2}\). The bound after Definition 4.2 gives mixing time at least \(1/\Phi=m=\Omega(n^2)\). Theorem 4.5, with \(\pi_{\min}=\frac{n-1}{2m}\), gives the upper bound \(O(m^2\ln n/\varepsilon^3)=O(n^4\ln n/\varepsilon^3)\). The truth is the lower bound. Computing \(\mathbf{a}(t)\) exactly from the worst start (by linearity the worst start is a vertex) gives these values of the \(\frac14\)-mixing time:

So \(t_{\rm mix}(1/4)\approx2n^2\). (These are numerical values. We have no proof of a matching \(O(n^2)\) upper bound with this chapter’s tools.)

(2) Let the clique be \(K\), let \(u\in K\), and let \(v\) be the pendant vertex, so \(m=\binom n2+1\). Take \(S\) with \(\pi(S)\le\frac12\) and put \(s=|S\cap K|\).

(i) \(s=0\). Then \(S=\{v\}\) and \(\Phi(S)=1\).

(ii) \(s\ge1\). From \(s(n-1)\le\mathrm{vol}(S)\le m\) we get \(s\le n/2\) (for \(n\ge4\)). At least \(s(n-s)\) edges leave \(S\), and \(\mathrm{vol}(S)\le s(n-1)+2\le s(n+1)\), so

\begin{equation*} \Phi(S)\ge\frac{s(n-s)}{s(n+1)}\ge\frac{n/2}{n+1}\ge\frac13 . \end{equation*}

So \(\Phi\ge\frac13\), and Theorem 4.5 with \(\pi_{\min}=\frac1{2m}\ge n^{-2}\) gives \(O(\ln n/\varepsilon^3)\). From \(v\) the walk is at \(u\) after one step and is essentially uniform on \(K\) after two, so \(\|\mathbf{a}(t)-\pi\|_1\approx4/t\). From the start \(v\) we also have \(a_v(t)\ge1/t\), which forces \(t\ge1/(\varepsilon+\frac1{2m})\). Exact computation gives \(t_{\rm mix}(1/4)=15,16,16,16,16\) and \(t_{\rm mix}(1/10)=37,39,40,40,40\) for \(n=10,20,40,80,160\), so the mixing time is \(\Theta(1/\varepsilon)\), about \(4/\varepsilon\) (the lower bound is proved above; the matching upper bound is numerical).

Both graphs have \(\Phi=\Omega(1)\) with high probability, so by Theorem 4.5 both mix in \(O(\ln n/\varepsilon^3)\) steps; both also need \(\Omega(\ln n/\ln\ln n)\). As printed, (1) sits exactly at the isolated-vertex threshold of Theorem 8.6, where the graph is disconnected (infinite mixing time) with probability \(\to1-e^{-1}\), so we treat the intended \(p=c\ln n/n\) with \(c>16\).

Throughout, for \(|S|=s\le n/2\) we use \(\min(\pi(S),\pi(\bar S))\le\pi(S)\), so \(\Phi(S)\ge e(S,\bar S)/\mathrm{vol}(S)\). For \(|S|>n/2\) we use \(\Phi(S)=\Phi(\bar S)\).

(1) Let \(\mu=(n-1)p\approx c\ln n\). By Theorem 12.3 (\(\delta=1\)) and Theorem 12.4 (\(\gamma=\frac12\)),

\begin{equation*} \begin{aligned} \Pr\big(\exists x:\ d_x\notin[\tfrac\mu2,2\mu]\big) &\le n\big((e/4)^{\mu}+e^{-\mu/8}\big)\to0 . \end{aligned} \end{equation*}

For fixed \(S\), \(e(S,\bar S)\sim\mathrm{Bin}(s(n-s),p)\) has mean at least \(\frac{snp}2\), so Theorem 12.4 and a union bound give

\begin{equation*} \sum_{s\le n/2}\binom ns e^{-snp/16}\le\sum_{s\ge1}n^{s(1-c/16)}\to0 . \end{equation*}

Hence with high probability every \(S\) with \(s\le n/2\) satisfies \(e(S,\bar S)\ge\frac{snp}{4}\) and \(\mathrm{vol}(S)\le2s\mu\), so

\begin{equation*} \Phi(S)\ge\frac{snp/4}{2s(n-1)p}\ge\frac18 . \end{equation*}

Also \(\pi_{\min}\ge\frac{\mu/2}{2n\mu}=\frac1{4n}\), so Theorem 4.5 gives \(O(\ln n/\varepsilon^3)\). In the other direction, \(\mathbf{a}(t)\) started at \(x\) lives on the ball \(B(x,t)\). That ball has \(|B|\le(2\mu+1)^t\) vertices, each of stationary probability at most \(\frac4n\), so \(\|\mathbf{a}(t)-\pi\|_1\ge2(1-4(2\mu+1)^t/n)\). Being \(\varepsilon\)-close therefore forces \(t=\Omega(\ln n/\ln\ln n)\).

(2) Let \(f(x)\) be the random endpoint chosen by \(x\), so \(m=2n\) and \(d_x=3+\mathrm{in}(x)\). Let \(S\) consist of \(k\) arcs of the circle, let \(\mathrm{in}(S)=\#\{x:f(x)\in S\}\), and let \(s_{\rm in}=\#\{x\in S:f(x)\in S\}\). Then

\begin{equation*} \begin{aligned} \Phi(S)&\ge\frac{e(S,\bar S)}{\mathrm{vol}(S)}\\ &=\frac{2k+(s-s_{\rm in})+(\mathrm{in}(S)-s_{\rm in})}{3s+\mathrm{in}(S)} . \end{aligned} \end{equation*}

Fix \(\beta=\delta=0.02\).

(i) If \(\mathrm{in}(S)>2s\), the numerator is at least \(\mathrm{in}(S)-s\), so \(\Phi(S)\ge\frac15\).

(ii) If \(\mathrm{in}(S)\le2s\) and either \(k\ge\delta s\) or \(s_{\rm in}\le(1-\beta)s\), the numerator is at least \(\beta s\), so \(\Phi(S)\ge\frac{\beta s}{5s}=\frac1{250}\).

(iii) Otherwise \(S\) has fewer than \(\delta s\) arcs and more than \((1-\beta)s\) of its own choices fall in \(S\). For fixed \(S\) this has probability at most \(\binom s{\beta s}(s/(n-1))^{(1-\beta)s}\le e^{1/2}\binom s{\beta s}(s/n)^{(1-\beta)s}\) (as \(s\le n/2\)); the constant \(e^{1/2}\) is dropped below. There are at most \(s\binom n{\delta s}\binom s{\delta s}\) sets with \(|S|=s\) and fewer than \(\delta s\) arcs, so

\begin{equation*} \begin{aligned} \Pr(\text{(iii) for some }S) &\le\sum_{s\le n/2}s\Big[(e/\beta)^{\beta}(e/\delta)^{2\delta} \big(\tfrac sn\big)^{1-\beta-\delta}\Big]^s\\ &=\sum_{s\le n/2}s\Big[1.343\big(\tfrac sn\big)^{0.96}\Big]^s\to0, \end{aligned} \end{equation*}

since the bracket is at most \(0.69\) for \(s\le n/2\) and is \(O(n^{-0.96})\) at \(s=1\). Numerically the sum is \(1.9\times10^{-4}\) at \(n=10^4\) and \(2.3\times10^{-6}\) at \(n=10^6\).

So with high probability \(\Phi\ge\frac1{250}\). With \(\pi_{\min}\ge\frac3{4n}\), Theorem 4.5 gives \(O(\ln n/\varepsilon^3)\). Every in-degree is at most \(\ln n\) with high probability, because \(\Pr(\mathrm{in}(x)\ge L)\le1/L!\). The ball argument of part 1 then again gives \(\Omega(\ln n/\ln\ln n)\).

\(\Phi\ge\frac1{2n}\), with equality for even \(n\), so the \(\varepsilon\)-mixing time is \(O(n^2\ln n/\varepsilon^3)\). As printed, parts 1 and 2 are lower bounds rather than exact minima, and part 1 also needs \(|\bar S|\ge n^2/4\) (otherwise \(S=V\) is a counterexample). We prove these corrected forms.

The chain is the one of Section 4.4.1: \(p_{ij}=\frac14\) between lattice neighbours, self-loops on the boundary, and \(\pi_i=\frac1{n^2}\). Call a row or column pure if it lies entirely in \(S\) or entirely in \(\bar S\), and mixed otherwise. A mixed line contains at least one boundary edge, and these edges are distinct for distinct lines.

(1) Suppose \(|S|,|\bar S|\ge n^2/4\). If no row or no column is pure, there are \(n\) mixed lines, giving at least \(n\) edges. Otherwise a pure row and a pure column meet, so all pure lines have the same type, say \(S\). Let \(a\) rows and \(b\) columns be pure. Then \(\bar S\) lies in the other \((n-a)(n-b)\) cells, and AM-GM gives

\begin{equation*} \#\text{edges}\ge(n-a)+(n-b)\ge2\sqrt{(n-a)(n-b)}\ge2\sqrt{|\bar S|}\ge n . \end{equation*}

(With type \(\bar S\), use \(|S|\ge n^2/4\).) The first \(\lfloor n/2\rfloor\) columns attain \(n\).

(2) Let \(s=|S|\le n^2/4\), and let \(R\) and \(C\) be the numbers of rows and columns that meet \(S\), so \(s\le RC\).

(i) No line is pure-\(S\). Then all \(R+C\) such lines are mixed, so \(\#\text{edges}\ge R+C\ge2\sqrt{RC}\ge2\sqrt s\).

(ii) Some row is pure-\(S\) (columns are symmetric). Then every column meets \(S\), and at most \(s/n\le n/4\) columns are pure. So \(\#\text{edges}\ge\frac{3n}4\ge\frac32\sqrt s\), using \(\sqrt s\le n/2\).

Either way \(\#\text{edges}\ge\sqrt s\ge\lfloor\sqrt s\rfloor\).

(3) Each boundary edge carries flow \(\frac1{n^2}\cdot\frac14\). For \(|S|\le n^2/2\), Definition 4.2 gives

\begin{equation*} \Phi(S)=\frac{\#\text{edges}(S,\bar S)}{4|S|}\ge \begin{cases} \dfrac{\sqrt{|S|}}{4|S|}=\dfrac1{4\sqrt{|S|}}\ge\dfrac1{2n}, & |S|\le n^2/4,\\[2ex] \dfrac{n}{4\cdot n^2/2}=\dfrac1{2n}, & n^2/4\le|S|\le n^2/2, \end{cases} \end{equation*}

and \(\Phi(S)=\Phi(\bar S)\) when \(|S|>n^2/2\).

(4) By part 3, \(\Phi\ge\frac1{2n}\). For even \(n\), the left \(n/2\) columns have \(n\) boundary edges and \(\pi(S)=\frac12\), so

\begin{equation*} \Phi=\frac{n/(4n^2)}{1/2}=\frac1{2n}. \end{equation*}

For odd \(n\), \(\Phi=\Theta(1/n)\).

(5) Here \(\pi_{\min}=n^{-2}\), so Theorem 4.5 gives

\begin{equation*} \begin{aligned} t_{\rm mix}(\varepsilon) &=O\!\left(\frac{\ln(n^2)}{(1/(2n))^2\varepsilon^3}\right) =O\!\left(\frac{n^2\ln n}{\varepsilon^3}\right). \end{aligned} \end{equation*}

\(\Phi=\Theta\!\left(\frac{1}{dn}\right)\) and the \(\epsilon\)-mixing time is \(O\!\left(d^3n^2\ln n/\epsilon^3\right)\).

As in Section 4.4 (A lattice in d-dimensions), take vertices \([n]^d=\{0,\dots,n-1\}^d\), \(N=n^d\), each lattice edge has \(p_{xy}=\frac1{2d}\), the rest is a self-loop, and \(\pi_x=n^{-d}\). Write \(\partial S\) for the set of lattice edges with exactly one end in \(S\). Since \(\pi_xp_{xy}=\frac{n^{-d}}{2d}\) on every such edge, for \(|S|\le N/2\)

\begin{equation*} \Phi(S)=\frac{|\partial S|\,n^{-d}/(2d)}{|S|\,n^{-d}}=\frac{|\partial S|}{2d\,|S|}. \end{equation*}

(1) Let \(C_1(S)\) replace \(S\cap\ell\), on every line \(\ell\) parallel to the \(x_1\)-axis, by the initial segment \(\{x_1<k_\ell\}\) of the same size \(k_\ell=|S\cap\ell|\). Then \(|C_1(S)|=|S|\), and edges are counted line by line: (i) \(x_1\)-edges inside \(\ell\): at least \(1\) if \(0<k_\ell<n\), otherwise \(0\); after compression exactly \(1\) or \(0\). (ii) \(x_j\)-edges (\(j\ne1\)) between neighbouring parallel lines \(\ell,\ell^{\prime}\): they join points with equal \(x_1\), so their number is \(|(S\cap\ell)\,\triangle\,(S\cap\ell^{\prime})|\ge|k_\ell-k_{\ell^{\prime}}|\), with equality for two initial segments. Hence \(|\partial C_1(S)|\le|\partial S|\).

(2) Apply \(C_1,\dots,C_d\) repeatedly. A compression that changes \(S\) strictly lowers the integer \(\sum_{x\in S}\sum_i x_i\ge0\), so the process stops at a set \(S^*\) fixed by every \(C_i\). Every axis-parallel line meets \(S^*\) in an initial segment, so \(x\in S^*\) and \(x_i>0\) imply \(x-e_i\in S^*\). That is, \(S^*\) is a down-set. Its slices \(T_t=\{y:(t,y)\in S^*\}\) therefore satisfy \(T_0\supseteq T_1\supseteq\cdots\supseteq T_{n-1}\). By part 1, \(|S^*|=|S|\) and \(|\partial S^*|\le|\partial S|\).

(3) For a down-set, the \(x_1\)-edges between slices \(t\) and \(t+1\) number \(|T_t\setminus T_{t+1}|\), so they contribute exactly \(|T_0|-|T_{n-1}|\), and the in-plane edges contribute \(\sum_t|\partial_{d-1}T_t|\). With the direction \(x_1\) fixed in advance, the claim is false. Take \(d=2\), \(n=10\) and the down-set with column sizes \((k_0,\dots,k_9)=(5,4,4,\dots,4)\), so \(|S|=41\le50\). Here \(T_0=\{0,\dots,4\}\) is already optimal in one dimension, and \(|\partial S|=(5-4)+10=11\). Filling the slices up to the shape of \(T_0\) gives column sizes \((5,5,5,5,5,5,5,5,1,0)\), with \(|\partial|=5+9=14>11\) (checked by direct count). Compressing along \(x_2\) instead leaves this \(S\) unchanged (four full rows and one cell). So step 3 needs a good choice of direction, and proving that such a direction always exists is the edge-isoperimetric theorem for the grid (Bollobas-Leader). The argument does not need it: part 5 bounds every \(S\) directly.

(4) We may end up with corner boxes \([0,a)^r\times[0,n)^{d-r}\) (up to permuting coordinates) with one partial layer. With \(r=d\) this is a sub-cube of side \(a\), for which \(|\partial S|=d\,a^{d-1}\). With \(r=1\) it is a slab of width \(a\), for which \(|\partial S|=n^{d-1}\).

(5) Every \(S\) with \(|S|\le N/2\), modified or not, satisfies \(|\partial S|\ge|S|/n\). Route each ordered pair \((u,v)\) along the dimension-ordered path, which corrects \(x_1\), then \(x_2\), and so on. The path uses the edge \(\{x,x+e_i\}\) only if \(v_j=x_j\) for \(j<i\), \(u_j=x_j\) for \(j>i\), and \(u_i,v_i\) lie on opposite sides of \(x_i+\frac12\). So at most

\begin{equation*} n^{i-1}\,n^{d-i}\cdot2(x_i+1)(n-1-x_i)\le n^{d-1}\cdot\frac{n^2}{2}=\frac{n^{d+1}}{2} \end{equation*}

paths use any one edge. Each of the \(|S||\bar S|\) paths from \(S\) to \(\bar S\) uses an edge of \(\partial S\), and \(|\bar S|\ge N/2\), so

\begin{equation*} |S|\frac{N}{2}\le|S||\bar S|\le|\partial S|\frac{n^{d+1}}{2} \;\Longrightarrow\;|\partial S|\ge\frac{|S|}{n}. \end{equation*}

(6) By the formula \(\Phi(S)=|\partial S|/(2d|S|)\) above, a sub-cube of side \(a\) has \(\Phi(S)=\frac{d\,a^{d-1}}{2d\,a^d}=\frac1{2a}\), and a slab of width \(a\le n/2\) has \(\Phi(S)=\frac{n^{d-1}}{2d\,a\,n^{d-1}}=\frac1{2da}\).

(7) Part 5 gives the lower bound and the half-slab (\(a=\lfloor n/2\rfloor\)) gives the upper bound:

\begin{equation*} \frac1{2dn}\le\Phi\le\frac{1}{2d\lfloor n/2\rfloor}, \qquad\text{so}\qquad \Phi=\Theta\!\left(\frac1{dn}\right) \end{equation*}

(for even \(n\) the upper bound is \(\frac1{dn}\)).

(8) By Theorem 4.5, which applies to this random walk on an undirected graph, with \(\pi_{\min}=n^{-d}\):

\begin{equation*} t_\epsilon=O\!\left(\frac{\ln(1/\pi_{\min})}{\Phi^2\epsilon^3}\right) =O\!\left(\frac{d\ln n\cdot 4d^2n^2}{\epsilon^3}\right) =O\!\left(\frac{d^3n^2\ln n}{\epsilon^3}\right). \end{equation*}

(1) The constant functions, and only those. Let \(g_x=\sum_yp_{xy}g_y\) at every vertex, with \(p_{xy}>0\) exactly on edges, and let \(S\) be the set where \(g\) attains its maximum \(M\). If \(S\ne V\), connectedness gives an edge \((x,y)\) with \(x\in S\) and \(y\notin S\), and then

\begin{equation*} M=g_x=\sum_yp_{xy}g_y<\sum_yp_{xy}M=M, \end{equation*}

which is a contradiction. So \(S=V\), which is the max-principle argument of Section 4.5 with an empty boundary. Every constant is clearly harmonic.

(2) Here \(p_{xy}=1/d_x\) for \(y\sim x\). Stationarity gives \(q_x=\sum_{y\sim x}q_y/d_y\), so

\begin{equation*} \frac{q_x}{d_x}=\sum_{y\sim x}\frac1{d_x}\cdot\frac{q_y}{d_y}=\sum_yp_{xy}\frac{q_y}{d_y}. \end{equation*}

This is the Example in Section 4.5 with all \(c_{xy}=1\).

(3) On a connected graph, parts 1 and 2 force every stationary \(q\) to have \(q_x/d_x=c\) for a constant \(c\), so \(q=c\,d\). The harmonic functions form the one-parameter family of constants, and the normalisation \(\sum_xq_x=1\) selects the single value \(c=1/\sum_xd_x\).

(4) With \(m\) edges, \(\sum_xd_x=2m\), so

\begin{equation*} q_x=\frac{d_x}{2m}. \end{equation*}

No. It works exactly for the connected reversible chains, that is, those with a \(\pi\) satisfying \(\pi_xp_{xy}=\pi_yp_{yx}\) for all \(x,y\).

(i) Necessity. A network chain has \(p_{xy}=c_{xy}/c_x\), so

\begin{equation*} c_xp_{xy}=c_{xy}=c_{yx}=c_yp_{yx}, \end{equation*}

and \(\pi_x=c_x/c_0\) is a reversible stationary distribution (Section 4.5, Lemma 4.3).

(ii) Sufficiency. Given a reversible \(\pi\), put \(c_{xy}=\pi_xp_{xy}\), which is symmetric by detailed balance and positive exactly on the edges with \(p_{xy}>0\). Then

\begin{equation*} c_x=\sum_y\pi_xp_{xy}=\pi_x, \qquad \frac{c_{xy}}{c_x}=p_{xy}, \end{equation*}

so the network with resistances \(r_{xy}=1/(\pi_xp_{xy})\) reproduces the chain. A self-loop \(p_{xx}>0\) becomes a resistor from \(x\) to itself, which carries no current.

A counterexample is the directed 3-cycle \(p_{12}=p_{23}=p_{31}=1\). It is connected with \(\pi=(\frac13,\frac13,\frac13)\), but \(\pi_1p_{12}=\frac13\ne0=\pi_2p_{21}\), so by (i) no electrical network has it as its random walk. More generally, any chain with \(p_{xy}>0=p_{yx}\) fails.

(1) \(\frac14\). By Section 4.5 the probability is the voltage \(v_x\) with \(v_1=1\) and \(v_5=0\). On a path, harmonicity means \(v_x=\frac12(v_{x-1}+v_{x+1})\), so \(v\) is linear:

\begin{equation*} v_x=\frac{5-x}{4},\qquad v_4=\frac14. \end{equation*}

(2) \(\frac25\). The swap \(2\leftrightarrow3\), \(4\leftrightarrow6\) is a graph automorphism fixing \(1\) and \(5\). By uniqueness of the harmonic function (Section 4.5), \(v_2=v_3=\alpha\) and \(v_4=v_6=\beta\). Vertices \(2,3,4,6\) all have degree \(4\), so

\begin{equation*} \begin{aligned} \alpha&=\tfrac14(1+\alpha+2\beta)\;\Rightarrow\;3\alpha=1+2\beta,\\ \beta&=\tfrac14(2\alpha+\beta+0)\;\Rightarrow\;3\beta=2\alpha, \end{aligned} \end{equation*}

which gives \(\beta=\frac25\) and \(\alpha=\frac35\). Hence the probability is \(v_4=\frac25\).

\begin{equation*} R_{\mathrm{eff}}=\frac{(R_1+R_3)R_2}{R_1+R_2+R_3}. \end{equation*}

Apply a voltage \(V\) across the terminals. Kirchhoff’s law at the node between \(R_1\) and \(R_3\) forces the same current \(i_1\) through both resistors, and Ohm’s law on each branch gives

\begin{equation*} V=i_1R_1+i_1R_3,\qquad V=i_2R_2. \end{equation*}

Kirchhoff’s law at a terminal gives the total current \(i=i_1+i_2\), so

\begin{equation*} \begin{aligned} i&=V\left(\frac1{R_1+R_3}+\frac1{R_2}\right)=V\,\frac{R_1+R_2+R_3}{(R_1+R_3)R_2},\\ R_{\mathrm{eff}}&=\frac Vi=\frac{(R_1+R_3)R_2}{R_1+R_2+R_3}. \end{aligned} \end{equation*}

(1) \(v_c=\frac47\) and \(v_d=\frac37\). The conductances are \(c_{ac}=1\), \(c_{cb}=\frac12\), \(c_{cd}=1\), \(c_{ad}=\frac12\), \(c_{db}=1\). Kirchhoff’s law at \(c\) and at \(d\) gives

\begin{equation*} \begin{aligned} (1-v_c)+(v_d-v_c)-\tfrac12v_c=0&\;\Rightarrow\;\tfrac52v_c=1+v_d,\\ \tfrac12(1-v_d)+(v_c-v_d)-v_d=0&\;\Rightarrow\;\tfrac52v_d=\tfrac12+v_c, \end{aligned} \end{equation*}

whose solution is \(v_c=\frac47\), \(v_d=\frac37\).

(2) By Ohm’s law \(i_{xy}=(v_x-v_y)c_{xy}\):

Kirchhoff’s law checks: at \(c\), \(\frac37=\frac17+\frac27\); at \(d\), \(\frac27+\frac17=\frac37\).

(3) The current into \(a\) is \(i_a=\frac37+\frac27=\frac57\), so \(r_{\mathrm{eff}}=v_a/i_a=\frac75\).

(4) Take \(p_{xy}=c_{xy}/c_x\) (Section 4.5), with \(c_a=c_b=\frac32\) and \(c_c=c_d=\frac52\):

(5) \(\frac47\) from \(c\) and \(\frac37\) from \(d\): these are the voltages of part 1, by the probabilistic interpretation of voltages in Section 4.5. (Check: \(\frac47=\frac25+\frac25\cdot\frac37\).)

(6) \(\frac15\). The current \(i_{cd}\) is the net number of \(c\to d\) traversals when the current into \(a\) is one walk, so we rescale by \(i_a\):

\begin{equation*} \frac{i_{cd}}{i_a}=\frac{1/7}{5/7}=\frac15. \end{equation*}

(7) \(\frac{11}{21}\). By Section 4.5, \(p_{\mathrm{escape}}=c_{\mathrm{eff}}/c_a\), so

\begin{equation*} 1-p_{\mathrm{escape}}=1-\frac{5/7}{3/2}=\frac{11}{21} =\sum_yp_{ay}v_y=\frac23\cdot\frac47+\frac13\cdot\frac37. \end{equation*}

Set \(v_a=1\) and \(v_b=0\). The voltages are harmonic with boundary \(\{a,b\}\) (Section 4.5), so by the maximum and minimum principles of Section 4.5, \(0\le v_y\le1\) for every vertex \(y\). Then

\begin{equation*} c_{\mathrm{eff}}=\frac{i_a}{v_a}=i_a=\sum_y(1-v_y)\,c_{ay}\le\sum_yc_{ay}=c_a, \end{equation*}

so \(c_{\mathrm{eff}}/c_a\le1\).

Among all flows that carry a current \(i_a\) into the network at \(a\) and out at \(b\) while obeying Kirchhoff’s law, the one that minimizes \(E\) is the actual current \(i\), which is the unique one that also obeys Ohm’s law, \(i_{xy}=(v_x-v_y)/r_{xy}\). Here a flow is a function \(f\) on ordered adjacent pairs with \(f_{yx}=-f_{xy}\), and Kirchhoff’s law means \(\sum_y f_{xy}\) equals \(i_a\) at \(a\), \(-i_a\) at \(b\), and \(0\) elsewhere.

Let \(f\) be any such flow and put \(d=f-i\). Then \(d_{yx}=-d_{xy}\) and \(\sum_y d_{xy}=0\) at every vertex \(x\), \(a\) and \(b\) included. Expanding the square gives

\begin{equation*} \begin{aligned} E(f)&=E(i)+\sum_{x,y} i_{xy}d_{xy}r_{xy}+E(d),\\ \sum_{x,y} i_{xy}r_{xy}d_{xy}&=\sum_{x,y}(v_x-v_y)d_{xy} =\sum_x v_x\sum_y d_{xy}+\sum_y v_y\sum_x d_{yx}=0 . \end{aligned} \end{equation*}

The first equality in the second line is Ohm’s law, and the next uses \(d_{xy}=-d_{yx}\). Hence

\begin{equation*} E(f)=E(i)+\tfrac12\sum_{x,y}d_{xy}^2r_{xy}\ \ge\ E(i), \end{equation*}

with equality iff \(d\equiv0\), since every \(r_{xy}>0\). So the actual current is the unique minimizer of the energy dissipation. \(\blacksquare\)

Pass a unit current from \(a\) to \(b\). The current then dissipates exactly \(r_{\mathrm{eff}}\). Indeed, Ohm’s law and the Kirchhoff conditions from Exercise 4.29 give

\begin{equation*} E(i)=\tfrac12\sum_{x,y}i_{xy}(v_x-v_y)=\sum_x v_x\sum_y i_{xy}=v_a-v_b=r_{\mathrm{eff}} . \end{equation*}

The last step holds because \(r_{\mathrm{eff}}=(v_a-v_b)/i_a\) and \(i_a=1\).

Now let \(r^{\prime}\le r\) edgewise, and let \(i\) and \(i^{\prime}\) be the unit currents in the networks with resistances \(r\) and \(r^{\prime}\). The Kirchhoff conditions do not involve the resistances, so \(i\) is also an admissible unit flow for the network \(r^{\prime}\). Thomson’s principle (Exercise 4.29) applied in that network gives

\begin{equation*} r^{\prime}_{\mathrm{eff}}=E_{r^{\prime}}(i^{\prime})\le E_{r^{\prime}}(i)=\tfrac12\sum_{x,y}i_{xy}^2r^{\prime}_{xy} \le\tfrac12\sum_{x,y}i_{xy}^2r_{xy}=E_r(i)=r_{\mathrm{eff}} . \end{equation*}

So reducing a resistor (or any set of resistors) cannot increase the effective resistance. Swapping the roles of \(r\) and \(r^{\prime}\) gives the second claim: if \(r\le r^{\prime}\), then \(r_{\mathrm{eff}}\le r^{\prime}_{\mathrm{eff}}\), so increasing a resistor cannot decrease the effective resistance. \(\blacksquare\)

\(h_{uv}=n-1\) on the cycle, and \(h_{uv}=(n-1)^2\) once the edge \((u,v)\) is removed.

The cycle has \(m=n\) unit resistors. Between \(u\) and \(v\) the edge \(uv\) (resistance \(1\)) is in parallel with the path of the other \(n-1\) edges (resistance \(n-1\)), so

\begin{equation*} r_{uv}=\frac{1\cdot(n-1)}{1+(n-1)}=\frac{n-1}{n},\qquad \operatorname{commute}(u,v)=2m\,r_{uv}=2(n-1) \end{equation*}

by Theorem 4.9 (the cycle is connected and all its edges have unit weight). The reflection of the cycle that swaps \(u\) and \(v\) is a graph automorphism, so \(h_{uv}=h_{vu}\), and each equals \(n-1\).

Removing \((u,v)\) leaves a path on \(n\) vertices with \(u\) and \(v\) at its two ends. Lemma 4.7, whose proof computes \(h_{1n}\) exactly, gives \(h_{uv}=(n-1)^2\).

\(h_{uv}=3,\ 5,\ 7\) for graphs (i), (ii) and (iii).

In each graph \(v\) hangs off \(u\) by a single edge, and the rest of the graph (all vertices except \(v\)) is connected, with \(m=1,2,3\) edges. Exercise 4.34 therefore applies and gives

\begin{equation*} h_{uv}=2m+1=3,\quad 5,\quad 7 . \end{equation*}

\(h_{24}=4\) in (a), \(8\) in (b) and \(2\) in (c). Going from (a) to (b) adds the edge \(\{1,2\}\) and increases \(h_{24}\). Going from (a) to (c) adds the edge \(\{2,4\}\) and decreases it.

A walk from \(2\) never leaves its connected component, so the isolated vertex \(1\) is irrelevant in (a) and (c).

(a) The walk runs on the path \(2-3-4\) of \(3\) vertices from one end to the other. Lemma 4.7 gives \(h_{24}=(3-1)^2=4\).

(b) The proof of Lemma 4.7 gives \(h_{i,i+1}=2i-1\) on the path \(1-2-3-4\), so

\begin{equation*} h_{24}=h_{23}+h_{34}=3+5=8 . \end{equation*}

(c) The walk runs on a triangle, the cycle of length \(3\), and \(2,4\) are adjacent on it. Exercise 4.31 gives \(h_{24}=3-1=2\).

The whole graph has \(m+1\) edges, and the walk from \(v\) moves to \(u\) on its first step, so \(h_{vu}=1\). Apply Theorem 4.9; its hypotheses hold because the graph is connected and undirected with unit resistances. The edge \((u,v)\) is the only edge at \(v\), so all current between \(u\) and \(v\) flows through it and \(r_{uv}=1\). Therefore

\begin{equation*} \begin{aligned} h_{uv}&=\operatorname{commute}(u,v)-h_{vu}=2(m+1)\,r_{uv}-1\\ &=2(m+1)-1=2m+1 . \end{aligned} \end{equation*}

\(\blacksquare\)

\(n-1\) steps, starting from any other vertex.

At each step the walk on \(K_n\) moves to a uniformly random one of the \(n-1\) vertices other than its current one. So whenever it is not yet at the target \(y\), its next step hits \(y\) with probability \(\frac1{n-1}\), independently of the past. The hitting time is therefore geometric with success probability \(\frac1{n-1}\), and

\begin{equation*} h_{xy}=\sum_{t\ge1}t\Bigl(1-\tfrac1{n-1}\Bigr)^{t-1}\tfrac1{n-1}=n-1 \qquad (x\ne y). \end{equation*}

The general solution is \(t(i) = c_1 2^i + c_2 3^i\) with \(c_1, c_2\) arbitrary constants, and two boundary conditions make it unique.

Trying \(t(i) = \lambda^i\) gives the characteristic equation

\begin{equation*} \lambda^{i}\left(\lambda^2 - 5\lambda + 6\right) = \lambda^{i}(\lambda - 2)(\lambda - 3) = 0, \end{equation*}

so \(2^i\) and \(3^i\) are solutions, and by linearity so is every combination \(c_1 2^i + c_2 3^i\). Conversely, the recurrence \(t(i+2) = 5t(i+1) - 6t(i)\) fixes the whole sequence once \(t(0), t(1)\) are given, so the solution space has dimension \(2\). The two independent solutions \(2^i, 3^i\) therefore span it.

Two values at distinct indices \(i \ne j\) pin down \(c_1, c_2\), because

\begin{equation*} \begin{aligned} \det\begin{pmatrix} 2^i & 3^i \\ 2^j & 3^j \end{pmatrix} &= 2^i 3^j - 3^i 2^j \\ &= 2^i 2^j\left(\left(\tfrac32\right)^j - \left(\tfrac32\right)^i\right) \ne 0 . \end{aligned} \end{equation*}

For example, \(t(0) = a\) and \(t(1) = b\) give \(t(i) = (3a - b)\,2^i + (b - 2a)\,3^i\). So \(\boxed{t(i) = c_1 2^i + c_2 3^i,\ \text{two boundary conditions}}\).

The answers are: (1) \(t(i) = \sum_{j=1}^{k} c_j \lambda_j^i\); (2) \(t(i) = \sum_j p_j(i)\lambda_j^i\) with \(\deg p_j < m_j\); (3) \(k\); (4) \(f\) plus the homogeneous solutions. The book prints \(a_{k-i}t^{k-1}\), which should read \(a_{k-1}t^{k-1}\). Throughout, write \(p(\lambda) = \sum_{s=0}^k a_s \lambda^s\) and assume \(a_k \ne 0\) and \(a_0 \ne 0\), so the order really is \(k\) and every root is nonzero.

(1) “Distinct roots” means \(p\) has \(k\) distinct roots \(\lambda_1,\dots,\lambda_k\) (so \(r = k\)). Each \(\lambda_j^i\) is a solution because \(\sum_s a_s \lambda_j^{i+s} = \lambda_j^i p(\lambda_j) = 0\). They are linearly independent: the Vandermonde matrix \((\lambda_j^{i})_{0\le i<k,\,1\le j\le k}\) is nonsingular. By the dimension count in (3) they form a basis, so

\begin{equation*} \boxed{t(i) = c_1\lambda_1^i + c_2\lambda_2^i + \cdots + c_k\lambda_k^i .} \end{equation*}

(2) Suppose the distinct roots are \(\lambda_1,\dots,\lambda_r\) with multiplicities \(m_1+\cdots+m_r = k\). If \(\lambda\) is a root of multiplicity \(m\) and \(0 \le l < m\), then \(i^l\lambda^i\) is a solution. With \(D = d/d\lambda\), \(\sum_s a_s s^q\lambda^s = ((\lambda D)^q p)(\lambda)\), which vanishes for \(q < m\). Hence

\begin{equation*} \begin{aligned} \sum_{s} a_s (i+s)^l \lambda^{i+s} &= \lambda^i \sum_{q=0}^{l}\binom{l}{q} i^{\,l-q}\sum_s a_s s^q\lambda^s \\ &= \lambda^i \sum_{q=0}^{l}\binom{l}{q} i^{\,l-q}\,\big((\lambda D)^q p\big)(\lambda) = 0 . \end{aligned} \end{equation*}

These \(k\) sequences \(i^l\lambda_j^i\) (\(0 \le l < m_j\)) are linearly independent. Their generating functions \(\sum_{i\ge0} i^l\lambda_j^i z^i\) are rational, with a pole of order exactly \(l+1\) at \(z = 1/\lambda_j\) and no other pole, so a nontrivial combination cannot be identically \(0\). By the dimension count in (3) they form a basis:

\begin{equation*} \boxed{t(i) = \sum_{j=1}^{r} \big(c_{j,0} + c_{j,1} i + \cdots + c_{j,m_j-1} i^{m_j-1}\big)\lambda_j^i .} \end{equation*}

(3) Since \(a_k \ne 0\), the recurrence \(t(i+k) = -a_k^{-1}\sum_{s<k} a_s t(i+s)\) determines the sequence from \(t(0), \dots, t(k-1)\), and those \(k\) values can be anything. The map \(t \mapsto (t(0),\dots,t(k-1))\) is therefore a linear bijection onto \(\mathbb{R}^k\), so the solution space has \(\boxed{\dim = k}\).

(4) If \(L t = g\) and \(L f = g\), where \(L\) is the linear difference operator, then \(L(t - f) = 0\). So the full solution set is the affine space

\begin{equation*} \boxed{\{\, f + h : h \text{ solves the homogeneous equation} \,\}}, \end{equation*}

that is, \(f(i)\) plus the general form from (2). It has dimension \(k\).

Use Theorem 4.9, \(\operatorname{commute}(x,y) = 2m\,r_{xy}\). One added edge can lower \(r_{xy}\) a lot, or leave \(r_{xy}\) unchanged while raising \(m\).

(i) Decrease. Take the path \(x - z - y\), so \(m = 2\) and \(r_{xy} = 2\). Adding the edge \(\{x,y\}\) gives a triangle with \(m = 3\) and \(r_{xy} = 1 \parallel 2 = \tfrac23\):

\begin{equation*} \operatorname{commute}(x,y) = 2\cdot 2\cdot 2 = 8 \;\longrightarrow\; 2\cdot 3\cdot \tfrac23 = 4 . \end{equation*}

(ii) Increase. Take the path \(x - y - z - w\), so \(m = 3\) and \(r_{xy} = 1\). Adding the edge \(\{y,w\}\) keeps \(\{x,y\}\) a bridge. All current from \(x\) to \(y\) still flows through that one resistor, so \(r_{xy} = 1\), but now \(m = 4\):

\begin{equation*} \operatorname{commute}(x,y) = 2\cdot 3\cdot 1 = 6 \;\longrightarrow\; 2\cdot 4\cdot 1 = 8 . \end{equation*}

Both graphs are connected before and after the edge is added, as Theorem 4.9 requires.

The answers are (1) \(n\) and (2) \(nH_n = n\sum_{k=1}^n \frac1k \approx n\ln n\). This is the coupon collector problem of Section 4.6.

(1) Each draw is \(1\) independently with probability \(p = 1/n\), so the number of draws \(T\) is geometric with mean \(1/p\):

\begin{equation*} \mathbb{E}[T] = \sum_{t\ge1} t\,(1-p)^{t-1}p = \frac1p = \boxed{n}. \end{equation*}

(2) Let \(T_k\) be the number of draws taken while exactly \(k-1\) distinct integers have been seen. Each of those draws is new with probability \(\frac{n-k+1}{n}\), so \(T_k\) is geometric with mean \(\frac{n}{n-k+1}\). By linearity,

\begin{equation*} \mathbb{E}\Big[\sum_{k=1}^n T_k\Big] = \sum_{k=1}^{n}\frac{n}{n-k+1} = n\sum_{j=1}^{n}\frac1j = \boxed{nH_n = n\ln n + O(n)}. \end{equation*}

In every case the return time to \(A\) is \(1/\pi_A = 2m/d_A = 2m\), because \(d_A = 1\). The return time is the reciprocal of the stationary probability (Section 4.8), and \(\pi_x = d_x/2m\) for a walk on a connected undirected graph (Section 4.4, via Lemma 4.3). Each graph is connected and has \(d_A = 1\).

(a) \(m = n-1\), so the return time is \(\boxed{2(n-1) = 2m}\).

(b) \(m = 1 + (n-2) = n-1\), so the return time is \(\boxed{2(n-1) = 2m}\).

(c) \(m = 1 + \binom{n-1}{2}\), so the return time is

\begin{equation*} \boxed{2m = 2 + (n-1)(n-2) = n^2 - 3n + 4}. \end{equation*}

The return time is \(\boxed{2m}\), assuming the replacement graph is connected so that the whole walk is connected.

The whole graph has \(m\) edges and \(d_A = 1\). As in Exercise 4.40, the return time is \(1/\pi_A = 2m/d_A = 2m\).

Method (2): the first step goes from \(A\) to \(B\). By Exercise 4.34, applied with \(m-1\) edges in the attached subgraph, \(h_{BA} = 2(m-1)+1\). So the return time is \(1 + h_{BA} = 2m\).

The expected length is \(\boxed{2m}\), assuming the whole graph is connected.

The whole graph has \(m\) edges and \(d_A = d\), so one return to \(A\) takes expected time \(1/\pi_A = 2m/d\) (Sections 4.4 and 4.8). Write \(R_1, \dots, R_d\) for the lengths of the successive excursions from \(A\) back to \(A\). By the Markov property at each return, every \(R_j\) has the law of a single return time. Hence

\begin{equation*} \mathbb{E}\Big[\sum_{j=1}^{d} R_j\Big] = d\cdot\frac{2m}{d} = 2m . \end{equation*}

The eigenvalues are \(0\) and \(2\). Let the component be the edge \(\{u,v\}\), so \(d_u = d_v = 1\) and \(u,v\) have no other neighbours. Ordering \(u,v\) first, \(L\) is block diagonal with the block

\begin{equation*} L_{\{u,v\}} = \begin{pmatrix} 1 & -1 \\ -1 & 1 \end{pmatrix}, \qquad L(e_u + e_v) = 0, \qquad L(e_u - e_v) = 2\,(e_u - e_v), \end{equation*}

since rows of \(L\) outside \(\{u,v\}\) vanish on vectors supported on \(\{u,v\}\). Hence \(0\) and \(2\) are eigenvalues of \(L\), with eigenvectors \(e_u + e_v\) and \(e_u - e_v\).

His measure assigns importance \(0\) to every edge. For the random walk on an undirected graph with \(m\) edges the stationary probability is \(p_i = d_i/2m\) (Section 4.4, via Lemma 4.3), so for every edge \((i,j)\)

\begin{equation*} \frac{1}{d_i}\, p_i = \frac{1}{2m} = \frac{1}{d_j}\, p_j, \qquad\text{hence}\qquad \left| \frac{1}{d_i} p_i - \frac{1}{d_j} p_j \right| = 0 . \end{equation*}

(On a disconnected graph \(p_i = \lambda_c d_i/2m_c\) on each component \(c\), and the same cancellation holds.) This is detailed balance: in stationarity every edge is traversed equally often in both directions, so his net flow carries no information about the edge.

Let \(X_t, Y_t\) be the walks (both stepping at each time \(t\)) and \(D_t = X_t - Y_t\); they meet at time \(t \ge 1\) iff \(D_t = 0\). Each increment of \(D\) is \(\xi - \eta\) with \(\xi,\eta\) independent and uniform on \(\{\pm e_1, \pm e_2\}\), so \(D_t\) lies in \(\{(x,y) : x+y \text{ even}\}\), and the bijection \((x,y) \mapsto \big(\tfrac{x+y}{2}, \tfrac{x-y}{2}\big)\) onto \(\mathbb{Z}^2\) turns \(D\) into a walk \(D^{\prime}\) on \(\mathbb{Z}^2\) with

\begin{equation*} P(\text{stay}) = \tfrac{4}{16}, \qquad P(\pm e_1) = P(\pm e_2) = \tfrac{2}{16}, \qquad P(\pm e_1 \pm e_2) = \tfrac{1}{16}. \end{equation*}

Deleting the holding steps (geometric, hence a.s. finite) leaves the random walk on the electrical network \(\mathbb{Z}^2\) with conductance \(2\) on lattice edges and \(1\) on the diagonal edges of each square (transition probabilities \(2/12\) and \(1/12\), Section 4.5); so \(c_a = 12\) at every vertex and it suffices that the escape probability \(p_{\text{escape}} = c_{\text{eff}}/c_a\) of this network tends to \(0\) (Section 4.5).

Let \(S_k = \{z : \|z\|_\infty = k\}\) and short each \(S_k\), which can only decrease \(r_{\text{eff}}\) (Section 4.7). Every edge joins shells differing by at most \(1\), so the shorted network is a series chain \(S_0, S_1, \dots, S_n = b\) whose \(k\)-th link is the conductance \(C_k\) of all edges from \(S_k\) to \(S_{k+1}\). On the side of the square \(S_k\) a non-corner vertex has one outward lattice edge and two outward diagonals, and each of the \(4\) corners has two and three, so

\begin{equation*} \begin{aligned} C_k &= 2\big[4(2k-1) + 4\cdot 2\big] + \big[4(2k-1)\cdot 2 + 4\cdot 3\big] = 32k + 12, \\ r_{\text{eff}} &\ge \sum_{k=0}^{n-1} \frac{1}{32k+12} \longrightarrow \infty \quad (n \to \infty), \end{aligned} \end{equation*}

(also correct for \(k = 0\)). Hence \(p_{\text{escape}} = 1/(12\, r_{\text{eff}}) \to 0\), so \(D^{\prime}\), and therefore \(D\), returns to \(0\) with probability one: the two walks meet with probability one. \(\blacksquare\)

Yes, with probability one. Let \(A_t, B_t\) be the two counts after \(t\) flips each and \(D_t = (A_t - B_t)/2 \in \mathbb{Z}\). Each round changes \(A_t - B_t\) by \(0\) (probability \(\tfrac12\)) or by \(\pm 2\) (probability \(\tfrac14\) each), so

\begin{equation*} P(D_{t+1} - D_t = 0) = \tfrac12, \qquad P(D_{t+1} - D_t = \pm 1) = \tfrac14 . \end{equation*}

Removing the holding steps (each run of them is geometric, hence a.s. finite) leaves the simple random walk on \(\mathbb{Z}\), whose escape probability \(1/(2 r_{\text{eff}})\) tends to \(0\) because \(r_{\text{eff}} = n/2 \to \infty\) (Section 4.7, one dimension). So \(D\) returns to \(0\) with probability one, i.e. at some time \(t \ge 1\) the counts agree.

The escape probability is \(0\): the walk returns to the origin with probability one. Every vertex has degree \(c_a = 8\) (four lattice and four diagonal unit resistors), so \(p_{\text{escape}} = 1/(8\, r_{\text{eff}})\) (Section 4.5). Short each square shell \(S_k = \{z : \|z\|_\infty = k\}\), which only decreases \(r_{\text{eff}}\) (Section 4.7); since every edge joins shells differing by at most \(1\), this gives a series chain \(S_0, S_1, \dots, S_n\) whose \(k\)-th link consists of \(C_k\) unit resistors in parallel. A non-corner vertex of \(S_k\) has \(3\) edges into \(S_{k+1}\) and each of the \(4\) corners has \(5\), so

\begin{equation*} \begin{aligned} C_k &= 3 \cdot 4(2k-1) + 5 \cdot 4 = 24k + 8 \qquad (k \ge 0), \\ r_{\text{eff}} &\ge \sum_{k=0}^{n-1} \frac{1}{24k+8} \ge \frac{1}{24} \ln n \longrightarrow \infty . \end{aligned} \end{equation*}

Hence \(p_{\text{escape}} \le 3/\ln n \to 0\).

The simulation gives \(p_{\text{escape}} \approx 0.66\). Following Section 4.5, the escape probability for the box \(\{\|z\|_\infty < n\}\) is the probability that a walk from the origin reaches \(\|z\|_\infty = n\) before returning; the program estimates it from \(10^5\) walks for each \(n\) (seed \(0\)).

import numpy as np
rng = np.random.default_rng(0)
steps = np.vstack([np.eye(3, dtype=int), -np.eye(3, dtype=int)])
def escape(n, walks):
    pos = steps[rng.integers(0, 6, walks)]       # first step away from 0
    alive = np.ones(walks, bool); esc = 0
    while alive.any():
        idx = np.flatnonzero(alive)
        pos[idx] += steps[rng.integers(0, 6, idx.size)]
        out = np.abs(pos[idx]).max(axis=1) >= n  # reached the boundary
        back = ~pos[idx].any(axis=1)             # returned to the origin
        esc += out.sum(); alive[idx[out | back]] = False
    return esc / walks
for n in [2, 5, 10, 20, 40, 80]:
    print(n, escape(n, 100000))

The estimates (standard error about \(0.0015\) each) are

Solving the harmonic equations exactly for small boxes gives \(p_{\text{escape}}(2) = 17/22 \approx 0.7727\) and \(p_{\text{escape}}(5) \approx 0.6982\), within one standard error of the simulation. The estimates decrease in \(n\) and level off near \(0.66\) (the classical value is \(1 - 0.3405 = 0.6595\)), inside the bounds \(\tfrac16 \le p_{\text{escape}} \le \tfrac56\) of Section 4.7.

\(p_{\text{escape}} = \tfrac12\). Make each edge a unit resistor and, as in Section 4.5, merge all vertices at depth \(n\) into \(b\). Level \(k\) has \(2^k\) vertices, and by symmetry all vertices of one level have the same voltage, so shorting each level leaves \(r_{\text{eff}}\) unchanged; the link between levels \(k-1\) and \(k\) is then \(2^k\) unit resistors in parallel. With \(c_a = 2\) at the root,

\begin{equation*} \begin{aligned} r_{\text{eff}}(n) &= \sum_{k=1}^{n} 2^{-k} = 1 - 2^{-n}, \\ p_{\text{escape}}(n) &= \frac{c_{\text{eff}}}{c_a} = \frac{1}{2\,(1 - 2^{-n})} \longrightarrow \frac12 . \end{aligned} \end{equation*}

Method (2): after the first step the depth performs a walk that moves down with probability \(\tfrac23\) and up with probability \(\tfrac13\), and by gambler’s ruin it ever reaches depth \(0\) from depth \(1\) with probability \((1/3)/(2/3) = \tfrac12\); so the walk never returns with probability \(\tfrac12\).

The escape probability is \(\boxed{1/2}\).

The walk is the random walk on the path \(0 - 1 - 2 - \cdots\) with conductance \(c_{i,i+1} = 2^i\) on edge \(\{i, i+1\}\): at \(i \ge 1\) this gives

\begin{equation*} p_{i,i+1} = \frac{2^i}{2^{i-1} + 2^i} = \frac23, \qquad p_{i,i-1} = \frac13, \end{equation*}

and at \(0\) the only edge forces \(p_{01} = 1\), with \(c_0 = c_{01} = 1\). As in Section 4.5, merge all vertices at distance \(\ge d\) from \(0\) into one vertex \(b\). The resistors \(1/c_{i,i+1} = 2^{-i}\) are in series, so

\begin{equation*} r_{\mathrm{eff}}(d) = \sum_{i=0}^{d-1} 2^{-i} = 2\left(1 - 2^{-d}\right) \to 2 , \end{equation*}

and the formula \(p_{\mathrm{escape}} = c_{\mathrm{eff}}/c_a\) of Section 4.5 gives

\begin{equation*} p_{\mathrm{escape}} = \lim_{d \to \infty} \frac{1}{c_0\, r_{\mathrm{eff}}(d)} = \frac{1}{1 \cdot 2} = \frac12 . \end{equation*}

Method (2): the first step goes to \(1\), and from \(1\) a walk with up-probability \(p = 2/3\) ever reaches \(0\) with probability \((q/p)^1 = 1/2\) (gambler’s ruin against an infinitely rich opponent), so it escapes with probability \(1 - 1/2 = 1/2\).

Undirected: \(\pi = \left(\pi_A, \pi_B, \pi_C, \pi_D, \pi_E\right) = \left(\tfrac16, \tfrac16, \tfrac14, \tfrac14, \tfrac16\right)\), every edge carries flow \(\tfrac1{12}\) in each direction. Directed: \(\pi = \left(\pi_A, \pi_B, \pi_C, \pi_D\right) = \left(\tfrac27, \tfrac17, \tfrac27, \tfrac27\right)\).

Undirected graph. The degrees are \(d_A = d_B = d_E = 2\), \(d_C = d_D = 3\), with \(m = 6\) edges, so by Section 4.4 (Lemma 4.3 applied to \(\pi_x p_{xy} = \pi_y p_{yx}\)) \(\pi_x = d_x/2m = d_x/12\). The flow along \(x \to y\) is

\begin{equation*} \pi_x p_{xy} = \frac{d_x}{12}\cdot\frac1{d_x} = \frac1{12} = \pi_y p_{yx} \qquad \text{for every edge } \{x,y\}. \end{equation*}

The condition is detailed balance: the flow through each edge is the same in both directions (so the net flow on every edge is \(0\)), and here it is even the same, \(1/2m\), on every edge.

Directed graph. The out-degrees are \(2, 1, 1, 1\), and \(\pi = \pi P\) reads

\begin{equation*} \begin{aligned} \pi_A &= \pi_D, & \pi_B &= \tfrac12\pi_A, \\ \pi_C &= \tfrac12\pi_A + \pi_B, & \pi_D &= \pi_C, \end{aligned} \end{equation*}

so \(\pi_A = \pi_C = \pi_D = 2\pi_B\) and normalising gives \(\pi = (2,1,2,2)/7\) (unique by Theorem 4.2, the graph being strongly connected). The edge flows \(\pi_x p_{xy}\) are

Edges carry flow in one direction only, so detailed balance fails; the condition is conservation of flow (Kirchhoff’s current law): at every vertex the total flow in equals the total flow out, both equal to \(\pi_x\) (e.g. at \(C\): \(1/7 + 1/7 = 2/7\)). This conservation also holds in the undirected graph, where it is implied by detailed balance.

A star of \(k\) pages raises the pagerank of its centre by almost exactly \(k\,r(1-r)/N\) (\(N = 200 + k\)), which is a factor of about \(30\) for a low-pagerank centre but only about \(5.5\) for a high-pagerank one at \(k = 100\), \(r = 0.15\); a self-loop multiplies pagerank by only \(1/\left(1 - \frac{1-r}{d+1}\right) \approx 1.06\) to \(1.13\) for out-degree \(d \approx 8\).

The program builds \(G(200, 8/199)\) with seed \(0\) (mean out-degree \(7.79\)), keeps the \(k\) new vertices in both versions (without the star they have no out-edges and so restart with probability \(1\), as in Section 4.8), and computes pagerank by the power method. The star centre \(v\) is taken to be the vertex of lowest pagerank (\(200\,\pi_v = 0.32\), in-degree \(3\), out-degree \(8\)) and of highest pagerank (\(200\,\pi_v = 2.07\), in-degree \(12\), out-degree \(7\)).

import numpy as np
rng = np.random.default_rng(0); n = 200
A0 = (rng.random((n, n)) < 8/199).astype(float); np.fill_diagonal(A0, 0)
def pagerank(A, r):                      # power method, restart prob r
    N = len(A); d = A.sum(1); P = np.full((N, N), 1/N)   # dangling: restart
    P[d > 0] = (1-r)*A[d > 0]/d[d > 0, None] + r/N
    p = np.ones(N)/N
    for _ in range(3000): p = p @ P
    return p
def graph(k, v, star, loop=False):       # k extra vertices, optional star/loop
    A = np.zeros((n+k, n+k)); A[:n, :n] = A0
    if star: A[n:, v] = 1
    if loop: A[v, v] += 1
    return A
pr = pagerank(A0, .15); lo, hi = pr.argmin(), pr.argmax()
for v in (lo, hi):
    for r in (.05, .15, .3, .5):
        gain = [pagerank(graph(k, v, 1), r)[v] / pagerank(graph(k, v, 0), r)[v]
                for k in (1, 5, 10, 20, 50, 100)]
        loop = pagerank(graph(0, v, 0, 1), r)[v] / pagerank(A0, r)[v]
        print(v, r, np.round(gain, 2), round(loop, 3))

The ratio \(\pi_v(\text{star})/\pi_v(\text{no star})\) came out as follows.

Each star vertex \(u\) has no in-edges, so \(\pi_u = r/N\) (measured: \(N\pi_u = 0.150\) at \(r = 0.15\)), and it pushes \((1-r)\pi_u\) into \(v\) per step:

\begin{equation*} \begin{aligned} &\Delta \pi_v \approx \frac{k\,r(1-r)}{N}, \\ &k = 10,\ r = 0.15:\quad k\,r(1-r) = 1.275, \\ &N\,\Delta\pi_v = 1.27 \text{ (low)},\quad 1.29 \text{ (high)} . \end{aligned} \end{equation*}

The absolute gain is thus the same for every centre, so the relative gain is large exactly when \(\pi_v\) is small, as the exercise warns; it grows linearly in \(k\), vanishes as \(r \to 0\) (the countermeasure of Section 4.8: rarer restarts) and per star vertex is largest near \(r = 1/2\). Other vertices change much less: at \(r = 0.15\) the out-neighbours of the low centre rise by a median \(10\%\) (\(k = 10\)) and \(80\%\) (\(k = 100\)), while the remaining vertices fall by a median \(1.5\%\) and \(9\%\), since the \(k\) restart sources that spread mass uniformly have become feeders of \(v\).

For the self-loop, the loop takes a share \(1/(d+1)\) of \(v\)’s out-edges, so with in-flow \(I\) essentially unchanged, \(\pi_v = I + (1-r)\pi_v/(d+1)\) and

\begin{equation*} \begin{aligned} &\frac{\pi_v(\text{loop})}{\pi_v} \approx \frac{1}{1 - \frac{1-r}{d+1}}, \\ &d = 8:\ 1.118,\ 1.104,\ 1.084,\ 1.059 \text{ for } r = 0.05, 0.15, 0.3, 0.5, \end{aligned} \end{equation*}

matching the low-centre column exactly; for \(d = 1\) the same formula gives the book’s factor \(1/(1 - 0.85/2) = 1.74\).

The star cuts the hitting time of its centre roughly as \(1/h \approx 1/h_0 + k\,r(1-r)/N\) (at \(r = 0.15\), \(k = 100\): \(672 \to 22.8\) for the low-pagerank centre, \(106 \to 19.3\) for the high one), whereas a self-loop at the centre leaves the hitting time exactly unchanged.

Here \(h_v = \frac1N\sum_x h_{xv}\) is the hitting time of Section 4.8, computed with restarts from the linear system \(h = \mathbf 1 + Q h\), where \(Q\) is the transition matrix with row and column \(v\) deleted; the graph, the \(k\) extra vertices and the two centres are those of Exercise 4.52 (same seed, reusing its A0 and graph).

def hitting(A, r, v):                  # mean over starts of E[time to hit v]
    N = len(A); d = A.sum(1); P = np.full((N, N), 1/N)
    P[d > 0] = (1-r)*A[d > 0]/d[d > 0, None] + r/N
    keep = np.arange(N) != v
    h = np.linalg.solve(np.eye(N-1) - P[np.ix_(keep, keep)], np.ones(N-1))
    return np.r_[h, 0].mean(), h[:n-1].mean()   # all N starts / original x != v
for v in (lo, hi):
    for r in (.05, .15, .3, .5):
        print(v, r, [(hitting(graph(k, v, 0), r, v), hitting(graph(k, v, 1), r, v))
                     for k in (1, 10, 100)],
              hitting(graph(0, v, 0, 1), r, v)[0] - hitting(A0, r, v)[0])

Entries are \(h_v\) without the star / with the star (averaged over all \(N\) starts) / with the star (averaged over the \(199\) original vertices other than \(v\)).

Without the star, a walk from the original vertices reaches \(v\) at rate about \(1/h_0\); the star adds a second route, since each step restarts with probability \(r\), lands in the star with probability \(k/N\) and then steps to \(v\) with probability \(1-r\). Adding the rates,

\begin{equation*} \begin{aligned} &\frac1{h} \approx \frac1{h_0} + \frac{k\,r(1-r)}{N}, \\ &\text{low}, r = 0.15, k = 100:\ \Big(\tfrac1{672.1} + \tfrac{0.1275}{3}\Big)^{-1} = 22.7 , \end{aligned} \end{equation*}

against the measured \(22.8\); the same formula gives \(59.0, 13.9, 11.7\) (low, \(r = 0.05, 0.3, 0.5\)) and \(37.7, 19.3, 12.8, 11.2\) (high) against \(59.3, 14.0, 11.7\) and \(37.7, 19.3, 12.9, 11.2\). So, as for pagerank in Exercise 4.52, the star helps a poorly linked centre far more (low, \(r = 0.15\): hitting time divided by \(29\), pagerank multiplied by \(30\)) than a well linked one (divided by \(5.5\); pagerank \(\times 5.5\)), and its effect grows with \(k\) and with \(r(1-r)\). The self-loop difference was \(0\) in every run, as it must be: a loop only changes row \(v\) of the transition matrix, and that row is deleted from the hitting-time system. This is the contrast Section 4.8 uses to detect spam: short cycles at \(v\) raise pagerank without lowering hitting time, while a star moves both.

By at most a factor \(\boxed{1/\left(1-(1-r)^2\right) = 1/\left(r(2-r)\right)}\), which is \(3.60\) at \(r = 0.15\), however many loops are added; any further gain comes from the new pages acting as a star, not from the loops.

Let page \(v\) have \(d\) out-links and receive pagerank flow \(J\) per step from the rest of the graph (restarts included), held fixed as in the self-loop calculation of Section 4.8. Add new pages \(u_1, \dots, u_m\) with edges \(v \to u_i\) and \(u_i \to v\). Each \(u_i\) is entered only by restart or from \(v\), so

\begin{equation*} \begin{aligned} \pi_{u_i} &= \frac rN + \frac{(1-r)\,\pi_v}{d+m}, \\ \pi_v &= J + (1-r)\sum_{i=1}^m \pi_{u_i} = J + \frac{m\,r(1-r)}{N} + (1-r)^2\frac{m}{d+m}\,\pi_v , \end{aligned} \end{equation*}

hence

\begin{equation*} \pi_v = \frac{J + m\,r(1-r)/N}{1 - (1-r)^2\,\frac{m}{d+m}} . \end{equation*}

The numerator is the pagerank \(v\) would get if the \(u_i\) only pointed to \(v\) (a star, Exercise 4.52); the loops multiply it by

\begin{equation*} \frac{1}{1 - (1-r)^2\frac{m}{d+m}} \ \le\ \frac1{1-(1-r)^2} = \frac{1}{r(2-r)} , \end{equation*}

with equality when \(v\) drops its \(d\) original out-links (the capture of Section 4.8) and in the limit \(m \to \infty\) otherwise. For \(r = 0.05, 0.15, 0.5\) the cap is \(10.26, 3.60, 1.33\), roughly half the factor \(1/r\) that unlimited self-loops would give (in Section 4.8, \(k \to \infty\) self-loops turn \(\pi_i = 0.85\,\pi_j p_{ji}\) into \(5.67\,\pi_j p_{ji}\), a factor \(6.67\)), since each excursion now costs two steps that must both avoid a restart. On the graph of Exercise 4.52 (low-pagerank centre, \(d = 8\), \(r = 0.15\)) the measured loop factor (against the plain star) was \(1.384\) for \(m = 5\) against the formula’s \(1.385\), and \(3.593\) against \(3.604\) once the original out-links were deleted.

(1) \(h_{1j} = n\) for every \(j \ne 1\) (and \(h_{11} = 0\)); (2) \(h_{1j} = j - 1\). The two notes are answered after these, the first one being false as printed.

(1) With self-loops each vertex has out-degree \(n\), so every step goes to a uniformly random vertex, independently of the past. The time \(T\) to reach \(j \ne 1\) is therefore geometric with success probability \(1/n\):

\begin{equation*} h_{1j} = \mathbb E\, T = \sum_{t \ge 1} t\left(1 - \tfrac1n\right)^{t-1}\tfrac1n = n . \end{equation*}

(2) On the cycle \(1 \to 2 \to \cdots \to n \to 1\) the walk is deterministic, reaching \(j\) at time exactly \(j - 1\), so \(h_{1j} = j - 1\), the largest being \(h_{1n} = n - 1\).

Full rank versus strong connectivity. Neither implies the other, so the note “full rank implies strongly connected” is false. The printed matrix is the adjacency matrix of the undirected path \(1 - 3 - 2\), which is strongly connected, yet rows \(1\) and \(2\) coincide, so its rank is \(2 < 3\) (its eigenvalues are \(\pm\sqrt2, 0\)). Conversely, two disjoint \(2\)-cycles have adjacency matrix

\begin{equation*} \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \oplus \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \qquad \det = 1 , \end{equation*}

which is of full rank, but the graph is not strongly connected.

Aperiodicity. Let \(P\) be the walk on a strongly connected graph, with eigenvalues ordered so that \(1 = \lambda_1 \ge |\lambda_2| \ge \cdots\); the proof of Lemma 4.1 shows the eigenspace of \(\lambda_1 = 1\) is spanned by \(\mathbf 1\). By the Perron-Frobenius theorem for irreducible nonnegative matrices, \(\lambda_1 = 1\) is in fact a simple root, and if the graph has period \(p\) then the eigenvalues of modulus \(1\) are exactly the \(p\)-th roots of unity \(e^{2\pi i k/p}\), \(k = 0, \dots, p-1\), each simple. Hence

\begin{equation*} \text{aperiodic} \iff p = 1 \iff \lambda_1 > |\lambda_2| . \end{equation*}

Read with the real ordering \(\lambda_1 > \lambda_2\), the statement is false: \(\lambda_1 = 1\) being simple, \(\lambda_1 > \lambda_2\) always holds, yet a single undirected edge has eigenvalues \(1, -1\) and period \(2\). The modulus is essential.

Parse the HTML into tags (not with a regular expression), take the href of every <a> and <area> tag (and the src of <frame> / <iframe>), resolve each against the page’s base URL, strip the #fragment, keep only http(s) links, and deduplicate; these URLs become the crawl’s out-edges of the page.

The page source shows hyperlinks as <a href=...>anchor text</a>; the details a crawler must get right are:

• Parsing: tag and attribute names are case-insensitive, attribute values may use single, double or no quotes, and real pages are malformed, so a tolerant HTML tokenizer is used rather than pattern matching on the text.
• Resolution: relative links (intro.html, ../index.html, /docs/, //host/path) are joined to the page URL, or to the URL in a <base href> tag if one is present, and the result is normalised (lower-case host, default port and #fragment dropped), so that one page is one vertex of the web graph.
• Filtering: mailto:, javascript: and similar schemes are not pages; links marked rel=nofollow, and every link of a page whose <meta name=robots> says nofollow, are not followed, and before fetching a target the crawler checks its site’s robots.txt.

A minimal extractor with Python’s standard-library tokenizer, run on a small sample page (no live fetch):

from html.parser import HTMLParser
from urllib.parse import urljoin, urldefrag
class Links(HTMLParser):
    def __init__(self, url):
        super().__init__(); self.base, self.out = url, []
    def handle_starttag(self, tag, attrs):
        a = dict(attrs)
        if tag == 'base' and a.get('href'): self.base = urljoin(self.base, a['href'])
        if tag in ('a', 'area', 'frame', 'iframe'):
            ref = a.get('href') or a.get('src')
            if ref and 'nofollow' not in (a.get('rel') or ''):
                u = urldefrag(urljoin(self.base, ref.strip()))[0]
                if u.startswith(('http://', 'https://')) and u not in self.out:
                    self.out.append(u)
page = '''<html><head><base href="/docs/"></head><body>
<a href="intro.html#s2">Intro</a> <A HREF='../index.html'>Home</A>
<a href="https://example.org/x?q=1">Ext</a> <a href="mailto:me@x.org">Mail</a>
<a href="javascript:void(0)">JS</a> <a rel="nofollow" href="ad.html">Ad</a>
<iframe src="frame.html"></iframe> <a href="intro.html">Intro again</a></body></html>'''
p = Links('https://www.site.edu/a/page.html'); p.feed(page); print(p.out)