Best approximation of sum of unit vectors by a smaller subset

Let $A$ denote the matrix whose columns are $v_1,\dots,v_N$. Then your problem is that of minimizing $\|s_N - Ax\|$ subject to the constraint that $x$ has $0,1$ entries and $\|x\| \leq \sqrt{n}$.

Removing the constraint that $x$ has $0,1$ entries leaves us with a much easier problem to deal with. I suspect that its solution will yield a useful heuristic.

If $A = U \Sigma V^T$ is an SVD and we make the substitutions $b = U^Ts_N$ and $y = V^Tx$, we are left with the simplified problem $$ \min \|\Sigma y - b\| \quad \text{s.t. } \quad \|y\| \leq \sqrt{n}. $$ This is easily solved with Lagrange multipliers. The squared objective and constraint functions have the forms $$ f(y) = \|\Sigma y - b\|^2 \implies \nabla f = 2 [\Sigma^2 y - \Sigma b] \\g(y) = \|y\|^2 \implies \nabla g = 2y $$ So, we have $$ \nabla f = \lambda \nabla g \implies \Sigma^2 y - \Sigma b = \lambda y \implies (\Sigma^2 - I)y = \lambda \Sigma b \implies y = \lambda(\Sigma^2 - I)^{-1}\Sigma b. $$ Note: this assumes that $A$ does not have $1$ as a singular value, which occurs with probability $1$. Plugging into the constraint yields $$ \|\lambda(\Sigma^2 - I)^{-1}\Sigma b\|^2 = n \implies \lambda = \pm \sqrt{\frac{n}{\|(\Sigma^2 - I)^{-1}\Sigma b\|^2}}, $$ which is simply to say that this solution for $y$ should be normalized to the radius-$\sqrt{n}$ sphere.

I'm not sure if this can be written in terms that remove the SVD. For what it's worth, though, we have $$ (\Sigma^2 - I)^{-1}\Sigma = V^T[(A^TA - I)^{-1}\sqrt{A^TA}]V. $$