I am trying to find a simple, ideally closed form formula for the (not necessarily unique) unit vector $\vec{x}$ minimizing total squared cosine distance from a collection of unit vectors $\vec{v_i}$. In other words I want to solve
$$\operatorname*{argmin}_{||\vec{x}|| = 1} \sum_i (1 - \vec{v_i} \cdot \vec{x})^2.$$
I have attempted to solve this problem using Lagrange multipliers with objective function
$$f(\vec{x}, \lambda) = \sum_i \left(1 - \frac{\vec{v_i} \cdot \vec{x}}{||\vec{x}||}\right)^2 + \lambda (1 - ||\vec{x}||^2).$$
But the system of equations that I get is not easy to work with:
$$\vec{v} = \left( \left(\sum_i \vec{v_i} \otimes \vec{v_i} \right) + ( \vec{v} \cdot \vec{x} - (\vec{v} \cdot \vec{x})^2 - \lambda ) I\right) \vec{x}$$ and $$||\vec{x}|| = 1$$
where $\vec{v} = \sum_i \vec{v_i}$. In contrast, without the constraint $||\vec{x}|| = 1$ the problem has a beautiful set of critical points:
$$\vec{v} = \left( \sum_i \vec{v_i} \otimes \vec{v_i} \right) \vec{x}.$$
At this point I am having trouble proceeding. I could stop here and use numerical methods to find approximate solutions but can this system of equations be further simplified or is there a better approach?
I would diagonalize the covariance matrix $\sum_i v_i^Tv_i$ first. That matrix is symmetric (and positive semidefinite), so it has an orthogonal eigenbasis. The eigendirections can be found by some standard numeric methods.
After the orthonormal basis transform, the matrix $\sum_i v_i^Tv_i$ is diagonal with consisting of the eigenvalues along the diagonal; let $c_j$ be the $j$th diagonal element. It is reasonable to compute the vector $\sum v_i$ as well; let $s=\sum_i v_i$.
After that, your system of equations will be $$ (c_j-\lambda)x_j=s_j \quad (j=1,2,\ldots,d); \qquad \sum_{j=1}^d x_i^2=1 $$ where $d$ is the dimension of the space we work in. Hence, we have to solve the $2d$-degree equation $$ \sum_{j=1}^d \left(\frac{s_j}{c_j-\lambda}\right)^2 =1 $$ for $\lambda$, then we obtain $$ x_j=\frac{s_j}{c_j-\lambda} \quad(j=1,\ldots,d). $$