Suppose that I have a set of $M$ objects, each with the mass of $m_i$ and at the position $\mathbf r_i$ from the point $\mathbf 0 \in\Bbb R^3$. Assuming also that the center of mass of this group of object is at the origin, i.e. $$ \sum_i m_i\mathbf r_i = \mathbf 0. $$
Suppose that this group of objects almost lie in a plane (in some certain sense). I want to find an "approximating normal vector" $\mathbf n \in\Bbb S^1$ to this approximated plane. My idea is that $\mathbf n $ should solve the following constrained minimization problem: $$\begin{align} &\text{minimize} &f(\mathbf n) = \sum_i m_i(\mathbf r_i\cdot \mathbf n)^2 \\ &\text{subject to} &g(\mathbf n) = |\mathbf n| = 1. \end{align}$$ The thought process is that $f(\mathbf n)$ measures the deviation of $\mathbf n$ from being the "true" normal vector to the plane. Of course, if all of our objects actually lie in a plane then the solution to this problem is the unique (up to a change of sign) normal vector to this plane.
Is there an explicit solution to this constrained minimization problem? Namely, can we write $\mathbf n$ as an algebraic expression of $m_i$'s and $\mathbf r_i$'s?
I tried to solve it using the Lagrange multiplier method but I don't know how to algebraically solve a system of quadratic polynomials in $3$ variables.
The reason for asking the above question is that, for my application, each $\mathbf r_i$ is actually a time-dependent function $\mathbf r_i(t)$ so I can't just solve the problem numerically for every single $t$ because it would take too long. Alternatively, an efficient numerical algorithm that would be suitable for the time-dependent version of this problem would be very welcomed as well (even with different weighting $|\cdot|^p$ instead of $(\cdot)^2$ up there).
The problem as posted is only geometrical so the mass associated to each point is irrelevant. Now given a set of points $p_k,\ k=1,\cdots,n$ and given a plane $\Pi\to (p-p_0)\cdot \vec n = 0$ we have $(p_k-p_0)\cdot \vec n = \delta_k$ as a deviation from the exact inclusion $p_k\in \Pi$ so a good procedure to find such a plane is by solving the minimization problem
$$ \min_{p_0,\vec n}\sum_{k=1}^n\delta_k^2 = \min_{p_0,\vec n}\sum_{k=1}^n\left((p_k-p_0)\cdot \vec n\right)^2 \ \ \text{s. t.}\ \ ||\vec n|| = 1 $$
which is the typical quadratic regression problem.