Signal processing and optimization, by hand

Companion to the notes on Chapter 13 of Amari's Information Geometry and Its Applications. Back to the notes. Eight small experiments, one for each claim of the chapter that is worth seeing move. They use the same mathematics as code/signal_processing.py; the default readouts of widgets 1, 2, 3 (top table), 4, 5, 6 and 7 are the numbers that script prints. Widgets 3 (the simulation button) and 8 draw their own seeded random numbers, so their simulated numbers agree with the script only up to Monte Carlo error; widget 4 starts from the script's own source samples.

1. Is the Stiefel manifold stable under the principal-component flows?

Take $k=2$ orthonormal vectors in $\mathbb R^6$, collected in $W$ with $W^{\mathsf T}W=I$, and let the flow act on a slightly larger matrix $M$ that has been pushed off the manifold. The defect $\lVert M^{\mathsf T}M-I\rVert$ tells whether the flow brings $M$ back, leaves it where it is, or throws it away. Oja's flow (13.31) is stable; the ascent flow (13.34) for minor components blows up; the flow (13.35) of Chen et al. conserves $M^{\mathsf T}M$ exactly, so the defect neither grows nor decays; the sign of the second term of (13.36) as printed is wrong, and the printed flow shrinks $M$ to zero, while the sign-corrected version conserves $M^{\mathsf T}M$ like (13.35); (13.40) and (13.41) pull $M^{\mathsf T}M$ back to $D^2$. Change the push (slider) and compare the last column of the readout, the distance of the final subspace from the principal and the minor subspace of $V_X$.

2. Second-order statistics leave a rotation free; independence picks it

Two independent uniform sources are mixed and then whitened by PCA, as in the book's example (13.44)–(13.46). The result is the source pair up to a rotation by some angle, and no covariance can tell which. Turn the rotation: the outputs stay uncorrelated for every angle, but they are independent only at $0^\circ$ and $90^\circ$. The mutual information of the pair has the closed form $I=2\ln a+b/a$ with $a=\max(\lvert\cos\theta\rvert,\lvert\sin\theta\rvert)$ and $b=\min(\ldots)$ (maximum $1-\ln2$ at $45^\circ$). The fourth-order statistics see the angle: the excess kurtosis of $y_1$ is $k(\cos^4\theta+\sin^4\theta)$ and the correlation of the squares is nonzero. The button puts the angle where the principal components of the book's example sit. Switch the sources to Gaussian and the picture, and every number, stops depending on the angle: this is the indefiniteness of rotation of (13.19).

3. When does the separating point attract? The sign of one eigenvalue

For the natural-gradient rule $\dot C=(I-\varphi(y)y^{\mathsf T})C$ the Hessian of the loss at the separating point splits into $1\times1$ blocks $E[\varphi_i'(s_i)s_i^2]+1$ and $2\times2$ blocks $\left[\begin{smallmatrix}m_i\sigma_j^2&1\\1&m_j\sigma_i^2\end{smallmatrix}\right]$ (with $m_i=E\varphi_i'(s_i)$, sources rescaled so that $E[\varphi(s)s]=1$). With $\kappa_i=m_i\sigma_i^2$ the point is a minimum exactly when $\kappa_1\kappa_2>1$; the weak direction is the antisymmetric pair, a rotation of the outputs. For $\varphi=y^3$ this reads $(3+k_1)(3+k_2)<9$ in terms of the excess kurtoses. Pick two sources and a nonlinearity: a uniform pair is stable with $y^3$ and unstable with $\tanh$; a uniform and a Laplace source are unstable with $y^3$; a uniform and a sech source (kurtosis 2) sit exactly on the boundary. The button runs 3000 paired online chains started at $C=I\pm\delta v$ along the weakest eigenvector and compares the mean response at $t=2$ with $e^{-h_{\min}t}$ (the step $\varepsilon=0.005$ leaves an $O(\varepsilon)$ gap of a few percent, which the script removes by extrapolation).

4. The natural gradient does not care what the mixing matrix is

