# 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$

• why do you need this? – Aksakal May 11 '17 at 2:16
• This seems to be what the Amelia package does, in the context of doing multiple imputation of missing data. Check out the Blackwell et al paper linked here: gking.harvard.edu/Amelia and take a look at the material around page 8 / equation 2. I've only used this for imputation and have not investigated this model in detail. – zbicyclist May 11 '17 at 3:22

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.

1. Compute $S = W'W$ (sounds rough but doable)
2. 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...)
3. 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}$