3
$\begingroup$

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

$$\epsilon(\sigma)=\sqrt{\Sigma_3}$$

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.

$\endgroup$
  • $\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 Oct 21 '14 at 1:39
  • $\begingroup$ Is there any motivation for why the $f$'s cannot be reused? $\endgroup$ – Kirill Oct 21 '14 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$ – DumpsterDoofus Oct 21 '14 at 12:58
3
$\begingroup$

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$.

$\endgroup$
  • $\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$ – DumpsterDoofus Oct 21 '14 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$ – DumpsterDoofus Oct 21 '14 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$ – DumpsterDoofus Oct 21 '14 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 Oct 21 '14 at 22:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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