Information Geometry · Study tools

InteractiveAmari, Ch. 1: Manifold, Divergence and Dually Flat Structure

Divergences and Legendre duality

A divergence measures how different two points are without having to be symmetric. Any strictly convex function produces one, and its Legendre transform produces a second coordinate system in which the same divergence reappears with its arguments swapped. This pair of coordinates is the backbone of the whole subject.

The idea

Divergence. A smooth function

D[\xi:\xi']

of two points with

D\ge 0

, equality only when

\xi=\xi'

, and, for nearby points,

D[\xi:\xi+d\xi] = \tfrac12 \sum_{i,j} g_{ij}(\xi)\, d\xi^i d\xi^j + O(|d\xi|^3)

with

g

positive definite. Symmetry and the triangle inequality are not required. The matrix

g

is a Riemannian metric.

Bregman divergence. For a smooth, strictly convex

\psi

on a convex domain,

D_\psi[\theta:\theta'] = \psi(\theta)-\psi(\theta')-\nabla\psi(\theta')\cdot(\theta-\theta'),

the height of the graph of

\psi

at

\theta

above the tangent plane at

\theta'

. Its metric is the Hessian,

g_{ij}=\partial_i\partial_j\psi

.

Legendre transform. Put

\eta=\nabla\psi(\theta)

. Strict convexity makes

\theta\mapsto\eta

one-to-one, and

\psi^*(\eta) = \sup_{\theta}\big(\theta\cdot\eta-\psi(\theta)\big) = \theta\cdot\eta-\psi(\theta), \qquad \theta = \nabla\psi^*(\eta).

\psi^*

is again convex, and its Hessian at

\eta

is the inverse of the Hessian of

\psi

at

\theta

.

Duality. For any

\theta

and

\eta'=\nabla\psi(\theta')

,

D_\psi[\theta:\theta'] = \psi(\theta)+\psi^*(\eta')-\theta\cdot\eta' = D_{\psi^*}[\eta':\eta].

The same number is a gap under

\psi

and a gap under

\psi^*

, with the order of the points reversed.

So a convex function gives two coordinate systems, $\theta$ and $\eta$ , and a divergence that looks like "half a squared distance" up close. The rest of Chapter 1 shows that straight lines in $\theta$ and straight lines in $\eta$ behave like two kinds of geodesic, and a Pythagorean theorem holds when they meet at a right angle. Below, the left plot is $\psi$ and the right plot is $\psi^*$ .

Explore

Convex function ψ(θ)

\psi(\theta)

. Drag the points

\theta

and

\theta'

sideways.

