Simultaneously finding least squares predictors of every feature based on all others Given set of $d$ datapoints and $f$ features in a $X_{f,d}$ data matrix, I'm trying to fit a least squares predictor of every feature based on every other feature.
In other words for every feature $i$ I want vector $\mathbf{b_i} = b_1, \ldots, b_{i-1}, 0, b_{i+1}, \ldots, b_f$ such that $\mathbf{b_i} X$ approximates $X_{i,\cdot}$ in least squares sense
In my case $f$ is around 10,000 and $d$ is millions so I need some kind of approximate/iterative solution. Is this a well known problem?
I tried minimizing following formulation by doing gradient descent, and setting setting diagonal entries of $B$ to zero at each iteration to satisfy the constraint, but getting it to converge was tricky and I feel like I'm reinventing the wheel
minimize
$\text{trace}(Y Y')$
where
$Y=X - B.X$ 
subject to
$B\odot I = 0$
 A: Maybe this is useful at reusing some stuff? 
Warning: I haven't thought too much about the numerical stability of this etc... I don't know if this at all makes sense on a problem of your magnitude. 
Some linear algebra math
Let $S = \begin{bmatrix} X & \mathbf{z}\end{bmatrix}'\begin{bmatrix} X & \mathbf{z} \end{bmatrix}$. If $ S = \begin{bmatrix} S_{xx} & \mathbf{s}_{xz} \\ \mathbf{s}_{xz}' & s_{zz} \end{bmatrix} $ and $ S^{-1} = \begin{bmatrix} A & \mathbf{u} \\ \mathbf{u}'  &  c \end{bmatrix} $ then:
$$S_{xx}^{-1} = A - \frac{\mathbf{u}\mathbf{u}'}{c}$$
And projection of $\mathbf{z}$ onto columns of $X$ is given by:
\begin{align*}
(X'X)^{-1} X'\mathbf{z} &= S^{-1}_{xx} \mathbf{s}_{xz}\\ 
\end{align*}
Hence maybe you could do the following?
You have some 1 million by 10,000 data matrix W.


*

*Compute $S = W'W$ (sounds rough but doable)

*Compute $S^{-1}$ using something like singular value decomposition (I'd guess takes roughly 5 minutes if 10,000 by 10,000?). You might have to do something with too small of singular values? (Does pseudoinverse work for this? I'd have to think...)

*Iterate over $i=1,\ldots,10000$ and


*

*Take the $i$th row and column out of $S^{-1}$ to get matrix $A_{-i}$

*Use the $i$th column to get vector $\mathbf{u}_{i}$ (which is missing $i$th row) and scalar $c_i$ (which is $i$th row of $i$th column)

*Compute $\left( A_{-i} - \frac{\mathbf{u}_{i} \mathbf{u}_i'}{c_i} \right) S_{X_{-i}\mathbf{x}_i} $


A: When $x$ have non-singular second-moment matix, we can get coefficients matrix $B$ from $(X'X)^{-1}$ as follows:
Let $D2$ be the inverse second moment with off-diagonal terms set to zero
$$D2=[(X'X)^{-1}]_d$$
Then coefficients $\{a_j\}$ for the model $x_i=a_0 x_0 + \ldots a_n x_n$ can be read from $i$th row of following matrix
$$B=I-(X'X)^{-1} D2^{-1}$$
https://colab.research.google.com/drive/1tXmS6kZcJ8hsfvzKaZSCvXf7iNRT25UG#scrollTo=mTMmob1gJ9Zc
