Range of $a$ such that $w \leftarrow w-a x \langle w, x \rangle$ converges almost surely? Edit Sep 19 this answer on Mathoverflow matches simulation results
Suppose $x_i$ come from 2-d standard Normal centered at 0. What is the range of $a$ for which the following iteration converges almost surely?
$$w_{i+1} = w_i-a x_i \langle w_i, x_i \rangle$$
Mean norm of $w$ diverges for $a>2/3$ but this is vacuous because variance also diverges. The answer is somewhere between 0.9 and 1.1.
How do I get a more accurate estimate?

Here's a way to verify that $a=2$ converges for 1D, which doesn't generalize to 2D
Note that our process defines the following recursion on $w^2$
$$\log w_{k+1}^2 = \log (1-\alpha x^2)^2 + \log w_k^2$$
Here we see that $\log w_{k+1}^2$ is a sum of $k$ independent random variables, tending towards $-\infty$ as $k$ grows whenever $E[\log (1-\alpha x^2)^2]<0$, from Central Limit Theorem.
Since $E[\log (1-\alpha x^2)^2]\approx 0$ for $\alpha \approx
2.42129$, this also forms the critical rate, and can be verified in simulations by looking at median value of trajectories for $\alpha$ slightly above, and slightly below critical

Edit Aug 24
Expectation in answer has a closed form solution in terms of exponential integral, can be obtained in Mathematica. Plotting it for $d=2$ there are two critical points, for $a=0.121312$ and for $a=0.918611$

 A: Your iterations can be rewritten as $w_n \leftarrow w_{n-1} - axx^\top w$, or equivalently
\begin{equation*}
w_n \leftarrow (I - axx^\top)w_{n-1}
\end{equation*}
and so can be seen to be a form of power iteration. Therefore, your iterations converge to the dominant eigenvector of $I - axx^\top$. It is straightforward to calculate that for $x = (x_1, x_2)$, the eigenvalues are $\lambda_1 = 1$ and $\lambda_2 = 1 - ax_1^2 - ax_2^2$ with corresponding eigenvectors $(-x_2/x_1, 1)$ and $(x_1/x_2, 1)$ respectively.
If $\lambda_1 = 1$ is the dominant eigenvalue, then your iterations will converge, since $w_n \rightarrow c\cdot(-x_2/x_1, 1)$ for some finite constant c (that depends on the initial $w_0$). Next, note that $\lambda_2 \leq 1$, and so is the dominant eigenvalue if and only if $\lambda_2 < -1$. This happens iff $x_1^2 + x_2^2 > 2/a$, which occurs with probability $\Pr(\chi^2_2 > 2/a)$. But in this event, $|\lambda_2| > 1$, and so $|w_n| \rightarrow \infty$ as $n\rightarrow \infty$ and your iterations diverge.
In summary, your iteration scheme diverges with probability $\Pr(\chi^2_2 > 2/a)$ for almost all initializations of $w_0$ (with the exception of $w_0 \in \operatorname{span}(\{(-x_2/x_1, 1)\})$, which of course happens with probability $0$).

Below is some Python 3 code that can help verify these findings. One can try various values of $a$ and find that the expected number of divergences is fairly close to the expected number from the chi-square distribution.
from scipy.stats import multivariate_normal
from scipy.stats import chi2
import numpy as np

np.random.seed(0)

a = 0.9
num_diverged = 0
for i in range(1000):
    x = multivariate_normal.rvs([0, 0])

    w = 20 * multivariate_normal.rvs([0, 0]) # a random initialization for w

    for _ in range(1000):
        w = w - a * w.dot(x) * x
        #w = w/(w[0]**2 + w[1]**2)**0.5  # Uncomment to see normalized w's

    print()
    print("Final w:              ", w)
    if np.isnan(sum(w)) or abs(w[0]) > 100 or abs(w[1]) > 100:
        num_diverged += 1
    print("Normalized w:         ", w/np.linalg.norm(w))
    if 2 <= a * x[0]**2 + a* x[1]**2:
        x = np.array([x[0]/x[1], 1])
    else:
        x = np.array([-x[1]/x[0], 1])
    print("Expected normalized w:", x/np.linalg.norm(x))

print()
print("Number diverged:", num_diverged)
print("Expected diverged:", chi2.sf(2/a, 2)*1000)

Partial output:
Number diverged: 335
Expected diverged: 329.1929878079054

A: $\newcommand{\E}{\operatorname{E}}$
This is a new answer for the updated question, where we have a sequence of Gaussians instead of a single $x$.
We have that $w_n = \left[\prod_{i=1}^n (I - ax_ix_i^\top)\right]w_0$, and want to examine the behavior as $n\rightarrow \infty$. We're going to use the trick that an infinite product converges to zero iff the product of squares converges to zero:
\begin{align*}
\lim_{n\rightarrow\infty}\prod_{i=1}^n(I - ax_ix_i^\top)^2 
&= \lim_{n\rightarrow\infty}\exp\left(n\cdot \frac{1}{n}\sum_{i=1}^n\log(I - ax_ix_i^\top)^2\right) \\
&\overset{a.s.}{\longrightarrow}\lim_{n\rightarrow\infty}\exp\left(n\cdot \E\left[\log(I - ax_ix_i^\top)^2\right]\right)
\end{align*}
Above, we use the matrix exponential and matrix logarithm. Using the spectral decomposition
\begin{equation*}
(I - ax_ix_i^\top)^2 = \begin{bmatrix}
x_{i,2}/x_{i,1} & x_{i, 1}/x_{i,2} \\ 1 & 1
\end{bmatrix}
\begin{bmatrix}
1 & 0 \\ 0 & (1-a||x_i||^2)^2
\end{bmatrix}
\begin{bmatrix}
-\frac{x_{1,i}x_{2,i}}{||x_i||^2} & \frac{x_{1,i}^2}{||x_i||^2} \\ \frac{x_{1,i}x_{2,i}}{||x_i||^2} & \frac{x_{2,i}^2}{||x_i||^2} 
\end{bmatrix}
\end{equation*}
we can calculate that
\begin{equation*}
\E\left[\log(I - ax_ix_i^\top)^2\right] = \begin{bmatrix}
\E\left[\frac{x_{1,i}^2}{||x_i||^2}\cdot\log(1-a||x_i||^2)^2\right] & 0 \\ 0 & \E\left[\frac{x_{2,i}^2}{||x_i||^2}\cdot\log(1-a||x_i||^2)^2\right]
\end{bmatrix}
\end{equation*}
Thus, we are concerned with when the following matrix has limit zero:
\begin{equation*}
\lim_{n\rightarrow\infty}\begin{bmatrix}
\exp\left(\E\left[\frac{x_{1,i}^2}{||x_i||^2}\cdot\log(1-a||x_i||^2)^2\right]\right) & 0 \\ 0 & \exp\left(\E\left[\frac{x_{2,i}^2}{||x_i||^2}\cdot\log(1-a||x_i||^2)^2\right]\right)
\end{bmatrix}^n
\end{equation*}
and this of course happens iff it has norm less than 1. That is,
\begin{equation*}
\sqrt{2\cdot \left[\exp\left(\E\left[\frac{x_{1,i}^2}{||x_i||^2}\cdot\log(1-a||x_i||^2)^2\right]\right)\right]^2} < 1
\end{equation*}
Some numerical integration suggests that the critical value is roughly $a \approx 0.92$.
