Apologies in advance if this has already been asked before (I suspect it has, but I'm not experienced to know what to call it, or how to classify this problem). Given a set of $m$ points in space,

$$\mathbf{X}=\left( \begin{array}{cccc} x_1 & x_2 & \text{...} & x_m \\ y_1 & y_2 & \text{...} & y_m \\ z_1 & z_2 & \text{...} & z_m \\ \end{array} \right)$$

and a set of $N$ numbers $S=\{f_1,...,f_N\}$ with $N>m$, I would like to find the $m$-element ordered subset $S_\text{opt}=\left\{f_{\sigma (1)},f_{\sigma (2)},\text{...},f_{\sigma (m)}\right\}$ of $S$ (where $\sigma:\{1,...,m\}\rightarrow\{1,...,N\}$ is an injective function) such that

$$\mathbf{Y}(\sigma)=\mathbf{X}+\left( \begin{array}{cccc} 0 & 0 & \text{...} & 0 \\ 0 & 0 & \text{...} & 0 \\ f_{\sigma (1)} & f_{\sigma (2)} & \text{...} & f_{\sigma (m)} \\ \end{array} \right)$$

"most closely" resembles a plane. To be precise, I would like $\sigma$ to minimize the RMS error between the plane of best linear fit and the set of points $\mathbf{Y}(\sigma)$.

I know from linear algebra that given a set of points in 3-space $\mathbf{Y}(\sigma)$, the RMS error of the plane of best fit is given by


where $\Sigma_3$ is the 3rd singular value of $\mathbf{Y}(\sigma)-\mu(\mathbf{Y})$ (where the singular values are ordered from largest to smallest), where $\mu(\mathbf{Y})$ is the centroid of the points $\mathbf{Y}(\sigma)$.

I also know that the number of possible choices of $\sigma$ is given by $$\frac{N!}{(N-m)!}$$ which means that brute-force enumeration is out of the question (for the application I am looking to use this in, $N=47$ and $m=19$, giving 848255586288593837691617280000 combinations).

  • Is there some method of reducing the size of the computation space for this sort of problem? Initially I briefly thought that it might be a simple linear programming problem, since we're finding an optimal choice of members of $S$ to make $\mathbf{Y}=\mathbf{X}+S_\text{opt}$ as flat as possible, but obviously the problem isn't linear... right?

  • What general class of optimization does this sort of procedure fall under? I am somewhat inexperienced in scientific computing, and thus have no idea what heading this falls under.

If it is any help, the initial set $\mathbf{X}$ already closely resembles points in an $xy$ plane, and the $f_{\sigma(k)}$ are designed to be tiny $z$-axis height corrections to make the resulting set of points become more closely planar.

  • $\begingroup$ Does the plane of best linear fit mean ordinary least squares? And is the RMS error computed as $\sum_k |z_k'-z_k|^2$? $\endgroup$
    – Kirill
    Commented Oct 21, 2014 at 1:39
  • $\begingroup$ Is there any motivation for why the $f$'s cannot be reused? $\endgroup$
    – Kirill
    Commented Oct 21, 2014 at 2:10
  • $\begingroup$ @Kirill: In the case I am looking at, the $f_k$ physically correspond to a set of metal shims which have been precisely measured as to their thicknesses. In a larger-scale engineering application, you would simply measure a large enough number of shims and bin them so that you can simply choose $f_k=g_k$ (in the notation of your below answer) and this whole problem would be trivial, but for the moment there are only around 40 precisely measured samples, so finding optimal assignment of shims to locations on the device is somewhat trickier. $\endgroup$ Commented Oct 21, 2014 at 12:58

1 Answer 1


This looks like a very programming-contest-type problem. Let the OLS-fit plane be $z=\alpha x+\beta y$, so that we would like to minimize $$ \sum_k \big(z_k + f_{\sigma(k)} - \alpha x_k - \beta y_k\big)^2, $$ over $\sigma$. After computing $\alpha,\beta$, and by defining $g_k$ appropriately, this is the same as minimizing $$ \sum_k \big(f_{\sigma(k)} - g_k\big)^2 $$ over $\sigma$, where $g_k$ is a given set of real numbers.

The point of writing it this way is that this makes it clear that this is the optimal assignment problem, in a graph where each node $g_j$ is connected to each node $f_i$ with an edge that has weight (or "cost") $C_{i,j}=|f_i-g_j|^2$. This is a maximum matching problem in a complete weighted bipartite graph. It can be solved with the Hungarian algorithm (see also this tutorial on TopCoder).

P.S. If I misunderstood the problem description, this exact same approach works so long as the costs are additive, and then $C_{i,j}$ is the cost of assigning $f_i$ to the point $j$.

  • $\begingroup$ Thank you! Yes, before going to sleep I figured that a good approximation would be to use the plane is very nearly $xy$, and thus minimizing the RMS of the corrected plane is approximately the same as minimizing the residuals of the $z$-components of the initial best fit plus the height corrections, which gives the $\sum_k \big(z_k + f_{\sigma(k)} - \alpha x_k - \beta y_k\big)^2$. The one thing I don't get is that according to the Wiki article on the Hungarian algorithm, it requires the number of tasks and agents to be equal, but in my case, $N>m$. Does that wreck things? $\endgroup$ Commented Oct 21, 2014 at 12:53
  • $\begingroup$ In effect, this is a maximum matching problem in a complete weighted bipartite graph $K_{N,m}$ where $N\neq m$ (I think?). $\endgroup$ Commented Oct 21, 2014 at 13:13
  • $\begingroup$ I was able to code this in a single line of Mathematica code (using FindIndependentEdgeSet), and it worked correctly! Thank you. $\endgroup$ Commented Oct 21, 2014 at 19:35
  • $\begingroup$ @DumpsterDoofus The Mathematica one-linear is certainly a very clean way to do this; thanks for teaching me about that function. $\endgroup$
    – Kirill
    Commented Oct 21, 2014 at 22:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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