If you use constrained optimization, you do get linear equations. I set $g(x)=||x||^2$. For every $x \in \mathbb{R}^3$, $$\nabla f(x) = \sum_i 2(x \cdot n_i)n_i \text{ and } \nabla g(x) = 2x.$$ You look at (unit) vectors $x$ such that $\nabla f(x) = \lambda \nabla g(x)$ for some real number $\lambda$, namely you look at eigenvectors of the symmetric endomorphism $u$ given by $$u(x) = \sum_i (x \cdot n_i)n_i.$$ Actually, using an orthogonal basis of eigenvectors of $u$, you get that the minimum of $f(x) = x \cdot u(x)$ over all $x \in \mathbb{S}_2$ is achieved when $x$ is an eigenvector associated to the least eigenvalue of $u$.

If you use constrained optimization, you do get linear equations. I set $g(x)=||x||^2$. For every $x \in \mathbb{R}^3$, $$\nabla f(x) = \sum_i 2(x \cdot n_i)n_i \text{ and } \nabla g(x) = 2x.$$ You look at (unit) vectors $x$ such that $\nabla f(x) = \lambda \nabla g(x)$ for some real number $\lambda$, namely you look at eigenvectors of the symmetric endomorphism $u$ given by $$u(x) = \sum_i (x \cdot n_i)n_i.$$ Actually, using an orthogonal basis of eigenvectors of $u$, you get that the minimum of $f(x) = x \cdot u(x)$ over all $x \in \mathbb{S}_2$ is achieved when $x$ is an eigenvector associated to the least eigenvalue of $u$.

This can be proved directly. Call $\lambda_1 \le \lambda_2 \le \lambda_3$ the eigenvalues of $u$ and call $(e_1,e_2,e_3)$ an orthonormal basis of associated eigenvectors. Then for every $x = \xi_1e_1+\xi_2e_2+\xi_1e_3 \in \mathbb{R}^3$, $$f(x) = x \cdot u(x) = \sum_{i=1}^3 \lambda_i\xi^2 \ge \sum_{i=1}^3 \lambda_1\xi^2 = \lambda_1||x||^2,$$and equality holds when $x$ is a multiple of $e_1$.