Tonga: low rank approximation of the natural gradient, question regarding Le Roux et al. 2007

This question concerns an implementation of the topmoumoute natural gradient (tonga) algorithm as described in page 5 in the paper Le Roux et al 2007 http://research.microsoft.com/pubs/64644/tonga.pdf.

I understand that the basic idea is to augment stochastic gradient ascent with the covariance of the stochastic gradient estimates. Basically, the natural gradient approach multiplies a stochastic gradient with the inverse of the covariance of the gradient estimates in order to weight each component of the gradient by the variance of this component. We prefer moving into directions that show less variance during the stochastic gradient estimates:

$ng \propto C^{-1} g$

Since updating and inverting the covariance in an online optimisation setting is costly, the authors describe a cheap'' approximate update algorithm as described on page 5 as:

$C_t = \gamma \hat{C}_{t-1} + g_tg_t^T$ where $C_{t-1}$ is the low rank approximation at time step t-1. Writing $\hat{C}_{t} = X_tX_t^T$ with $X_t =\sqrt{\gamma} X_{t-1}\ \ g_t]$ they use an iterative update rule for the gram matrix $X_t^T X_T = G_t =$

($\gamma G_{t-1}$      $\sqrt{\gamma} X^T_{t-1}g_t$

$\sqrt{\gamma} g^T_t X_{t-1}$      $g_t^Tg_t$)

They then state To keep a low-rank estimate of $\hat{C}_{t} = X_tX_t^T$, we can compute its eigendecomposition and keep only the first k eigenvectors. This can be made low cost using its relation to that of the Gram matrix

$G_t= X_t^T X_T$:

$G_t = VDV^T$

$C_t = (X_tVD^{-\frac12})D(X_tVD^{-\frac12})^T$''

Because it's cheaper than updating and decomposing G at every step, they then suggest that you should update X for several steps using $C_{t+b} = X_{t+b}X_{t+b}^T$ with $X_{t+b} = \left[\gamma U_t, \ \gamma^{\frac{b-1}{2}g_{t+1}},...\ \ \gamma^{-\frac12}g_{t+b-1}, \ \gamma^{\frac{t+b}{2}g_{t+b}}\right]$

I can see why you can get $C_t$ from $G_t$ using the eigendecomposition. But I'm unsure about their update rule for X. The authors don't explain where U is coming from. I assume (by notation) this is the first k eigenvectors of $C_t$, correct? But if so, why would the formula for $X_t$ be a good approximation for $X_t$? When I implement this update rule, $X_t$ does not seem to be a good approximation of the ''real'' $X_t$ (that you would get from $X_t = [\sqrt{\gamma} X_{t-1}\ \ g_t]$) at all. So why should I then be able to get a good approximation of $C^{-1}$ from this (I don't)? The authors are also not quite clear about how they keep $G_t$ from growing (The size oft he matrix $G_t$ increases at each iterative update). I assume they replace $G_t$ by $\hat V\hat D\hat V^T$ with $\hat V$ and $\hat D$ being the first k components of the eigendecomposition?

So in summary:

• I tried implementing this update rule, but I'm not getting good results and am unsure my implementation is correct

• why should the update rule for $X_t$ be reasonable? Is $U_t$ really the first k eigenvectors of $C_t$? (Clearly I cannot let $X_t$ and $G_t$ grow for each observed gradient $g_t$)

• This is far fetched, but has anyone implemented this low rank approximate update of the covariance before and has some code to share so I can compare it to my implementation?

As a simple example if I simulate having 15 gradients in matlab like:

X = rand(15, 5);
c = cov(X);
e = eig(c);
%if the eigenvectos of C would give a good approximation of the original , this should give an appoximation of c, no?:
c_r = e'*e; %no resemblence to c


So I'm quite certainly doing it wrong, I guess U might actually not be the eigenvectors of C, but then what is U? Any suggestions or references would be most welcome!

(sorry about the terrible layout, looks like only a subset of latex is supported, no arrays for matriced and embedded latex doesn't look particularly good, all the formulas are much more readable on page 5 of the referenced paper :) Also, is this considered off-topic here? It's really more related to optimisation and machine-learning...)

So it turns out the first assumption was actually correct: U is indeed the first k eigenvectors of C, that we calculate from G by means of the eigendecompposition $(X_tVD^{-\frac12})D^{1/2} = X_tV$.