Point closest to a set four of lines in 3D

Given four lines in $3D$ (represented as a couple of points), I want to find the point in space which minimizes the sum of distances between this point and every line.

I'm trying to find a way to formulate this as a problem of Least Square, but I'm not quite sure as to how I should. I'm currently trying to use the definition of distance provided at: http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html

Any ideas?

-
I'm trying to avoid using Calculus because I need to use this computation in an algorithm and I don't have access to any stable optimization libraries for C++. – Prateek May 2 '11 at 10:29
@Dennis, is it clear that minimizing the sum of squared distances is the same as minimizing the sum of distances? – lhf May 2 '11 at 11:16
@Dennis: minimizing the sum of the squares is not the same as minimizing the sum of distances. Take three coplanar lines that meet in a point and another line in a plane $d$ away which passes above the intersection point. To minimize the sum of distances you choose the intersection point and find $d$. To minimize the sum of squares, take a point $a$ from the intersection to the other line. The sum of squares is minimized at $a=d/4$ and find $3d^2/4$ compared to $d^2$ for the intersection. – Ross Millikan May 2 '11 at 13:37
That's true, Sorry. – Dennis Gulko May 2 '11 at 17:34
My instinct would be to set up the problem using the formula to which you link in full generality (lots of constants/variables), bash it over the head with something like Mathematica (using calculus), and see if a concise formula for what you want comes out the other side. I don't have Mathematica handy to try it at the moment, though. – Isaac May 4 '11 at 0:22

As pointed out by the others, minimising the sum of distances is different from minimising the sum of squared distances. The former requires the use of an optimisation package. The latter is a least square problem that is fairly well studied with free solvers in abundance. Specifically:

1. Let each straight line, say $L_i$, be specfied by two points $v_i$ and $w_i$ on it.
2. Compute the unit vector $u_i = (v_i-w_i)/\Vert v_i-w_i\Vert$. This is the direction vector of $L_i$. The projection matrix $I-u_iu_i^T$ will project every vector $x\in\mathbb{R}^3$ to the plane $P_i$ that passes through the origin and orthogonal to $L_i$.
3. In particular, the intersection of $L_i$ and $P_i$ is given by $p_i = (I-u_iu_i^T)w_i$ and the distance from a point $x\in\mathbb{R}^3$ to the line $L_i$ is $\Vert (I-u_iu_i^T)x - p_i\Vert$.
4. Hence the sum of squared distances is $\sum_{i=1}^4 \Vert (I-u_iu_i^T)x - p_i\Vert^2$, which can be expanded into $x^T Ax - 2b^T x + c$ with $A = \sum_i(I-u_iu_i^T)$ (a 3x3 matrix), $b = \sum_i p_i$ (a 3-vector) and $c=\sum_i p_i^Tp_i$ (a scalar).
5. Now come the pretty standard stuffs. The minimiser of $x^T Ax - 2b^T x + c$ is the least square solution to $Ax=b$, which in theory is $x = A^+ b$, where $A^+$ denotes the Moore-Penrose pseudoinverse of $A$.
6. In practice, $A$ is usually nonsingular; in this case, $A^+$ is just $A^{-1}$. If you are sure that $A$ is nonsingular, you may simply invert the matrix and calculate $x$ as $A^{-1}b$ or (for better numerical stability) solve the equation by QR factorisation.
7. However, $A$ can at times be singular or nearly singular. If you really want to play safe, you may solve the least square problem by applying singular value decomposition (SVD) or by using the conjugate gradient method. For C++ users, the Eigen3 library already contains a SVD solver. See the section "Least squares solving" in the official tutorial.
-
More specifically: you can use QR decomposition for solving $\mathbf A\mathbf x=\mathbf b$ in point 5, or use SVD (with appropriate thresholding) if it turns out your $\mathbf A$ is rather ill-conditioned. – J. M. Aug 3 '11 at 8:24
@user1551, I am also working on the similar problem. So, I would like to know two things. (1) If we are using the sum of the distances to find the closest point, then do we need to use function minimizes? So, please tell me a simple minimizer that suit for this problem. (2) If we are using least square solutions, then should we the use sum of the square distances? Is this least square solution steps are differ from other least square solutions steps (for e.g. with the steps to be followed in fitting a plane for given points). Do they different or same way? – niro Sep 1 '11 at 20:44
@J. M., I would like to clarify 2 things from you. May be the things what I am asking would be stupid questions, but I want to clear this. Hope you may comment me. What is this threshold and how we can formulate it. Can we get the minimum Eigen value and its respective Eigen vector as the best solution then still do we need this threshold? – niro Sep 1 '11 at 20:47
@g_niro: I have converted your answers to comments. Answers should be reserved for posts that answer the question. But because you do not have 50 reputation points yet, you can only comment on your own questions and answers. So, you didn't do anything wrong; the "add comment" button will only appear for you once you gain 50 points. Here is an explanation of reputation points. – Zev Chonoles Sep 2 '11 at 2:58
@g_niro: (1) To minimise the sum of distances, you need a nonsmooth minimisation package because the distance function is not differentiable. You may try the Nelder-Mead method implemented in the GNU Scientific Library. (2) Linear least square problems arise in various contexts. The interim steps for deriving a solution can be quite different, but up to some point, you always end up with minimising some quadratic form $x^T Ax - 2b^T x + c$ (with $A$ symmetric). Hence the solution forms are similar and same set of computational techniques (QR, SVD, CG, etc.) are applicable. – user1551 Sep 3 '11 at 6:56

I made a program in Mathematica for calculating the point coordinates. The result is a large algebraic formula. I uploaded it to ideone for you.

Here is the program, in case you have Mathematica at hand:

(*Load package*)
Needs["VectorAnalysis"]
(*Define four lines, by specifying 2 points in each one*)
Table[p[i, j] = {x[i, j], y[i, j], z[i, j]}, {i, 4}, {j, 2}];

(*Define the target point*)
p0 = {x0, y0, z0};

(*Define a Norm function // using Std norm squared here*)
norm[a_] := a[[1]]^2 + a[[2]]^2 + a[[3]]^2

(*Define a function for the distance from line i to point v
used http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html (11) *)
d[i_, v_] :=  norm[Cross[(v - p[i, 1]), (v - p[i, 2])]]/norm[p[i, 2] - p[i, 1]]

(*Define a function for the sum of distances*)
dt[p_] := Sum[d[i, p], {i, 4}]

`