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$


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