Write $C=WA$. For the natural gradient $W\leftarrow W+\varepsilon(I-\varphi(y)y^{\mathsf T})W$ the update of $C$ is $C\leftarrow C+\varepsilon(I-\varphi(Cs)(Cs)^{\mathsf T})C$ and contains $A$ nowhere: with the same source samples and the same $C(0)$, two different mixing matrices give the same trajectory of $C$, to rounding error (relative gradient, Cardoso and Laheld). The ordinary gradient $W\leftarrow W+\varepsilon(W^{-\mathsf T}-\varphi(y)x^{\mathsf T})$ gives $\Delta C=\varepsilon(C^{-\mathsf T}-\varphi(y)s^{\mathsf T})A^{\mathsf T}A$: the speed, and eventually the stability, depend on $A$. Uniform sources, $\varphi=y^3$, $\varepsilon=0.01$. Move the first entry of $A_2$ and watch the ordinary gradient (orange, green) change while the natural gradient (blue, dashed) does not. (The default samples are the script's own, so the default readout equals the script's; the second button draws a fresh random stream here, for which the numbers differ but the identity of the two natural trajectories does not.)

5. Lee–Seung is exponential gradient with a particular step

The multiplicative update (13.124)–(13.125) is presented as the analogue of the exponential gradient (13.121). It is exactly an exponential-gradient step $A\leftarrow A\,e^{-\eta_{it}(\partial L/\partial A)_{it}}$ whose rate is $\eta_{it}=1/\mathrm{logmean}(N_{it},P_{it})$ (with $\partial L/\partial A=P-N$), so it is not (13.121) for any single learning constant. It decreases the loss at every iteration, while exponential gradient with a fixed $\eta$ does so only for small $\eta$. The data are the matrices of the script ($X=AS$ with $A$ of size $6\times3$ and $S$ of size $3\times40$, 300 iterations from a common start).

6. Least angle: the inactive axes can be the closer ones

Along the lasso path the dual coordinates satisfy $\eta_A=-\lambda s_A$ on the active set, so $\lvert\langle\dot\theta,e_i\rangle\rvert$ is the same for all active $i$ (Theorem 13.4 (1)). The book adds that the inactive axes are at larger angles. With $d=G_{AA}^{-1}s_A$ the inactive inner products are $g_j=(G_{jA}d)$, per unit inner product of the active ones, and nothing keeps $\lvert g_j\rvert\le1$: that is the irrepresentable condition of lasso theory, and it fails often. Three unit columns in $\mathbb R^3$: the first two have correlation $\rho$ (set below), the third lies almost in their span. The table lists every segment of the exact path for the data of the script; at $\rho=-0.9$ the segment with active set $\{1,2\}$ has an inactive axis with $\lvert g_3\rvert=4.27$ and makes an angle of $17^\circ$ with the path, against $77^\circ$ for the active axes.

7. The central path is straight in the dual coordinates

For a linear program over a polygon, minimise $c\cdot\theta$ with the barrier $\psi=-\sum_k\log(A_k\theta-b_k)$: the minimiser of $tc\cdot\theta+\psi$ satisfies $\nabla\psi(\theta^*(t))=-tc$, a straight line in $\eta=\nabla\psi$ (an m-geodesic from the centre) and a curved path in $\theta$. The dual vector $y_k=1/(ts_k)$ is feasible and the gap is exactly $m/t$. The natural-gradient (affine-scaling) step $\theta\leftarrow\theta-\varepsilon G^{-1}c$ follows this path from the analytic centre with an error that is first order in $\varepsilon$: the path's curvature in the $\theta$ chart is what limits the step. Move $t$ along the path and change the step size; the table gives the error of the Euler iterates at $t=2$.

8. The Hyvärinen score is efficient only for the Gaussian

The score-matching estimator solves $\sum_i s(x_i,\theta)=0$ with $s=\partial_\theta(\ddot l+\tfrac12\dot l^2)$, needs no normalising constant and is consistent, but its asymptotic covariance $K^{-1}VK^{-1}$ exceeds the Cramér–Rao bound unless $s$ lies in the span of the score. For $p\propto\exp(-\theta\lvert x\rvert^k)$ the efficiency is a ratio of Gamma functions: $1$ at $k=2$ (Gaussian), $0.8214$ at $k=3$ (the book's Example 13.1, whose printed $\psi$ and $s$ have sign slips) and smaller for larger $k$. The simulation draws $\theta x^k\sim\mathrm{Gamma}(1/k)$ exactly, with $\hat\theta_{\rm MLE}=N\theta/(k\sum y_i)$ and $\hat\theta_{\rm H}=\frac{k-1}{k}\sum x_i^{k-2}/\sum x_i^{2k-2}$.