D_\psi[\theta:\theta']

D_\psi[\theta':\theta]

\psi^*(\eta)

, the same two points in the dual coordinate

\eta=\psi'(\theta)

. You can drag here too.

Worked example

Take $\psi(\theta)=e^\theta$ on $\R$ . Find the dual coordinate and the dual function, write $D_\psi$ in terms of $\lambda=e^\theta$ , and compute both orders for $\lambda=1,\ \lambda'=2$ .

What are $\eta$ and $\psi^*(\eta)$ ?

$\eta=\psi'(\theta)=e^\theta=\lambda$ , so $\eta>0$ and $\theta=\log\eta$ . Then

$\psi^*(\eta)=\theta\eta-\psi(\theta)=\eta\log\eta-\eta.$
Write $D_\psi[\theta:\theta']$ in terms of $\lambda,\lambda'$ .

$D_\psi[\theta:\theta'] = e^\theta-e^{\theta'}-e^{\theta'}(\theta-\theta') = \lambda-\lambda'+\lambda'\log\frac{\lambda'}{\lambda}.$

This is the KL divergence from the Poisson distribution with mean $\lambda'$ to the one with mean $\lambda$ . ( $e^\theta$ is the log-partition function of the Poisson family, which is why the explorer labels it that way.)
Check the duality formula using $\psi^*$ .

$D_{\psi^*}[\lambda':\lambda] = \psi^*(\lambda')-\psi^*(\lambda)-(\log\lambda)(\lambda'-\lambda) = \lambda'\log\lambda'-\lambda'-\lambda\log\lambda+\lambda-\lambda'\log\lambda+\lambda\log\lambda,$

which simplifies to $\lambda-\lambda'+\lambda'\log(\lambda'/\lambda)$ , the same expression.
Evaluate both orders for $\lambda=1,\lambda'=2$ .

$D_\psi[\theta:\theta'] = 1-2+2\log2\approx 0.386$ and $D_\psi[\theta':\theta]=2-1+1\cdot\log\tfrac12\approx 0.307$ . The two orders differ. Choose $e^\theta$ in the explorer: it starts at $\theta=0$ , $\theta'=\log 2\approx0.69$ , and shows both gaps.

Check your understanding

For $\psi(\theta)=\tfrac12|\theta|^2$ on $\R^n$ , what is $D_\psi[\theta:\theta']$ ?

$\tfrac12|\theta-\theta'|^2$
$|\theta-\theta'|$ A divergence behaves like half a squared distance near the diagonal, not like a distance. Expand the definition.
$\theta\cdot\theta'$ That can be negative, and a divergence never is.

\tfrac12|\theta|^2-\tfrac12|\theta'|^2-\theta'\cdot(\theta-\theta') = \tfrac12|\theta-\theta'|^2

. Here

\eta=\theta

and

\psi^*=\psi

: the Euclidean case is self-dual, which is why its divergence is symmetric.

Which property does a Bregman divergence fail to have in general?

$D_\psi[\theta:\theta']\ge 0$ That always holds: a convex function lies above each of its tangent planes.
Convexity in the first argument $\theta$ As a function of $\theta$ it is $\psi$ minus an affine function, so it is convex.
Symmetry, $D_\psi[\theta:\theta']=D_\psi[\theta':\theta]$

The gap at

\theta

under the tangent at

\theta'

and the gap at

\theta'

under the tangent at

\theta

are different lengths unless

\psi

curves the same way everywhere between them. In the explorer the red and gold gaps differ for every

\psi

except

\theta^2/2

. It is also usually not convex in its second argument.

$\psi$ has second derivative $\psi''(\theta)=g$ at some $\theta$ . What is $(\psi^*)''$ at $\eta=\psi'(\theta)$ ?

$1/g$
$g$ Test it on $\psi=e^\theta$ : $\psi''=e^\theta=\eta$ , while $\psi^*=\eta\log\eta-\eta$ has $(\psi^*)''=1/\eta$ .
$-g$ $\psi^*$ is a supremum of affine functions of $\eta$ , so it is convex too.

(\psi^*)'(\eta)=\theta(\eta)

is the inverse function of

\eta(\theta)=\psi'(\theta)

, so its derivative is

1/\psi''(\theta)

. In several variables the Hessians are inverse matrices: the metric written in

\eta

-coordinates is the inverse of the metric in

\theta

-coordinates.

Practice

The Bernoulli potentialA Legendre pair and a KL divergence

Let $\psi(\theta)=\log(1+e^\theta)$ . Find $\eta$ and $\psi^*(\eta)$ , and show that $D_\psi[\theta:\theta']$ is the KL divergence $\sum_x p'(x)\log\frac{p'(x)}{p(x)}$ between the Bernoulli distributions with success probabilities $\eta'$ and $\eta$ .

Hint

$\eta$ is the logistic function of $\theta$ , and $\psi^*$ turns out to be minus an entropy. Use $\psi^*$ and the duality formula rather than $\psi$ directly.

Solution

$\eta=e^\theta/(1+e^\theta)$ , so $\theta=\log\frac{\eta}{1-\eta}$ and $\psi(\theta)=-\log(1-\eta)$ . Then $\psi^*(\eta)=\theta\eta-\psi=\eta\log\eta+(1-\eta)\log(1-\eta)$ . By duality

D_\psi[\theta:\theta'] = \psi(\theta)+\psi^*(\eta')-\theta\eta' = -\log(1-\eta)+\eta'\log\eta'+(1-\eta')\log(1-\eta')-\eta'\log\tfrac{\eta}{1-\eta},

which rearranges to $\eta'\log\frac{\eta'}{\eta}+(1-\eta')\log\frac{1-\eta'}{1-\eta}$ .

The three-point identityWhere the Pythagorean theorem comes from

Writing $\eta=\nabla\psi(\theta)$ and similarly for primed points, show that for any three points

D_\psi[\theta:\theta''] = D_\psi[\theta:\theta'] + D_\psi[\theta':\theta''] - (\theta-\theta')\cdot(\eta''-\eta').

When does the last term vanish?

Solution

Expand all three divergences with the definition. The $\psi(\theta')$ terms cancel, and collecting the inner products leaves exactly $-(\theta-\theta')\cdot(\eta''-\eta')$ . It vanishes when the vector from $\theta'$ to $\theta$ , written in $\theta$ -coordinates, is orthogonal to the vector from $\theta'$ to $\theta''$ , written in $\eta$ -coordinates. Because $\theta$ and $\eta$ are dual, this pairing is the Riemannian inner product at $\theta'$ . That is the generalized Pythagorean theorem.

InteractiveAmari, Ch. 2: Exponential Families and Mixture Families of Probability Distributions

Exponential and mixture families

There are two natural ways to interpolate between distributions: mix them, or mix their logarithms. These are the two kinds of straight line in information geometry, the m-geodesic and the e-geodesic. When one meets the other at a right angle, KL divergences add like squared lengths.

The idea

Exponential family.

p(x,\theta)=\exp\big(\theta\cdot t(x)-\psi(\theta)\big)

with respect to a fixed base measure. The log-partition function

\psi

is convex,

\eta=\nabla\psi(\theta)=\mathbb E_\theta[t(x)]

, and

\nabla^2\psi(\theta)=\mathrm{Cov}_\theta[t(x)]

is the Fisher information. So Module 1 applies with

\psi

the log-partition function, and its Bregman divergence is a KL divergence:

D_\psi[\theta:\theta']=\mathrm{KL}[p_{\theta'}:p_\theta]

, where

\mathrm{KL}[p:q]=\sum_x p(x)\log\frac{p(x)}{q(x)}

.

Mixture family.

p(x,\eta)=\sum_i \eta_i\, p_i(x)+\big(1-\sum_i\eta_i\big)p_0(x)

: linear in its parameters. The distributions on

n+1

outcomes, the simplex

S_n

, are both:

\theta^i=\log(p_i/p_0)

makes it exponential and

\eta_i=p_i

makes it a mixture family.

Two geodesics. Between

p

and

q

, the m-geodesic is

p_t=(1-t)p+tq

. The e-geodesic is

\log p_t=(1-t)\log p+t\log q-c(t)

, that is

p_t\propto p^{1-t}q^t

. A set of distributions is m-flat if it contains the m-geodesic between any two of its points, and e-flat for e-geodesics.

Pythagorean theorem. If the m-geodesic from

P

to

Q

and the e-geodesic from

Q

to

R

are orthogonal at

Q

in the Fisher metric, then

\mathrm{KL}[P:R]=\mathrm{KL}[P:Q]+\mathrm{KL}[Q:R].

Projection theorem. Let

M

be m-flat and

R

any distribution. The point

Q\in M

minimizing

\mathrm{KL}[Q:R]

is unique, and it is where an e-geodesic from

R

meets

M

orthogonally (the e-projection). For every

P\in M

,

\mathrm{KL}[P:R]=\mathrm{KL}[P:Q]+\mathrm{KL}[Q:R]

.

The maximum entropy principle is the special case where $R$ is uniform: $\mathrm{KL}[Q:u]=\log(n+1)-H(Q)$ , so maximizing entropy under linear constraints is e-projecting the uniform distribution, and the answer has exponential form.

Explore

The simplex of distributions on three outcomes; each corner is a certain outcome. m-geodesic e-geodesic

Show

Worked example

Find the distribution on $\{0,1,2\}$ with mean $10/7$ that has the largest entropy, and check the Pythagorean relation for another distribution with the same mean.

Before computing: what form must the answer have, and why?

The distributions with mean $10/7$ form an m-flat set $M$ (a linear constraint). Maximizing entropy over $M$ is e-projecting the uniform distribution $u$ onto $M$ , so the answer lies on an e-geodesic from $u$ that meets $M$ orthogonally. The e-geodesics through $u$ are $q_x\propto e^{\lambda v(x)}$ , and orthogonality to every direction inside $M$ (changes $\delta p$ with $\sum\delta p_x=0$ and $\sum x\,\delta p_x=0$ ) forces $v$ to be affine in $x$ . So $q_x\propto e^{\lambda x}$ , that is, $q=(1,r,r^2)/(1+r+r^2)$ with $r=e^\lambda$ .
Solve for $r$ .

The mean is $\dfrac{r+2r^2}{1+r+r^2}=\dfrac{10}{7}$ , so $7r+14r^2=10+10r+10r^2$ , that is $4r^2-3r-10=0$ , with positive root $r=2$ . Hence $q=(1,2,4)/7$ .
Compute $H(q)$ .

$\log q_x = x\log2-\log7$ , so $H(q)=\log7-\tfrac{10}{7}\log2\approx 0.956$ nats (against $\log 3\approx1.099$ for the uniform distribution).
Take $p=(2/7,\,0,\,5/7)$ , which also has mean $10/7$ . Show $\mathrm{KL}[p:u]=\mathrm{KL}[p:q]+\mathrm{KL}[q:u]$ .

Since $\log q_x$ is affine in $x$ , $\sum_x p_x\log q_x$ depends on $p$ only through its mean, so it equals $\sum_x q_x\log q_x=-H(q)$ . Therefore $\mathrm{KL}[p:q]=-H(p)+H(q)$ , and

$\mathrm{KL}[p:q]+\mathrm{KL}[q:u] = \big(H(q)-H(p)\big)+\big(\log3-H(q)\big) = \log 3-H(p)=\mathrm{KL}[p:u].$

Numerically $H(p)\approx0.598$ , so $\mathrm{KL}[p:q]\approx0.357$ . The step "depends only through the mean" is exactly the orthogonality in the theorem.

Check your understanding

$p=(0.8,0.1,0.1)$ and $q=(0.1,0.1,0.8)$ . Which is the midpoint of the e-geodesic between them?

$(0.45,\,0.10,\,0.45)$ That is the m-midpoint, the average of the probabilities.
$(0.425,\,0.150,\,0.425)$
$(0.5,\,0,\,0.5)$ An e-geodesic between distributions with full support never reaches the boundary.

The e-midpoint is the normalized geometric mean:

\sqrt{pq}=(0.283,0.1,0.283)

, which sums to

0.666

; dividing gives

(0.425,0.150,0.425)

. The geometric mean punishes disagreement, so the outcome on which

p

and

q

agree gains weight.

$Q$ is the e-projection of $R$ onto an m-flat set $M$ , and $P$ is another point of $M$ . Which is guaranteed?

$\mathrm{KL}[P:R]\ge\mathrm{KL}[Q:R]$ , with equality only if $P=Q$
$\mathrm{KL}[R:P]\ge\mathrm{KL}[R:Q]$ The order matters. The e-projection minimizes $\mathrm{KL}[\,\cdot:R]$ , with the point of $M$ first. Minimizing $\mathrm{KL}[R:\cdot\,]$ is a different projection.
$\mathrm{KL}[P:Q]=\mathrm{KL}[Q:P]$ KL is not symmetric, even between points of a flat set.

Pythagoras gives

\mathrm{KL}[P:R]=\mathrm{KL}[P:Q]+\mathrm{KL}[Q:R]

and

\mathrm{KL}[P:Q]>0

for

P\ne Q

. That is why the e-projection is the unique minimizer. Switch the explorer to "Projection" and slide

P

along the line.

Which of these sets of distributions on three outcomes is both e-flat and m-flat?

$\{p : p_1 = 2p_3\}$
$\{p : p_1-p_3 = 0.2\}$ It is linear in $p$ , so m-flat. But the e-midpoint of $(0.5,0.2,0.3)$ and $(0.3,0.6,0.1)$ has $p_1-p_3\approx0.236$ .
$\{p : p_1p_3 = p_2^2\}$ In logarithms it reads $\log p_1+\log p_3-2\log p_2=0$ , which is linear in $\theta$ , so it is e-flat. Is the m-midpoint of $(\tfrac13,\tfrac13,\tfrac13)$ and $(\tfrac47,\tfrac27,\tfrac17)$ in it?

p_1=2p_3

is linear in

p

(m-flat), and it is equivalent to

\log p_1-\log p_3=\log2

, which is linear in

\theta^i=\log(p_i/p_0)

(e-flat). In the simplex picture it is a straight line through a corner.

Practice

The Bernoulli familyNatural and expectation coordinates

Write the Bernoulli distribution $p(x)=\eta^x(1-\eta)^{1-x}$ , $x\in\{0,1\}$ , as an exponential family. Find $\theta$ , $\psi(\theta)$ , and the Fisher information in the coordinate $\theta$ and in the coordinate $\eta$ .

Solution

$\log p = x\log\frac{\eta}{1-\eta}+\log(1-\eta)$ , so $t(x)=x$ , $\theta=\log\frac{\eta}{1-\eta}$ and $\psi(\theta)=\log(1+e^\theta)$ . In $\theta$ : $g=\psi''(\theta)=\eta(1-\eta)$ , the variance of $x$ . In $\eta$ : $g=1/(\eta(1-\eta))$ , the inverse, as Module 1 predicts.

Mutual information as a projectionProjecting onto independent distributions

Let $p(x,y)$ be a joint distribution on a finite set, with marginals $p_X,p_Y$ . Show that among all product distributions $q_1(x)q_2(y)$ (an e-flat set), $\mathrm{KL}[p:q_1q_2]$ is minimized by $q_1=p_X,\ q_2=p_Y$ , and the minimum is the mutual information $I(X;Y)$ .

Hint

Split $\log q_1(x)q_2(y)$ into its two terms and sum out the other variable in each.

Solution

$\mathrm{KL}[p:q_1q_2]=-H(X,Y)-\sum_x p_X(x)\log q_1(x)-\sum_y p_Y(y)\log q_2(y)$ $= -H(X,Y)+H(X)+\mathrm{KL}[p_X:q_1]+H(Y)+\mathrm{KL}[p_Y:q_2]$ . The two KL terms vanish exactly when $q_1=p_X$ , $q_2=p_Y$ , leaving $H(X)+H(Y)-H(X,Y)=I(X;Y)$ . Here the second argument is projected onto an e-flat set (an m-projection), the dual of the explorer's situation.

InteractiveAmari, Ch. 3: Invariant Geometry of Manifold of Probability Distributions

Invariance and the Fisher metric

How should we measure the difference between two distributions? Ask only that merging outcomes can never make them easier to tell apart, and that a merge which loses no information changes nothing. Those two demands single out the f-divergences, and they force the metric to be the Fisher metric, up to a constant factor.

The idea

f-divergence. For a convex

f

on

(0,\infty)

with

f(1)=0

and

f''(1)=1

,

D_f[p:q]=\sum_i p_i\, f\!\left(\frac{q_i}{p_i}\right).

f(u)=-\log u

gives

\mathrm{KL}[p:q]

;

f(u)=u\log u

gives

\mathrm{KL}[q:p]

;

f(u)=2(\sqrt u-1)^2

gives

2\sum_i(\sqrt{p_i}-\sqrt{q_i})^2

.

Coarse graining. Merge the outcomes into blocks

A_1,\dots,A_m

:

\bar p_j=\sum_{i\in A_j}p_i

. The merge is sufficient for

p,q

if

q_i/p_i

is constant on each block: then nothing that distinguishes

p

from

q

is lost.

Information monotonicity.

D_f[\bar p:\bar q]\le D_f[p:q]

for every coarse graining, with equality when it is sufficient. (Jensen's inequality inside each block.)

Uniqueness. On three or more outcomes, a divergence of the form

\sum_i F(p_i,q_i)

that is monotone in this sense is an f-divergence (up to a constant factor). A Riemannian metric on

S_n

that is invariant under sufficient embeddings is a constant multiple of the Fisher metric

g_{ij}(\xi)=\mathbb E_\xi\big[\partial_i\log p(x,\xi)\ \partial_j\log p(x,\xi)\big], \qquad ds^2=\sum_i\frac{dp_i^2}{p_i}\ \text{ on } S_n.

Every f-divergence with $f''(1)=1$ has this same metric, $D_f[p:p+dp]\approx\tfrac12\sum_i dp_i^2/p_i$ . The substitution $\xi_i=2\sqrt{p_i}$ turns $ds^2$ into $\sum d\xi_i^2$ and $S_n$ into part of a sphere of radius 2, so the Fisher distance is $d(p,q)=2\arccos\sum_i\sqrt{p_iq_i}$ . Below, compare the level sets of $\mathrm{KL}[P:\cdot\,]$ , $\mathrm{KL}[\,\cdot:P]$ and $\tfrac12 d^2$ : they agree near $P$ and part company further away.

Explore

Drag P and Q. Curves are the levels 0.01, 0.03, 0.08, 0.18, 0.35, 0.7, 1.4 of the chosen divergence from P.

Level sets of

Worked example

Let $p=(0.5,0.3,0.2)$ and $q=(0.2,0.3,0.5)$ . Merge outcomes 2 and 3 and compare KL divergences before and after; then change $q$ so that the merge loses nothing.

Compute $\mathrm{KL}[p:q]$ .

$0.5\log2.5+0.3\log1+0.2\log0.4\approx0.458+0-0.183=0.275$ .
Compute $\mathrm{KL}[\bar p:\bar q]$ for $\bar p=(0.5,0.5)$ , $\bar q=(0.2,0.8)$ .

$0.5\log2.5+0.5\log0.625\approx0.458-0.235=0.223$ . Smaller, as monotonicity demands.
Where exactly is the loss? Use the ratios $q_i/p_i$ .

On the merged block the ratios are $q_2/p_2=1$ and $q_3/p_3=2.5$ . They differ, so the merged block mixes evidence that pointed in different directions. The loss is $\bar p_2$ times a KL divergence of the within-block conditionals: $0.5\,\mathrm{KL}[(0.6,0.4):(0.375,0.625)]\approx0.052$ , and indeed $0.275-0.223=0.052$ .
Keep $q_1=0.2$ and choose $q_2,q_3$ so the merge is sufficient. What is $\mathrm{KL}[p:q]$ now?

We need $q_2/q_3=p_2/p_3=3/2$ with $q_2+q_3=0.8$ , so $q=(0.2,0.48,0.32)$ . Then $\mathrm{KL}[p:q]=0.5\log2.5+0.3\log0.625+0.2\log0.625\approx0.223$ , equal to the merged value. The "Make the merge sufficient" button does this to Q.

Check your understanding

Which of these is not an f-divergence on the simplex?

$\mathrm{KL}[q:p]$ Take $f(u)=u\log u$ : $\sum p_i\frac{q_i}{p_i}\log\frac{q_i}{p_i}=\sum q_i\log\frac{q_i}{p_i}$ .
$2\sum_i(\sqrt{p_i}-\sqrt{q_i})^2$ Take $f(u)=2(\sqrt u-1)^2$ .
$\sum_i(p_i-q_i)^2$

Split one outcome into

k

equally likely copies. That refinement is sufficient, so an invariant divergence must not change. But the squared Euclidean distance replaces

(a-b)^2

by

k(a/k-b/k)^2=(a-b)^2/k

. An f-divergence depends only on the ratios

q_i/p_i

, which splitting leaves alone.

$p=(0.4,0.4,0.2)$ . After merging outcomes 2 and 3, for which $q$ is $\mathrm{KL}[p:q]$ unchanged?

$q=(0.2,0.4,0.4)$ The ratios on the merged block are $1$ and $2$ : the merge throws away the difference between them.
$q=(0.1,0.6,0.3)$
$q=(0.4,0.2,0.4)$ Ratios $0.5$ and $2$ on the merged block. Not sufficient.

Here

q_2/p_2=q_3/p_3=1.5

, so within the merged block

q

is just a rescaled copy of

p

: the merged outcome carries the same evidence as the two separate ones. Equality in monotonicity is exactly this condition.

All f-divergences with $f''(1)=1$ agree with $\tfrac12ds^2$ to second order. So when $Q$ is close to $P$ , how big is $\mathrm{KL}[P:Q]-\mathrm{KL}[Q:P]$ ?

Third order in the distance between them
Second order, comparable to the divergences themselvesBoth orders have the same second-order term $\tfrac12 ds^2$ , so it cancels in the difference.
Exactly zeroThe cubic terms differ. Drag Q close to P and watch the two KL readings: they agree to more and more digits, but not exactly.

The asymmetry of KL is a third-order effect: near the diagonal every invariant divergence looks like half the squared Fisher distance. (The third-order term is where the dual connections of Chapter 6 come from.)

Practice

Every f-divergence gives the Fisher metricA Taylor expansion

Show that $D_f[p:p+dp]=\tfrac12 f''(1)\sum_i dp_i^2/p_i+O(|dp|^3)$ for $p$ and $p+dp$ in the simplex.

Hint

Write $q_i/p_i=1+u_i$ with $u_i=dp_i/p_i$ , expand $f$ around 1, and use $\sum_i dp_i=0$ .

Solution

$f(1+u)=f'(1)u+\tfrac12f''(1)u^2+O(u^3)$ . Summing $p_if(1+u_i)$ : the linear term is $f'(1)\sum_i dp_i=0$ , and the quadratic term is $\tfrac12f''(1)\sum_i p_i u_i^2=\tfrac12 f''(1)\sum_i dp_i^2/p_i$ .

Fisher distance through the sphereA closed form

Show that $\xi_i=2\sqrt{p_i}$ sends the simplex into the sphere $|\xi|=2$ and turns $\sum dp_i^2/p_i$ into $\sum d\xi_i^2$ . Deduce $d(p,q)=2\arccos\sum_i\sqrt{p_iq_i}$ , and evaluate it for $p=(\tfrac12,\tfrac12)$ , $q=(0.9,0.1)$ .

Solution

$\sum\xi_i^2=4\sum p_i=4$ and $d\xi_i=dp_i/\sqrt{p_i}$ , so $d\xi_i^2=dp_i^2/p_i$ . Distances on a sphere of radius 2 are $2\times$ (angle), and $\cos(\text{angle})=\xi\cdot\xi'/4=\sum\sqrt{p_iq_i}$ . (The image is the positive part of the sphere, which contains the great-circle arc between any two of its points, so this is the distance within the simplex too.) Here $\sum\sqrt{p_iq_i}=\sqrt{0.45}+\sqrt{0.05}\approx0.894$ , so $d\approx2\times0.464=0.927$ .

InteractiveAmari, Ch. 4: α-Geometry, Tsallis q-Entropy and Positive-Definite Matrices

α-divergences and α-geodesics

The m-geodesic averages probabilities; the e-geodesic averages their logarithms. Average their powers instead and you get a whole family of curves, means and divergences indexed by one number $\alpha$ , with $\alpha$ and $-\alpha$ always dual to each other.

The idea

α-representation.

\ell_\alpha(u)=\dfrac{2}{1-\alpha}\,u^{(1-\alpha)/2}

for

\alpha\ne1

, and

\ell_1(u)=\log u

. So

\alpha=-1

is

u

itself (m),

\alpha=1

is

\log u

(e), and

\alpha=0

is

2\sqrt u

.

α-divergence. For positive vectors

m,n

and

\alpha\ne\pm1

,

D_\alpha[m:n]=\frac{4}{1-\alpha^2}\sum_i\Big(\tfrac{1-\alpha}{2}m_i+\tfrac{1+\alpha}{2}n_i-m_i^{(1-\alpha)/2}\,n_i^{(1+\alpha)/2}\Big).

As

\alpha\to-1

it tends to

\sum_i\big(m_i\log\frac{m_i}{n_i}-m_i+n_i\big)

, the KL divergence, and as

\alpha\to1

to the same with

m,n

swapped. Always

D_\alpha[m:n]=D_{-\alpha}[n:m]

. It is the f-divergence with

f(u)=\frac{4}{1-\alpha^2}\big(1-u^{(1+\alpha)/2}\big)+\frac{2}{1-\alpha}(u-1)

.

α-geodesic. On positive measures it is the straight line in

\ell_\alpha

:

m_i(t)^{(1-\alpha)/2}=(1-t)\,m_i^{(1-\alpha)/2}+t\,n_i^{(1-\alpha)/2}

. On the simplex, rescale

p(t)

to sum to 1. For

\alpha=0

this is the great-circle arc of Module 3.

Pythagoras, with a correction. Positive measures are α-flat for every

\alpha

, and the Pythagorean theorem holds exactly there. The simplex is flat only for

\alpha=\pm1

. For other

\alpha

, when the α-geodesic

pq

and the

(-\alpha)

-geodesic

qr

are orthogonal,

D_\alpha[p:r]=D_\alpha[p:q]+D_\alpha[q:r]-\tfrac{1-\alpha^2}{4}\,D_\alpha[p:q]\,D_\alpha[q:r],

the spherical correction for a space of constant curvature.

The same exponents give the α-mean of positive numbers, $m_\alpha(a,b)=\ell_\alpha^{-1}\big(\tfrac12\ell_\alpha(a)+\tfrac12\ell_\alpha(b)\big)$ : arithmetic for $\alpha=-1$ , geometric for $\alpha=1$ , harmonic for $\alpha=3$ . Tsallis's q-logarithm $\ln_q u=(u^{1-q}-1)/(1-q)$ is the α-representation shifted by a constant, with $q=(1+\alpha)/2$ . The chapter also builds invariant divergences on positive-definite matrices (Practice below).

Explore

α-geodesic for reference, the m- and e-geodesics (α = −1, 1). Drag P and Q.

α

Jump to:

Worked example

Compute the α-mean of $a=1$ and $b=4$ for several $\alpha$ , and find its limits as $\alpha\to\pm\infty$ .

Write $m_\alpha(a,b)$ explicitly for $\alpha\ne1$ .

With $r=(1-\alpha)/2$ , the constant $2/(1-\alpha)$ cancels and $m_\alpha(a,b)=\Big(\dfrac{a^r+b^r}{2}\Big)^{1/r}$ : a power mean with exponent $r$ .
Evaluate it for $\alpha=-1,0,1,3$ .

$\alpha=-1$ ( $r=1$ ): $2.5$ . $\alpha=0$ ( $r=\tfrac12$ ): $\big(\tfrac{1+2}{2}\big)^2=2.25$ . $\alpha=1$ (the limit $r\to0$ ): $\sqrt{1\cdot4}=2$ . $\alpha=3$ ( $r=-1$ ): $\big(\tfrac{1+1/4}{2}\big)^{-1}=1.6$ .
What happens as $\alpha\to-\infty$ and $\alpha\to+\infty$ ?

$r\to+\infty$ : the power mean tends to $\max(a,b)=4$ . $r\to-\infty$ : it tends to $\min(a,b)=1$ . So $\alpha$ slides the mean continuously and monotonically from optimistic (max) to pessimistic (min).
Connect this to the explorer.

Before rescaling, the α-geodesic midpoint of $P$ and $Q$ is their α-mean taken coordinate by coordinate. For large $\alpha$ each coordinate is pulled toward the smaller of $p_i,q_i$ , so after rescaling the midpoint gives more weight to outcomes on which $P$ and $Q$ agree. For very negative $\alpha$ the larger value wins and the midpoint spreads toward both of their favourite outcomes. Put $P$ and $Q$ near different corners and move α from −3 to 3.

Check your understanding

For which $\alpha$ is $D_\alpha$ symmetric in its two arguments?

$\alpha=0$
$\alpha=1$ $D_1[p:q]=\mathrm{KL}[q:p]$ , which is not symmetric.
Every $\alpha$ Use the duality $D_\alpha[p:q]=D_{-\alpha}[q:p]$ : symmetry needs $D_\alpha=D_{-\alpha}$ .

Duality gives

D_\alpha[p:q]=D_{-\alpha}[q:p]

, so

D_0

is self-dual and symmetric. It equals

2\sum(\sqrt{p_i}-\sqrt{q_i})^2

on the simplex. Set

\alpha=0

above: the two readings coincide.

For which $\alpha$ does the plain Pythagorean relation $D_\alpha[p:r]=D_\alpha[p:q]+D_\alpha[q:r]$ hold exactly on the simplex (under the orthogonality condition)?

Only $\alpha=\pm1$
Only $\alpha=0$ At $\alpha=0$ the correction factor $(1-\alpha^2)/4$ is largest: the simplex is a piece of a sphere for that geometry.
Every $\alpha$ That is true for positive measures, which are α-flat for every $\alpha$ . The simplex is a curved slice of them.

The correction term carries the factor

1-\alpha^2

, which vanishes only at

\alpha=\pm1

: the e- and m-geometries, in which the simplex is dually flat.

As $\alpha$ increases from $-1$ to $1$ , the α-mean of 1 and 4…

decreases from 2.5 to 2
increases from 2 to 2.5 $\alpha=-1$ is the arithmetic mean, the larger one.
stays at 2.5Only when $a=b$ do all the means agree.

The exponent

r=(1-\alpha)/2

falls from 1 to 0, and power means increase with their exponent. On the way,

\alpha=0

gives

2.25

.

Practice

α = 0 is a function of the Fisher distanceHellinger and the sphere

Show that on the simplex $D_0[p:q]=8\sin^2(d/4)$ , where $d=2\arccos\sum_i\sqrt{p_iq_i}$ is the Fisher distance.

Solution

For probability vectors the first two terms of $D_0$ sum to 1, so $D_0=4\big(1-\sum\sqrt{p_iq_i}\big)=4(1-\cos(d/2))=8\sin^2(d/4)$ . So $D_0$ is an increasing function of the Fisher distance, and $D_0\approx d^2/2$ for small $d$ .

A divergence for covariance matricesInvariance under change of basis

On positive-definite $n\times n$ matrices take $\psi(P)=-\log\det P$ . Show that its Bregman divergence is $D[P:Q]=\mathrm{tr}(PQ^{-1})-\log\det(PQ^{-1})-n$ , that it is unchanged by $P\mapsto A^{\mathsf T}PA$ , $Q\mapsto A^{\mathsf T}QA$ for invertible $A$ , and what it becomes for diagonal matrices.

Hint

The gradient of $\log\det$ at $Q$ is $Q^{-1}$ , in the sense $d\log\det Q=\mathrm{tr}(Q^{-1}dQ)$ .

Solution

$D=-\log\det P+\log\det Q+\mathrm{tr}\big(Q^{-1}(P-Q)\big)$ , which is the stated formula. Under the change of basis $PQ^{-1}\mapsto A^{\mathsf T}PQ^{-1}A^{-\mathsf T}$ , a similar matrix, so trace and determinant are unchanged. For diagonal matrices it is $\sum_i\big(p_i/q_i-\log(p_i/q_i)-1\big)$ , a sum of Itakura–Saito divergences. It is also twice the KL divergence between the zero-mean Gaussians with covariances $P$ and $Q$ .

InteractiveAmari, Ch. 5: Elements of Differential Geometry

The Fisher metric on Gaussians

The normal distributions $N(\mu,\sigma^2)$ form a surface, the upper half-plane. With the Fisher metric it is a hyperbolic plane: distances blow up as $\sigma\to0$ , and the shortest path between two narrow distributions with different means goes through wider ones.

The idea

Riemannian metric. An inner product

g_{ij}(\xi)

on each tangent space. A curve

\xi(t)

has length

\int\sqrt{g_{ij}\dot\xi^i\dot\xi^j}\,dt

, and the distance between points is the infimum of lengths. The length of a curve does not depend on the coordinates used to describe it.

Connection and geodesics. An affine connection says how to compare tangent vectors at nearby points (covariant derivative, parallel transport). A geodesic is a curve whose velocity is parallel along itself. The Levi-Civita connection is the unique torsion-free connection that preserves

g

; its geodesics are locally the shortest curves. Curvature measures how parallel transport around a small loop rotates vectors.

The normal family. For

N(\mu,\sigma^2)

the Fisher metric is

ds^2=\frac{d\mu^2+2\,d\sigma^2}{\sigma^2}.

With

u=\mu/\sqrt2

this is

2\,(du^2+d\sigma^2)/\sigma^2

: twice the Poincaré half-plane metric, a surface of constant curvature

-\tfrac12

.

So geodesics are the hyperbolic ones: vertical lines, and half-circles centred on $\sigma=0$ in the $(u,\sigma)$ plane, which are half-ellipses $(\mu-c)^2/2+\sigma^2=R^2$ in $(\mu,\sigma)$ . The distance is

d=\sqrt2\,\operatorname{arccosh}\!\Big(1+\frac{(\mu_1-\mu_2)^2/2+(\sigma_1-\sigma_2)^2}{2\sigma_1\sigma_2}\Big).

Explore

The

(\mu,\sigma)

half-plane. geodesic, with dots at equal Fisher spacing straight segment. Small ellipses are Fisher circles of equal radius. Drag A and B.

The densities of A (blue), B (gold) and the evenly spaced distributions along the geodesic.

Worked example

Derive the Fisher metric of $N(\mu,\sigma^2)$ , then find the Fisher distance from $N(-1,1)$ to $N(1,1)$ and compare it with the path that keeps $\sigma=1$ .

Compute the two score functions.

$\log p=-\log\sigma-\frac{(x-\mu)^2}{2\sigma^2}+\text{const}$ . With $z=(x-\mu)/\sigma$ , $\partial_\mu\log p=z/\sigma$ and $\partial_\sigma\log p=(z^2-1)/\sigma$ .
Take expectations of their products.

$z$ is standard normal: $\mathbb E z^2=1$ , $\mathbb E z(z^2-1)=0$ , $\mathbb E(z^2-1)^2=\mathbb Ez^4-2\mathbb Ez^2+1=2$ . So $g_{\mu\mu}=1/\sigma^2$ , $g_{\mu\sigma}=0$ , $g_{\sigma\sigma}=2/\sigma^2$ .
Find the length of the horizontal path $\sigma=1$ from $\mu=-1$ to $\mu=1$ , and the geodesic distance.

Along $\sigma=1$ , $ds=d\mu$ , so the length is 2. The formula gives $d=\sqrt2\operatorname{arccosh}(1+\tfrac{2}{2})=\sqrt2\operatorname{arccosh}2\approx1.862$ . The horizontal path is not a geodesic.
Which distributions does the geodesic pass through?

In $(u,\sigma)$ the endpoints are $(\mp1/\sqrt2,1)$ , so the circle is centred at $u=0$ with radius $\sqrt{1/2+1}\approx1.225$ . The geodesic rises to $N(0,1.225^2)$ halfway: it is cheaper to widen the distribution, shift it, and narrow it again than to slide a narrow one sideways. Put A and B at these points above.

Check your understanding

In the Fisher metric, which is farther from $N(0,1)$ : $N(0,0.5^2)$ or $N(0,2^2)$ ?

$N(0,2^2)$ , because $\sigma$ changes by 1 rather than 0.5The metric is $\sqrt2\,d\sigma/\sigma$ along vertical lines: what counts is the ratio of the $\sigma$ 's, not their difference.
$N(0,0.5^2)$ , because narrow distributions are far apartNarrowness matters, but here both changes are by the same factor of 2.
They are equally far

Along

\mu=0

the distance is

\int\sqrt2\,d\sigma/\sigma=\sqrt2\,|\log(\sigma_2/\sigma_1)|

, which is

\sqrt2\log2\approx0.980

for both. Doubling

x

is an invertible change of units, so it preserves the Fisher geometry, and it maps

N(0,0.5^2)\mapsto N(0,1)

and

N(0,1)\mapsto N(0,2^2)

.

The geodesic between $N(-a,\sigma^2)$ and $N(a,\sigma^2)$ passes through…

wider distributions than the endpoints
only distributions with the same $\sigma$ Drag A and B to equal heights: the straight segment is longer than the geodesic.
narrower distributions than the endpointsDistances grow like $1/\sigma$ as $\sigma$ shrinks, so dipping down is expensive.

Horizontal movement costs

d\mu/\sigma

, so it is cheaper at larger

\sigma

. The geodesic climbs to where moving sideways is cheap, like a great-circle flight route bending toward the pole.

We switch coordinates from $(\mu,\sigma)$ to $(\mu,v)$ with $v=\sigma^2$ . What happens to the Fisher distance between two given normal distributions?

Nothing; only the formula for $g_{ij}$ changes
It changes, since $g_{vv}\ne g_{\sigma\sigma}$ The components change, but so do the coordinate differentials, and $ds^2$ is the same quadratic form.
It is multiplied by 2Lengths of curves are defined without reference to coordinates.

d\sigma=dv/(2\sqrt v)

, so

2\,d\sigma^2/\sigma^2=dv^2/(2v^2)

:

g_{vv}=1/(2v^2)

, a different number, but the same

ds^2

. The Fisher metric is a tensor on the space of distributions, not a property of a parametrization.

Practice

The curvature is −1/2Scaling a metric

Assuming the Poincaré metric $(du^2+d\sigma^2)/\sigma^2$ has Gaussian curvature $-1$ , show that the Fisher metric of the normal family has curvature $-\tfrac12$ .

Solution

Multiplying a metric by a constant $c$ multiplies all lengths by $\sqrt c$ and areas by $c$ . Gaussian curvature is angle excess per unit area (Gauss–Bonnet), and angles do not change, so $K$ is divided by $c$ . With $c=2$ , $K=-\tfrac12$ .

A one-dimensional familyDistance in the exponential family of rates

For the exponential distributions $p(x;\lambda)=\lambda e^{-\lambda x}$ , $x>0$ , compute the Fisher information $g(\lambda)$ and the Fisher distance between rates $\lambda_1$ and $\lambda_2$ .

Solution

$\partial_\lambda\log p=1/\lambda-x$ , with variance $\mathrm{Var}(x)=1/\lambda^2$ . So $g=1/\lambda^2$ and $d=\int_{\lambda_1}^{\lambda_2}d\lambda/\lambda=|\log(\lambda_2/\lambda_1)|$ . Again only ratios matter, because changing the unit of $x$ is invertible and multiplies every rate by the same factor.

InteractiveAmari, Ch. 6: Dual Affine Connections and Dually Flat Manifold

Dual connections

A statistical manifold carries more than a metric: it carries two notions of "straight", the e- and m-connections, and they are dual with respect to the Fisher metric. On the Gaussian half-plane you can see all three geodesics at once and test the Pythagorean theorem with a protractor.

The idea

Dual connections.

\nabla

and

\nabla^*

are dual with respect to

g

if

X\langle Y,Z\rangle=\langle\nabla_XY,Z\rangle+\langle Y,\nabla^*_XZ\rangle

. Equivalently, transporting one vector by

\nabla

and another by

\nabla^*

keeps their inner product constant. A connection is self-dual exactly when it preserves

g

, like the Levi-Civita connection.

The α-connections. With the cubic (Amari–Chentsov) tensor

T_{ijk}=\mathbb E[\partial_i\ell\,\partial_j\ell\,\partial_k\ell]

,

\ell=\log p

, define

\Gamma^{(\alpha)}_{ijk}=\Gamma^{(0)}_{ijk}-\tfrac\alpha2T_{ijk}

, where

\Gamma^{(0)}

is Levi-Civita. Then

\nabla^{(\alpha)}

and

\nabla^{(-\alpha)}

are dual;

\alpha=1

is the e-connection,

\alpha=-1

the m-connection, and their average is Levi-Civita. Any divergence produces such a pair: its second derivatives give

g

, its third derivatives give the two connections.

Dually flat manifolds. If

\nabla

is flat with affine coordinates

\theta

, then

\nabla^*

is flat too, with affine coordinates

\eta=\nabla\psi(\theta)

for a convex potential

\psi

. The canonical divergence is the Bregman divergence of Module 1, and the Pythagorean theorem holds. Exponential families are dually flat:

\theta

natural parameters,

\eta

expectation parameters.

For $N(\mu,\sigma^2)$ : $\theta=\big(\mu/\sigma^2,\,-1/(2\sigma^2)\big)$ , $\eta=(\mu,\ \mu^2+\sigma^2)$ . The e-geodesic is a straight line in $\theta$ , the m-geodesic a straight line in $\eta$ , and the Levi-Civita geodesic is the hyperbolic one from Module 5. Mixing the two systems also works: coordinates such as $(\eta_1,\theta_2)$ cut a dually flat manifold into two families of flat slices that meet orthogonally.

Explore

e-geodesic m-geodesic Levi-Civita geodesic. Axes

\mu

(across) and

\sigma

(up).

Show

Worked example

Write the normal family in natural parameters, find its potential $\psi$ , check that $\nabla\psi$ gives $\eta$ , and find the dual potential $\psi^*$ .

Put $\log p$ in the form $\theta\cdot t(x)-\psi(\theta)$ .

$\log p=\frac{\mu}{\sigma^2}x-\frac{1}{2\sigma^2}x^2-\frac{\mu^2}{2\sigma^2}-\tfrac12\log(2\pi\sigma^2)$ . So $t=(x,x^2)$ , $\theta_1=\mu/\sigma^2$ , $\theta_2=-1/(2\sigma^2)$ , and in terms of $\theta$

$\psi(\theta)=-\frac{\theta_1^2}{4\theta_2}+\tfrac12\log\Big(-\frac{\pi}{\theta_2}\Big).$
Differentiate $\psi$ .

$\partial_1\psi=-\theta_1/(2\theta_2)=\mu$ and $\partial_2\psi=\theta_1^2/(4\theta_2^2)-1/(2\theta_2)=\mu^2+\sigma^2$ . These are $\mathbb E[x]$ and $\mathbb E[x^2]$ , as they must be.
Compute $\psi^*(\eta)=\theta\cdot\eta-\psi(\theta)$ .

$\theta\cdot\eta=\frac{\mu^2}{\sigma^2}-\frac{\mu^2+\sigma^2}{2\sigma^2}=\frac{\mu^2}{2\sigma^2}-\frac12$ , and $\psi=\frac{\mu^2}{2\sigma^2}+\frac12\log(2\pi\sigma^2)$ . So $\psi^*=-\tfrac12\log(2\pi e\sigma^2)$ , minus the entropy of $N(\mu,\sigma^2)$ , with $\sigma^2=\eta_2-\eta_1^2$ .
Why are the e- and m-geodesics between $N(-1,1)$ and $N(1,1)$ so different?

Both endpoints have $\theta_2=-\tfrac12$ , so the e-geodesic keeps $\theta_2$ and hence $\sigma=1$ fixed: it is the horizontal segment. Both have $\eta_2=2$ , so the m-geodesic keeps $\mu^2+\sigma^2=2$ : it is an arc of the circle of radius $\sqrt2$ , reaching $\sigma=\sqrt2$ at $\mu=0$ . The Levi-Civita geodesic lies between, peaking at $\sigma\approx1.225$ .

Check your understanding

What is the average $\tfrac12(\nabla^{(1)}+\nabla^{(-1)})$ of the e- and m-connections?

The Levi-Civita connection of the Fisher metric
The flat connection of the $\theta$ coordinatesThat is the e-connection itself, one of the two being averaged.
It is not a connectionAn average of two connections (weights summing to 1) is again a connection.

From

\Gamma^{(\pm1)}=\Gamma^{(0)}\mp\tfrac12T

, the cubic terms cancel in the average. In the explorer the blue curve always runs between the gold and the green ones.

$A=N(-1,1)$ , $B=N(1,1)$ . Which geodesic from $A$ to $B$ keeps $\sigma$ constant?

The e-geodesic
The m-geodesicThe m-geodesic averages $\eta_2=\mu^2+\sigma^2$ . In the middle $\mu^2$ is smaller, so $\sigma^2$ must be larger.
The Levi-Civita geodesicThat one bends upward (Module 5).

\theta_2=-1/(2\sigma^2)

is the same at both ends, and the e-geodesic is linear in

\theta

. Place A and B at the same height to see it.

The m-geodesic from $A$ to $B$ and the e-geodesic from $B$ to $C$ are orthogonal at $B$ . Which identity follows?

$\mathrm{KL}[A:C]=\mathrm{KL}[A:B]+\mathrm{KL}[B:C]$
$\mathrm{KL}[C:A]=\mathrm{KL}[C:B]+\mathrm{KL}[B:A]$ That version needs the roles swapped: an m-geodesic from $C$ to $B$ orthogonal to an e-geodesic from $B$ to $A$ .
The Fisher distances satisfy $d_{AC}^2=d_{AB}^2+d_{BC}^2$ The half-plane is curved ( $K=-\tfrac12$ ), so the Riemannian Pythagoras is only approximate, for small triangles.

The m-side carries the first argument of KL, the e-side the second. In the explorer's triangle mode the defect is exactly

-(\theta_C-\theta_B)\cdot(\eta_A-\eta_B)

, the pairing that measures the angle at

B

.

Practice

Dual transport keeps inner productsBiorthogonal bases

In a dually flat manifold, let $X$ have constant components in the $\theta$ basis $\partial/\partial\theta^i$ (e-parallel) and $Y$ constant components in the $\eta$ basis $\partial/\partial\eta_j$ (m-parallel). Show that $\langle X,Y\rangle$ is constant.

Hint

Compute $\langle\partial/\partial\theta^i,\partial/\partial\eta_j\rangle$ using $\partial\eta_k/\partial\theta^i=g_{ik}$ .

Solution

$\partial/\partial\theta^i=\sum_k g_{ik}\,\partial/\partial\eta_k$ , and the $\eta$ -basis has inner products $g^{kj}$ (the inverse matrix). So $\langle\partial_{\theta^i},\partial_{\eta_j}\rangle=\sum_k g_{ik}g^{kj}=\delta_i^j$ . With $X=a^i\partial_{\theta^i}$ , $Y=b_j\partial_{\eta_j}$ , $\langle X,Y\rangle=\sum_ia^ib_i$ , which is constant when $a$ and $b$ are.

Mixed coordinates for two binary variablesMarginals and interaction

For $p(x,y)$ on $\{0,1\}^2$ , use $\eta_1=P(x=1)$ , $\eta_2=P(y=1)$ and $\theta_{12}=\log\frac{p_{11}p_{00}}{p_{10}p_{01}}$ . Show that independence is $\theta_{12}=0$ , and explain why projecting $p$ onto the independent distributions (minimizing $\mathrm{KL}[p:q]$ ) keeps $\eta_1,\eta_2$ fixed.

Solution

$p$ is a product exactly when its odds ratio is 1. Write $\log p$ in the exponential form $\theta_1x+\theta_2y+\theta_{12}xy-\psi$ ; then $\theta_{12}$ is the log odds ratio, and $\{\theta_{12}=0\}$ is e-flat. The sets of fixed marginals $\{\eta_1,\eta_2\text{ fixed}\}$ are m-flat (linear constraints), and in mixed coordinates the two families meet orthogonally. The m-projection runs along an m-geodesic orthogonal to the e-flat target, so it stays in the slice with $p$ 's marginals and lands at $p_X\,p_Y$ , matching the mutual-information problem of Module 2.

SimulationAmari, Ch. 7: Asymptotic Theory of Statistical Inference

Efficient estimation

Data from an exponential family land on a point of its expectation space. A model is a curve in that space, and an estimator is a rule for sliding the observed point back onto the curve. Whether the estimator wastes information is decided by one angle: the one at which it slides.

The idea

Curved exponential family. A model

M=\{p(x,u)\}

sitting inside an exponential family

S

as a curve (or submanifold)

\eta(u)

. After

N

observations the average of the sufficient statistics,

\hat\eta=\frac1N\sum_k t(x_k)

, is a point of

S

, usually not on

M

. It contains everything the data say about

u

.

Estimators and ancillary families. An estimator is a map

\hat u=f(\hat\eta)

. Its level sets

A(u)=\{\eta: f(\eta)=u\}

form the ancillary family. The estimator is consistent when each

A(u)

passes through

\eta(u)

.

First-order efficiency. Asymptotically

N\,\mathrm{Var}(\hat u)\ge 1/g(u)

, with

g

the Fisher information (Cramér–Rao), and a consistent estimator attains

1/g(u)

exactly when

A(u)

is orthogonal to

M

at

\eta(u)

in the Fisher metric.

The MLE is an m-projection.

\frac1N\sum_k\log p(x_k,u)=\theta(u)\cdot\hat\eta-\psi(\theta(u))

, so maximizing likelihood is minimizing

\mathrm{KL}[p_{\hat\eta}:p_u]

over

M

. The minimizer is where an m-geodesic from

\hat\eta

meets

M

orthogonally, so the MLE's ancillary family is m-flat and orthogonal: the MLE is efficient.

Higher-order terms (bias, and the $1/N^2$ part of the variance) depend on how curved $M$ and $A(u)$ are, measured by the e- and m-embedding curvatures; the MLE corrected for bias is second-order efficient. Below the model is $N(u,u^2)$ , whose standard deviation equals its mean, inside the normal family with $\eta=(\mathbb E x,\mathbb E x^2)$ . True value $u=1$ .

Explore

Axes

\eta_1=\bar x

(across) and

\eta_2=\overline{x^2}

(up). model

\eta(u)=(u,2u^2)

the estimator's ancillary lines simulated

\hat\eta

's. The grey region

\eta_2<\eta_1^2

is impossible.

Estimator

Sample size N

Worked example

For $N(u,u^2)$ , $u>0$ : find the MLE, the Fisher information, and the asymptotic efficiency of $\bar x$ .

Write the average log-likelihood in terms of $m_1=\bar x$ and $m_2=\overline{x^2}$ .

$\tfrac1N\ell(u) = -\log u-\frac{m_2-2um_1+u^2}{2u^2}+\text{const} = -\log u-\frac{m_2}{2u^2}+\frac{m_1}{u}+\text{const}.$

Only $(m_1,m_2)=\hat\eta$ enters: it is sufficient.
Solve the likelihood equation.

$-\frac1u+\frac{m_2}{u^3}-\frac{m_1}{u^2}=0$ gives $u^2+m_1u-m_2=0$ , so $\hat u=\tfrac12\big(-m_1+\sqrt{m_1^2+4m_2}\big)$ . For fixed $u$ the solutions form the line $\eta_2=u\,\eta_1+u^2$ through $\eta(u)=(u,2u^2)$ : the MLE's ancillary family is a family of straight lines, as the theory says (m-flat).
Compute the Fisher information $g(u)$ .

Use the normal-family metric $(d\mu^2+2d\sigma^2)/\sigma^2$ from Module 5 along $\mu=\sigma=u$ : $g(u)=(1+2)/u^2=3/u^2$ . So the best possible asymptotic variance is $u^2/(3N)$ .
What is the efficiency of $\bar x$ ?

$\mathrm{Var}(\bar x)=u^2/N$ , so its efficiency is $\frac{u^2/3N}{u^2/N}=\tfrac13$ . It ignores the information about $u$ carried by the spread of the data. Its ancillary lines are vertical, at a Fisher angle of about $35^\circ$ to the model instead of $90^\circ$ , and $\sin^2 35.26^\circ=\tfrac13$ .

Check your understanding

An estimator's ancillary line $A(u)$ is the set of observed points it maps to $u$ . What makes the estimator consistent?

Each $A(u)$ passes through the model point $\eta(u)$
Each $A(u)$ is orthogonal to the modelThat is efficiency. A consistent estimator can cross at any angle.
The lines $A(u)$ are straightStraightness (m-flatness) matters for higher-order properties, not for consistency.

As

N\to\infty

,

\hat\eta\to\eta(u_0)

by the law of large numbers, so

\hat u\to f(\eta(u_0))

, which is

u_0

exactly when

\eta(u_0)\in A(u_0)

. All three estimators in the explorer are consistent; only one is efficient.

In a full exponential family (not curved), what is the MLE of the expectation parameter $\eta$ ?

$\hat\eta$ itself, which attains the Cramér–Rao bound for every $N$
It has to be found numericallyMaximizing $\theta\cdot\hat\eta-\psi(\theta)$ gives $\nabla\psi(\theta)=\hat\eta$ directly.
$\hat\eta$ , but it is biased $\mathbb E\hat\eta=\mathbb E\,t(x)=\eta$ .

\hat\eta

is unbiased with covariance

\mathrm{Cov}(t)/N

, and the Fisher information in

\eta

-coordinates is

\mathrm{Cov}(t)^{-1}

. So it is exactly efficient. Curving the model is what makes estimation subtle.

$\sqrt{\overline{x^2}/2}$ has efficiency $8/9$ here. How many observations does it need to match the MLE with $N=800$ ?

900
711That inverts the ratio. A less efficient estimator needs more data.
800Both are consistent, but their variances differ by a constant factor at every large $N$ .

Asymptotic variances are

u^2/(3N)

for the MLE and

3u^2/(8N)

for this estimator. Equal variance needs

N'=800\times\frac{3/8}{1/3}=900

. Efficiency is a ratio of sample sizes.

Practice

Efficiency from the angleA Fisher inner product

At $u=1$ , $(x,x^2)$ has covariance $\Sigma=\begin{pmatrix}1&2\\2&6\end{pmatrix}$ , and the Fisher metric in $\eta$ -coordinates is $\Sigma^{-1}$ . Compute the angle between the model tangent $(1,4)$ and the vertical direction $(0,1)$ , and check that $\sin^2$ of it is $\tfrac13$ .

Solution

$\Sigma^{-1}=\tfrac12\begin{pmatrix}6&-2\\-2&1\end{pmatrix}$ . Then $\langle t,t\rangle=3$ , $\langle a,a\rangle=\tfrac12$ , $\langle t,a\rangle=\tfrac12(-2+4)=1$ . So $\cos^2=\frac{1}{3\cdot\frac12}=\tfrac23$ and $\sin^2=\tfrac13$ . The same computation with $(1,1)$ , the MLE direction, gives $\langle t,a\rangle=\tfrac12(6-2-8+4)=0$ .

The log-odds of a coinEfficiency in another parametrization

With $N$ tosses of a coin with $P(\text{heads})=p$ , the MLE is $\hat p$ , the frequency of heads. Find the asymptotic variance of the log-odds $\hat\theta=\log\frac{\hat p}{1-\hat p}$ , and compare it with the Cramér–Rao bound for $\theta$ .

Solution

$d\theta/dp=1/(p(1-p))$ , so by the delta method $\mathrm{Var}(\hat\theta)\approx\frac{p(1-p)/N}{(p(1-p))^2}=\frac{1}{Np(1-p)}$ . The Fisher information for $\theta$ is $\psi''(\theta)=p(1-p)$ , so this is the bound: the MLE stays efficient under reparametrization. (Only asymptotically: $\hat\theta$ is infinite with positive probability for every finite $N$ .)

Animated lessonAmari, Ch. 8: Estimation in the Presence of Hidden Variables

The EM algorithm

When some variables are hidden, maximum likelihood gets hard, but two easy projections remain. EM alternates them: an e-projection fills in the hidden part as best it can, an m-projection refits the model, and the likelihood of the observed data can never go down.

The idea

Two manifolds. Let

x

be observed and

h

hidden. The model manifold is

M=\{p(x,h;\xi)\}

. The data manifold is

D=\{q(x,h): q(x)=\hat p(x)\}

, all joint distributions whose

x

-marginal is the empirical distribution.

D

is m-flat: it is cut out by linear constraints.

ML is the distance between them. By the chain rule for KL,

\mathrm{KL}[q:p_\xi]=\mathrm{KL}[\hat p(x):p(x;\xi)]+\sum_x\hat p(x)\,\mathrm{KL}\big[q(h|x):p(h|x;\xi)\big].

Minimizing over

q\in D

kills the second term, and the first is

-\tfrac1N\log L(\xi)

plus a constant. So

\min_{q\in D,\ \xi}\mathrm{KL}[q:p_\xi]

is maximum likelihood.

E-step (e-projection onto D). For the current

\xi

, the minimizer over

q\in D

is

q(x,h)=\hat p(x)\,p(h|x;\xi)

: fill in the hidden variable by its posterior.

M-step (m-projection onto M). Minimize

\mathrm{KL}[q:p_\xi]

over

\xi

, that is, maximize

\mathbb E_q[\log p(x,h;\xi)]

. For an exponential family this matches expected sufficient statistics.

Monotonicity. Neither step increases

\mathrm{KL}[q:p_\xi]

, and right after each E-step it equals

-\tfrac1N\log L(\xi)+\text{const}

. So the log-likelihood never decreases. Typically it converges to a stationary point, which need not be the global maximum.

The animation fits a mixture $\tfrac12N(a,1)+\tfrac12N(b,1)$ with unknown means; the hidden variable is which component produced each point.

Watch

Data

The

(a,b)

plane: contours of the log-likelihood and the EM path. Drag the start point ξ₀.

component

a

component

b

; dots below are the data, coloured by their posterior.

The data, and a first guess $\xi_0=(a,b)$ for the two means.
E-step: give each point its posterior probability of coming from each component (e-projection onto $D$ ).
M-step: move each mean to the posterior-weighted average of the data (m-projection onto $M$ ).
Repeat. The path climbs the likelihood surface; every step goes uphill.
It settles at a stationary point of the likelihood. Try other starts, including one on the diagonal $a=b$ .

Worked example

Data $x=(-2,\ 0,\ 3)$ , model $\tfrac12N(a,1)+\tfrac12N(b,1)$ , start at $a=-1,\ b=1$ . Do one EM iteration by hand.

Find the posterior $r_k=P(h=b\mid x_k)$ as a function of $x_k$ .

$r=\dfrac{e^{-(x-b)^2/2}}{e^{-(x-a)^2/2}+e^{-(x-b)^2/2}}=\dfrac{1}{1+e^{((x-b)^2-(x-a)^2)/2}}$ . With $a=-1,b=1$ , $(x-b)^2-(x-a)^2=-4x$ , so $r=1/(1+e^{-2x})$ .
Evaluate the posteriors.

$r=(0.018,\ 0.5,\ 0.998)$ for $x=-2,0,3$ . The middle point is split evenly; the others are almost certain.
M-step: compute the new means.

$b'=\frac{\sum r_kx_k}{\sum r_k}=\frac{-0.036+0+2.993}{1.516}\approx1.951,\qquad a'=\frac{\sum(1-r_k)x_k}{\sum(1-r_k)}=\frac{-1.964+0.007}{1.484}\approx-1.318.$
Check that the log-likelihood went up.

$\sum_k\log\big(\tfrac12\phi(x_k-a)+\tfrac12\phi(x_k-b)\big)$ is about $-7.12$ at $(-1,1)$ and $-6.18$ at $(-1.318,1.951)$ . The guarantee comes from the geometry: the E-step made the divergence equal to minus the log-likelihood, and the M-step could only lower it.

Check your understanding

The data manifold $D=\{q(x,h): q(x)=\hat p(x)\}$ is…

m-flat
e-flatFixing a marginal is a linear condition on $q$ , not on $\log q$ . The normalized geometric mean of two members can have a different marginal.
neitherMix two members: $(1-t)q+tq'$ still has marginal $\hat p$ .

The marginal is linear in

q

, so mixtures of members are members. That is why the E-step is an e-projection onto

D

: the projection theorem of Module 2 needs the target to be m-flat when the projection is along an e-geodesic.

You start EM exactly on the diagonal, $a=b$ . What happens?

It stays on the diagonal and converges to $a=b=\bar x$ , a saddle point
It leaves the diagonal at once, because the likelihood is higher elsewhereEM is deterministic, and with $a=b$ every posterior is exactly $\tfrac12$ . Nothing breaks the symmetry.
It fails, because the M-step divides by zeroEach component gets total weight $N/2$ , so the weighted averages are defined.

With

a=b

, both components get posterior

\tfrac12

on every point, so both new means are

\bar x

. The point

(\bar x,\bar x)

is a stationary point of the likelihood. Along

a=\bar x+t,\ b=\bar x-t

the log-likelihood has second derivative

\sum_k\big((x_k-\bar x)^2-1\big)

, so it is a saddle, not a maximum, exactly when the sample variance exceeds 1, as it does for both data sets here. Drag the start onto the diagonal (it snaps) and play.

Which quantity is guaranteed never to decrease from one EM iteration to the next?

The log-likelihood of the observed data
The log-likelihood computed with the true hidden labelsEM never sees the true labels; nothing controls that quantity.
The closeness of $\xi$ to the true parameterEM can climb to a local maximum or saddle far from the truth.

After an E-step,

\mathrm{KL}[q:p_\xi]=-\tfrac1N\log L(\xi)+c

. The M-step lowers KL to some value

\ge-\tfrac1N\log L(\xi')+c

(by the chain rule), so

\log L(\xi')\ge\log L(\xi)

. The "Δ log-lik" reading in the animation is never negative.

Practice

The E-step is a projectionThe chain rule for KL

Show that $\mathrm{KL}[q:p]=\mathrm{KL}[q_X:p_X]+\sum_xq(x)\,\mathrm{KL}[q(\cdot|x):p(\cdot|x)]$ and conclude that over $q\in D$ the minimum of $\mathrm{KL}[q:p]$ is at $q(h|x)=p(h|x)$ .

Solution

Write $q(x,h)=q(x)q(h|x)$ and the same for $p$ ; the logarithm of the ratio splits into a marginal part and a conditional part, and summing gives the identity. On $D$ , $q_X=\hat p$ is fixed, so only the second term varies, and it is zero exactly when every conditional matches.

Unknown mixing weightExtending the M-step

Let the model be $(1-w)N(a,1)+wN(b,1)$ with $w$ unknown too. Given posteriors $r_k$ , find the M-step update for $w$ .

Solution

The expected complete-data log-likelihood contains $\sum_k\big(r_k\log w+(1-r_k)\log(1-w)\big)$ . Setting its derivative to zero gives $w'=\frac1N\sum_kr_k$ , the average posterior. The means update as before.

SimulationAmari, Ch. 9: Neyman-Scott Problem: Estimating Function and Semiparametric Statistical Model

Nuisance parameters

If every new observation brings its own new unknown, maximum likelihood can converge confidently to the wrong answer. An estimating function that is blind to the nuisance directions repairs it.

The idea

The Neyman–Scott problem. We observe pairs

(x_i,y_i)

,

i=1,\dots,n

, independent, with

x_i,y_i\sim N(\mu_i,\sigma^2)

. The parameter of interest is

\sigma^2

; the means

\mu_1,\mu_2,\dots

are nuisance parameters, one new one for each pair.

The MLE is inconsistent. Maximizing over all

\mu_i

gives

\hat\mu_i=(x_i+y_i)/2

and

\hat\sigma^2_{\mathrm{ML}}=\frac{1}{4n}\sum_i(x_i-y_i)^2\ \longrightarrow\ \frac{\sigma^2}{2}.

More data never fixes it: each

\hat\mu_i

is always based on two observations.

Estimating function. A function

y(x,\sigma^2)

with

\mathbb E_{\sigma^2,\mu}[y]=0

for every value of the nuisance. The estimate solves

\sum_iy(x_i,\hat\sigma^2)=0

, and it is consistent under mild conditions. Here

y=\tfrac12(x-y)^2-\sigma^2

works, giving

\tilde\sigma^2=\frac1{2n}\sum_i(x_i-y_i)^2

.

Geometry. Over each value of

\sigma^2

the model has a whole fibre of nuisance directions. A good estimating function is orthogonal (in the Fisher metric) to all of them. The best one is the efficient score: the score for

\sigma^2

minus its projection onto the nuisance directions. Here

x_i+y_i

carries all the information about

\mu_i

, and removing from the score its conditional expectation given

x_i+y_i

leaves a multiple of

\tfrac12(x_i-y_i)^2-\sigma^2

.

Explore

Nuisance means μᵢ

Estimates of

\sigma^2

(true value 1) as pairs arrive, log scale in

n

. MLE estimating function wrong model: one common mean

Each pair as a point

(x_i,y_i)

. The nuisance moves points along the diagonal;

\sigma

is the spread across it.

Worked example

Derive the MLE of $\sigma^2$ in the Neyman–Scott model and find its limit.

Maximize the log-likelihood over $\mu_i$ with $\sigma^2$ fixed.

$\ell=-n\log(2\pi\sigma^2)-\frac{1}{2\sigma^2}\sum_i\big((x_i-\mu_i)^2+(y_i-\mu_i)^2\big)$ . Each $\mu_i$ appears in one term, minimized at $\hat\mu_i=(x_i+y_i)/2$ , where $(x_i-\hat\mu_i)^2+(y_i-\hat\mu_i)^2=\tfrac12(x_i-y_i)^2$ .
Now maximize over $\sigma^2$ .

The profile log-likelihood is $-n\log\sigma^2-\frac{1}{4\sigma^2}\sum(x_i-y_i)^2+\text{const}$ . Setting the derivative to zero, $\hat\sigma^2=\frac{1}{4n}\sum(x_i-y_i)^2$ .
Find the limit as $n\to\infty$ .

$x_i-y_i\sim N(0,2\sigma^2)$ whatever $\mu_i$ is, so $\mathbb E(x_i-y_i)^2=2\sigma^2$ , and by the law of large numbers $\hat\sigma^2\to2\sigma^2/4=\sigma^2/2$ .
Why did the usual consistency theorem fail?

It assumes a fixed, finite number of parameters. Here the number of parameters, $n+1$ , grows with the data, so the $\hat\mu_i$ never settle down, and their errors bias $\hat\sigma^2$ by the same factor at every $n$ . The estimating function $\tfrac12(x-y)^2-\sigma^2$ avoids $\mu_i$ entirely.

Check your understanding

In the explorer, switch the nuisance means to "all equal to 0". The MLE still estimates a separate $\mu_i$ for each pair. Where does it converge now?

Still to $\sigma^2/2$
To $\sigma^2$ , since there is no variation in the means to confuse itThe estimator only sees $x_i-y_i$ , whose distribution does not depend on the $\mu_i$ at all.
To 0 $\mathbb E(x_i-y_i)^2=2\sigma^2>0$ .

The bias comes from fitting

n

free means, not from their values. What changes with the spread is the "common mean" estimator, which is wrong when the means vary and fine when they are all equal.

Suppose each $\mu_i$ is observed three times instead of twice. Where does the MLE of $\sigma^2$ converge?

$\tfrac23\sigma^2$
$\tfrac13\sigma^2$ The MLE divides the within-group sum of squares by $3n$ ; that sum has expectation $(3-1)n\sigma^2$ .
$\sigma^2$ More observations per nuisance parameter reduce the bias but do not remove it.

With

k

observations per group,

\hat\sigma^2=\frac{1}{kn}\sum_i\sum_j(x_{ij}-\bar x_i)^2\to\frac{k-1}{k}\sigma^2

. One degree of freedom per group is spent on estimating

\mu_i

.

What makes $y=\tfrac12(x-y)^2-\sigma^2$ a valid estimating function?

Its expectation is zero for every value of the nuisance $\mu_i$
It is the score function for $\sigma^2$ The score is $-\tfrac1{\sigma^2}+\frac{(x-\mu)^2+(y-\mu)^2}{2\sigma^4}$ , which involves the unknown $\mu$ .
Its expectation is zero when $\mu_i=0$ That would not be enough: we need unbiasedness whatever the $\mu_i$ are, since we never learn them.

x-y\sim N(0,2\sigma^2)

for every

\mu

, so

\mathbb E[\tfrac12(x-y)^2]=\sigma^2

regardless. Averaging over pairs then gives a consistent estimate, by the law of large numbers.

Practice

Four pairsThe two estimates side by side

The differences $x_i-y_i$ for four pairs are $0.4,\ -1.0,\ 1.6,\ 0.2$ . Compute $\hat\sigma^2_{\mathrm{ML}}$ and $\tilde\sigma^2$ .

Solution

$\sum d_i^2=0.16+1+2.56+0.04=3.76$ . So $\hat\sigma^2_{\mathrm{ML}}=3.76/16=0.235$ and $\tilde\sigma^2=3.76/8=0.47$ , always in the ratio 1 : 2.

Conditioning away the nuisanceA likelihood that works

Show that $d_i=x_i-y_i$ and $s_i=x_i+y_i$ are independent, that the distribution of $d_i$ does not involve $\mu_i$ , and that maximizing the likelihood of the $d_i$ alone gives $\tilde\sigma^2$ .

Solution

$(d_i,s_i)$ is a linear image of a Gaussian vector with $\mathrm{Cov}(d_i,s_i)=\sigma^2-\sigma^2=0$ , so they are jointly Gaussian and uncorrelated, hence independent. $d_i\sim N(0,2\sigma^2)$ . The MLE of the variance $2\sigma^2$ from the $d_i$ is $\frac1n\sum d_i^2$ , so $\sigma^2$ is estimated by $\frac1{2n}\sum d_i^2=\tilde\sigma^2$ .

Step-throughAmari, Ch. 11: Machine Learning

Clustering with divergences

k-means works for any Bregman divergence, not only squared distance. Surprisingly, the centre of a cluster is still its plain average and the boundaries between clusters are still straight lines. What the divergence changes is where those lines go.

The idea

Bregman centroid. For points

x_1,\dots,x_m

and any Bregman divergence,

\sum_kD_\psi[x_k:c]=\sum_kD_\psi[x_k:\bar x]+m\,D_\psi[\bar x:c],

so

\arg\min_c\sum_kD_\psi[x_k:c]=\bar x

, the arithmetic mean, whatever

\psi

is.

The other side.

\arg\min_c\sum_kD_\psi[c:x_k]

is the mean taken in the dual coordinates,

\nabla\psi^*\big(\tfrac1m\sum_k\nabla\psi(x_k)\big)

. For the generalized KL divergence that is the geometric mean, for Itakura–Saito the harmonic mean.

Straight boundaries.

D_\psi[x:c_1]-D_\psi[x:c_2]

is an affine function of

x

: the

\psi(x)

terms cancel. So the cells

\{x: D_\psi[x:c_j]\text{ smallest}\}

are convex polygons.

Bregman k-means. Alternate: assign each point to the centre with the smallest

D_\psi[x:c_j]

; move each centre to the mean of its points. Each step lowers

\sum_kD_\psi[x_k:c_{j(k)}]

or leaves it alone, so the algorithm stops after finitely many changes.

The chapter applies the same geometry to support vector machines, belief propagation, boosting and Bayesian learning; this module takes clustering. Below are positive 2-D data with multiplicative noise, so a divergence that measures relative error behaves very differently from squared distance. Per coordinate: squared distance $\tfrac12(x-c)^2$ ; generalized KL $x\log\frac xc-x+c$ ; Itakura–Saito $\frac xc-\log\frac xc-1$ .

Explore

Shaded cells: where each centre is closest in the chosen divergence. Drag the centres (crosses), or step the algorithm.

Divergence

Worked example

Two numbers, 1 and 4. Find both centroids for squared distance, generalized KL and Itakura–Saito.

Prove the centroid identity in the box above.

Expand: $\sum_kD[x_k:c]-\sum_kD[x_k:\bar x]=m\psi(\bar x)-m\psi(c)-\nabla\psi(c)\cdot\sum_k(x_k-c)+\nabla\psi(\bar x)\cdot\sum_k(x_k-\bar x)$ . The last sum is 0 and $\sum_k(x_k-c)=m(\bar x-c)$ , leaving $m\big(\psi(\bar x)-\psi(c)-\nabla\psi(c)\cdot(\bar x-c)\big)=mD[\bar x:c]$ .
What are the "data first" centroids of 1 and 4?

All three are $2.5$ , the arithmetic mean.
What are the "centre first" centroids?

Squared distance: $\nabla\psi(x)=x$ , so again $2.5$ . Generalized KL, $\psi=x\log x-x$ : $\nabla\psi=\log x$ , so the centroid is $\exp\big(\tfrac{\log1+\log4}{2}\big)=2$ . Itakura–Saito, $\psi=-\log x$ : $\nabla\psi=-1/x$ , so the centroid is $1/\big(\tfrac{1+1/4}{2}\big)=1.6$ .
Which side does the explorer use, and what difference does the divergence make there?

Data first, so every centre update is a plain average. The divergence enters only through the assignment: Itakura–Saito depends only on the ratio $x/c$ , so it hands a big outlying point to a big centre that is far away in absolute terms, while squared distance hands it to the nearest one.

Check your understanding

Under the generalized KL divergence $D[x:c]=\sum_j\big(x_j\log\frac{x_j}{c_j}-x_j+c_j\big)$ , which $c$ minimizes $\sum_kD[x_k:c]$ ?

The arithmetic mean of the $x_k$
The geometric meanThat minimizes $\sum_kD[c:x_k]$ , with the centre first.
The coordinatewise medianThe median minimizes absolute deviations, which are not a Bregman divergence.

The centroid identity holds for every Bregman divergence when the data are in the first slot. Directly:

\partial_{c_j}\sum_k(-x_{kj}\log c_j+c_j)=-\sum_kx_{kj}/c_j+m=0

, so

c_j=\bar x_j

.

With Itakura–Saito and data in the first slot, the boundary between two clusters is…

a straight line
a curve that bends around the smaller centreWrite out $D[x:c_1]-D[x:c_2]$ : the $\log x_j$ terms cancel and what is left is affine in $x$ .
a circleEven for squared distance the boundary between two centres is a line, the perpendicular bisector.

D[x:c_1]-D[x:c_2]=\sum_jx_j\big(\tfrac1{c_{1j}}-\tfrac1{c_{2j}}\big)+\sum_j\log\frac{c_{1j}}{c_{2j}}

, affine in

x

. The line is just not the perpendicular bisector. The shaded cells in the explorer are always polygons.

Why can each step of Bregman k-means only lower the objective $\sum_kD[x_k:c_{j(k)}]$ ?

Assignment picks the smallest divergence for each point, and the update picks the minimizing centre (the mean) for each cluster
Because Bregman divergences are symmetricThey are not (Module 1), and no symmetry is used.
Because the objective is convex in the centres and the assignments togetherIt is not jointly convex: k-means can stop at different local solutions from different starts. Try "New random start".

Each step minimizes the objective exactly over one block of variables with the other fixed: a coordinate-descent argument, the same alternating-minimization pattern as EM in Module 8.

Practice

The other-sided centroidUsing the dual divergence

Prove that $\arg\min_c\sum_kD_\psi[c:x_k]=\nabla\psi^*\big(\tfrac1m\sum_k\nabla\psi(x_k)\big)$ .

Hint

Module 1: $D_\psi[c:x]=D_{\psi^*}[\nabla\psi(x):\nabla\psi(c)]$ .

Solution

In dual coordinates the problem is $\min\sum_kD_{\psi^*}[y_k:\gamma]$ with $y_k=\nabla\psi(x_k)$ , $\gamma=\nabla\psi(c)$ . By the centroid identity for $\psi^*$ the minimizer is $\gamma=\bar y$ , so $c=\nabla\psi^*(\bar y)$ .

The KL boundaryAn explicit line

For the generalized KL divergence in two dimensions, find the boundary between centres $c_1=(1,1)$ and $c_2=(2,1)$ .

Solution

$D[x:c_1]-D[x:c_2]=\sum_jx_j\log\frac{c_{2j}}{c_{1j}}+\sum_j(c_{1j}-c_{2j})=x_1\log2-1$ . The boundary is the vertical line $x_1=1/\log2\approx1.443$ , not the bisector $x_1=1.5$ : it sits closer to the smaller centre, as relative error suggests.

SimulationAmari, Ch. 12: Natural Gradient Learning and Its Dynamics in Singular Regions

Natural gradient

Ordinary gradient descent depends on how you happen to parametrize a model. The natural gradient asks a different question: which small change of the distribution, measured by the Fisher metric, lowers the loss fastest? Its path does not depend on the parametrization.

The idea

Steepest descent in a metric. Minimizing

L(\xi+d\xi)

subject to

\tfrac12\,d\xi^{\mathsf T}G(\xi)\,d\xi=\varepsilon^2

gives

d\xi\propto-G^{-1}\nabla L

. The natural gradient is

\tilde\nabla L=G^{-1}\nabla L

, and learning is

\xi_{t+1}=\xi_t-\eta\,G^{-1}(\xi_t)\nabla L(\xi_t)

.

Invariance. Under a change of coordinates both

\nabla L

and

G

transform, and

G^{-1}\nabla L

transforms as a tangent vector. So the natural gradient flow traces the same curve of distributions in any coordinates; discrete steps agree to first order in

\eta

.

Efficiency. In online learning from one example at a time, natural gradient with step

\eta_t=1/t

is asymptotically Fisher-efficient: it achieves the Cramér–Rao bound of Module 7. Plain gradient descent generally does not.

Singular regions and plateaus. In a multilayer perceptron, different parameters can give the same function: when two hidden units coincide, or a unit's output weight is zero. On these sets the Fisher metric degenerates. Plain gradient learning gets trapped near them for long periods (plateaus), while the natural gradient, which rescales by

G^{-1}

, escapes them much faster.

Below, a normal model $N(\mu,\sigma^2)$ is fitted to the target $N(1,0.4^2)$ by minimizing $L=\mathrm{KL}[\text{target}:\text{model}]=\log\frac{\sigma}{0.4}+\frac{0.4^2+(\mu-1)^2}{2\sigma^2}-\frac12$ . With $G=\mathrm{diag}(1/\sigma^2,2/\sigma^2)$ the natural gradient step is $\Delta\mu=-\eta(\mu-1)$ , $\Delta\sigma=-\tfrac\eta2\big(\sigma-\frac{0.4^2+(\mu-1)^2}{\sigma}\big)$ .

Explore

Plain gradient descent in coordinates

Step size η

natural gradient plain gradient. Contours of

L

. Drag the start point; the target is the ring.

\log_{10}L

against iteration (first 150).

Worked example

Fit a Bernoulli model in its natural parameter $\theta$ (the log-odds) to data with frequency $\hat p=0.8$ , by minimizing the average negative log-likelihood $L(\theta)=\psi(\theta)-\theta\hat p$ , $\psi=\log(1+e^\theta)$ . Start at $\theta=0$ and use step size 1.

Compute $\nabla L$ and the Fisher information $G$ .

$L'(\theta)=\psi'(\theta)-\hat p=p-\hat p$ with $p=1/(1+e^{-\theta})$ , and $G=\psi''(\theta)=p(1-p)$ .
Take one plain and one natural gradient step from $\theta=0$ .

At $\theta=0$ : $p=0.5$ , $L'=-0.3$ , $G=0.25$ . Plain: $\theta=0.3$ , so $p=0.574$ . Natural: $\theta=0+0.3/0.25=1.2$ , so $p=0.769$ .
Take a second natural gradient step. How close is it to the answer $\theta^*=\log4\approx1.386$ ?

At $\theta=1.2$ : $p=0.7685$ , $L'=-0.0315$ , $G=0.1779$ , so $\theta=1.2+0.177=1.377$ , with $p=0.7985$ . The error went from 0.186 to 0.009: quadratic convergence.
Explain the speed.

For an exponential family in natural parameters, $L''(\theta)=\psi''(\theta)=G(\theta)$ : the Fisher information is the Hessian of the loss. So the natural gradient step with $\eta=1$ is exactly Newton's method here.

Check your understanding

From the same start, with a very small step, plain gradient descent on $L(\mu,\sigma)$ and on $L(\mu,\log\sigma)$ trace…

different paths of distributions, while the natural gradient traces the same path in both
the same path, since it is the same functionThe gradient's components change: $\partial L/\partial\log\sigma=\sigma\,\partial L/\partial\sigma$ , and a step in $\log\sigma$ moves $\sigma$ by $\sigma$ times as much again. The $\sigma$ -velocity is multiplied by $\sigma^2$ .
the same path, but at different speedsOnly the $\sigma$ component is rescaled, so the direction changes too.

Plain gradient descent secretly uses the Euclidean metric of whatever coordinates it is given. The natural gradient uses the Fisher metric, which belongs to the distributions. Toggle the coordinates in the explorer: only the red path moves.

Why does plain gradient descent in $(\mu,\sigma)$ need tiny steps when $\sigma$ is small?

The loss curvature in $\mu$ is about $1/\sigma^2$ , so a stable step must be of order $\sigma^2$
The loss is not differentiable at small $\sigma$ It is smooth for every $\sigma>0$ .
The gradient vanishes thereThe opposite: it is huge, like $1/\sigma^2$ and $1/\sigma^3$ .

Gradient descent on a quadratic with curvature

k

is stable only for

\eta<2/k

. The natural gradient multiplies by

G^{-1}

, which contains the factor

\sigma^2

, so its step adapts automatically.

What causes the long plateaus in gradient training of multilayer perceptrons?

Singular regions, where hidden units coincide or vanish and the Fisher metric degenerates
Local minimaA plateau is not a minimum: learning does eventually escape it, after a long slow phase.
A step size that is too smallLarger steps do not remove the plateau; the slow direction comes from the geometry near the singular set.

Near a singular region, many parameters represent nearly the same function, so the gradient in the escaping direction is tiny. In the metric that measures change of function, that direction is not short at all, and the natural gradient moves along it at normal speed.

Practice

A Gaussian meanThe natural gradient forgets σ

For $N(\mu,\sigma^2)$ with $\sigma$ known and loss $L(\mu)=\frac{1}{2\sigma^2}\cdot\frac1N\sum_k(x_k-\mu)^2$ , compute the plain and natural gradient steps in $\mu$ .

Solution

$L'(\mu)=(\mu-\bar x)/\sigma^2$ and $G=1/\sigma^2$ . Plain: $\Delta\mu=-\eta(\mu-\bar x)/\sigma^2$ , whose behaviour depends on $\sigma$ . Natural: $\Delta\mu=-\eta(\mu-\bar x)$ , the same for every $\sigma$ ; with $\eta=1$ it jumps to $\bar x$ at once.

Why it is invariantThe transformation rule

Let $\xi=f(\zeta)$ with Jacobian $J=\partial\xi/\partial\zeta$ . Show that $J\,G_\zeta^{-1}\nabla_\zeta L=G_\xi^{-1}\nabla_\xi L$ : the two natural gradients are the same tangent vector.

Solution

$\nabla_\zeta L=J^{\mathsf T}\nabla_\xi L$ and $G_\zeta=J^{\mathsf T}G_\xi J$ . So $G_\zeta^{-1}\nabla_\zeta L=J^{-1}G_\xi^{-1}J^{-\mathsf T}J^{\mathsf T}\nabla_\xi L=J^{-1}G_\xi^{-1}\nabla_\xi L$ ; multiply by $J$ . The plain gradient $\nabla_\zeta L$ instead becomes $JJ^{\mathsf T}\nabla_\xi L$ in $\xi$ -coordinates, which depends on $f$ .

SimulationAmari, Ch. 13: Signal Processing and Optimization

Independent component analysis

Two independent signals are mixed by an unknown matrix, and we see only the mixtures. ICA searches the group of invertible matrices for one that makes the outputs independent again. The natural gradient on that group makes the learning dynamics independent of how badly the signals were mixed.

The idea

The model.

x=As

, where

s

has independent, non-Gaussian components and

A

is invertible and unknown. We look for

W

such that

y=Wx

has independent components. At best this recovers

s

up to the order and scale of its components.

The loss. Choose model densities

q_i

for the sources. The negative log-likelihood per sample is

L(W)=-\log|\det W|-\sum_i\mathbb E\log q_i(y_i)

, with gradient

\nabla L=-W^{-\mathsf T}+\mathbb E[\varphi(y)x^{\mathsf T}]

, where

\varphi_i=-(\log q_i)'

.

Natural gradient on the matrix group. Use a metric on invertible matrices that is invariant under

W\mapsto WB

(it measures

dW

through

dW\,W^{-1}

). The natural gradient step is

\Delta W=\eta\big(I-\mathbb E[\varphi(y)y^{\mathsf T}]\big)W,

which also needs no matrix inverse.

Equivariance. Write

C=WA

, so that

y=Cs

. The natural gradient gives

\Delta C=\eta\big(I-\mathbb E[\varphi(Cs)(Cs)^{\mathsf T}]\big)C

:

A

has disappeared. The plain gradient gives

\Delta C=\eta\big(C^{-\mathsf T}-\mathbb E[\varphi(y)s^{\mathsf T}]\big)A^{\mathsf T}A

, which depends on

A

through

A^{\mathsf T}A

.

The chapter also treats PCA, non-negative matrix factorization, sparse signal processing, convex programming and game-theoretic scores; this module takes ICA. Below: 400 samples of two uniform sources (sub-Gaussian, so $\varphi(y)=y^3$ is a suitable choice), batch updates with $\eta=0.05$ , and both runs start from the same global matrix $C_0=W_0A$ .

Explore

Mixing matrix A

Scatter shows

The mixtures

x=As

. Independent uniform sources would fill an axis-aligned square.

Separation error of

C=WA

(0 means a scaled permutation), log scale. natural plain

Worked example

Derive the equivariance property, and find where the algorithm settles for uniform sources.

Turn the plain gradient into the natural gradient.

For this invariant metric the natural gradient is $\nabla L\,W^{\mathsf T}W$ . So $\Delta W=-\eta\big(-W^{-\mathsf T}+\mathbb E[\varphi(y)x^{\mathsf T}]\big)W^{\mathsf T}W=\eta\big(I-\mathbb E[\varphi(y)y^{\mathsf T}]\big)W$ , using $x^{\mathsf T}W^{\mathsf T}=y^{\mathsf T}$ .
Multiply on the right by $A$ . What do you get?

$\Delta C=\Delta W\,A=\eta\big(I-\mathbb E[\varphi(y)y^{\mathsf T}]\big)C$ , and $y=Cs$ . Everything is in terms of $C$ and $s$ : the trajectory of $C$ depends only on where it starts.
Where does it stop? Take $y=Ds$ with $D$ diagonal and $s_i$ uniform with unit variance.

Off the diagonal, $\mathbb E[y_i^3y_j]=\mathbb E[y_i^3]\,\mathbb E[y_j]=0$ by independence. On the diagonal we need $\mathbb E[y_i^4]=1$ : $d_i^4\,\mathbb E[s_i^4]=1$ . For uniform $s$ on $[-\sqrt3,\sqrt3]$ , $\mathbb E s^4=9/5$ , so $d_i=(5/9)^{1/4}\approx0.863$ . Look at the diagonal of $C$ in the explorer after a run (the sample of 400 moves it slightly).
And the plain gradient?

$\Delta C=\eta\big(C^{-\mathsf T}-\mathbb E[\varphi(y)s^{\mathsf T}]\big)A^{\mathsf T}A$ . Same fixed points, but the step is preconditioned by $A^{\mathsf T}A$ . A badly conditioned $A$ makes it slow in one direction and unstable in another, which is what the "severe" setting shows.

Check your understanding

Even with perfect data, ICA can recover the sources only up to…

their order and their scales (including sign)
a rotationA rotation of independent non-Gaussian sources makes them dependent again, and ICA can detect that.
nothing: they are recovered exactly $As=(AP^{-1}\Lambda^{-1})(\Lambda Ps)$ for any permutation $P$ and diagonal $\Lambda$ : the data cannot tell these apart.

Rescaling or reordering independent sources keeps them independent, and the mixing matrix absorbs the change. Everything else is identifiable when at most one source is Gaussian.

Why can ICA not separate two Gaussian sources?

A rotation of independent standard Gaussians is again a pair of independent standard Gaussians
Gaussian signals have no fourth momentThey have every moment; the fourth is $3\sigma^4$ .
The natural gradient is undefined for GaussiansThe algorithm runs; the problem is that the answer is not unique.

After whitening, every rotation gives independent outputs, so independence cannot single out the true one. Non-Gaussianity is what breaks the rotational symmetry. The uniform sources here fill a square, which is clearly not rotation-invariant.

Equivariance of the natural-gradient rule says that the trajectory of $C=WA$ …

depends only on its starting value and the sources, not on $A$
depends on $A$ only through its condition numberThat describes the plain gradient, whose dynamics carry the factor $A^{\mathsf T}A$ .
does not depend on the sourcesThe update averages $\varphi(y)y^{\mathsf T}$ over $y=Cs$ , so the sources matter.

In the explorer the blue error curve is identical for all three mixing matrices; the red one changes. Because both runs start at the same

C_0

, the only difference between the settings is

A

itself.

Practice

The gradient of the log-determinantWhere W^(−T) comes from

Show that $\partial\log|\det W|/\partial W_{ij}=(W^{-\mathsf T})_{ij}$ , and hence derive $\nabla L$ above.

Solution

$d\log|\det W|=\mathrm{tr}(W^{-1}dW)=\sum_{ij}(W^{-1})_{ji}\,dW_{ij}$ , so the gradient is $(W^{-1})^{\mathsf T}$ . For the other term, $\partial(-\log q_i(y_i))/\partial W_{ij}=\varphi_i(y_i)\,x_j$ , which gives $\mathbb E[\varphi(y)x^{\mathsf T}]$ .

Super-Gaussian sourcesChoosing the nonlinearity

For Laplace sources, density $\propto e^{-\sqrt2|s|}$ (unit variance), what is the score $\varphi=-(\log q)'$ if $q$ is the true density? Why do practical algorithms use $\tanh$ instead?

Solution

$-\log q=\sqrt2|s|+\text{const}$ , so $\varphi(s)=\sqrt2\,\mathrm{sign}(s)$ . It is discontinuous at 0, which makes the update rough; $\tanh(\beta y)$ is a smooth approximation of a sign function and has the same qualitative effect for super-Gaussian sources. (For sub-Gaussian sources, as in the explorer, a cubic such as $y^3$ is the usual choice.)