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


$\text{trace}(Y Y')$


$Y=X - B.X$

subject to

$B\odot I = 0$

  • 1
    $\begingroup$ why do you need this? $\endgroup$ – Aksakal May 11 '17 at 2:16
  • $\begingroup$ 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. $\endgroup$ – 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} $

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.