I think it should be possible to vectorize this loop (i.e. not use a Python `for` loop), with something like:
```python
diff = X.T - y
normdiff = norm(diff, axis=1)
mask = normdiff >= 1e-6
if mask.sum() < X.shape[1]:
    eta = 1.
diff = diff[mask, :]
normdiff = normdiff[mask]
res = np.sum(diff / normdiff, axis=0)
T_denom = np.sum(1 / normdiff, axis=0)
T_nom = np.sum(X / normdiff, axis=1)
```
(warning: code untested)