Re: N-Dimensional line
- To: mathgroup at smc.vnet.net
- Subject: [mg22815] Re: [mg22785] N-Dimensional line
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Fri, 31 Mar 2000 01:01:18 -0500 (EST)
- References: <200003250858.DAA24215@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
V Stokes wrote:
>
> Does anyone have a reference and/or Mathematica code for finding the
> line
> in N-dimensional space that minimizes the distances from this line to a
> set of
> points in this space?
>
> --Virgil
Your line may be parametrized as V1*t+V2 where V1 is a direction vector,
V2 a point through which the line passes, and t the scalar parameter.
Distance between a point P (regard it as a vector) and this line is the
same as the distance from (P-V2) to the subspace spanned by V1.
Distance from a point to a subspace is found by taking the projection
onto that subspace and subtracting that from the coordinates for the
point. The projection of a point Q onto the (one dimensional) subspace
spanned by V1 is just (Q.V1)*Q/(Q.Q).
With all this we can now put together a function to minimize. First some
auxiliaries.
norm[x_] := Sqrt[x.x]
proj[x_,y_] := (x.y)*y/(y.y)
perp[x_,y_] := x - proj[x,y]
Let's pick a dimension and size of our data set. Also get variables for
the direction vector and point-on-line coordinates.
dim = 2;
npoints = 20;
vec = Array[x,dim];
disp = Array[p,dim];
We want to minimize:
Sum[norm[perp[points[[j]]-disp, vec]], {j,npoints}]
Let's take a concrete example.
SeedRandom[1111];
vv = Table[Random[Real,{-1,1}], {dim}];
dd = Table[Random[Real,{-1,1}], {dim}];
Now we'll make up points that are "near" the line through dd in the
direction of vv.
points = Table[Random[Real,{-10,10}], {npoints}] *
Table[Table[Random[Real,{-.01,.01}],{dim}]+vv,{npoints}] +
Table[Table[Random[Real,{-.01,.01}],{dim}]+dd,{npoints}];
We'll form variable lists to use in FindMinimum.
vars = Join[vec,disp];
varseq = Apply[Sequence,Transpose[{vars,Table[Random[],{2*dim}]}]];
Now we compute values that define our line. Note that these are no
uniquely determined. No matter, any set suffices.
{min,vals} = FindMinimum[
Evaluate[Sum[norm[perp[points[[j]]-disp, vec]], {j,npoints}]],
Evaluate[varseq], MaxIterations->200]
One can check that the result is reasonable in this simple case just by
plotting.
Show[ListPlot[points,DisplayFunction->Identity],
ParametricPlot[(vec/.vals) * t + (disp/.vals), {t,-10,10},
DisplayFunction->Identity],
DisplayFunction->$DisplayFunction]
Daniel Lichtblau
Wolfram Research
- References:
- N-Dimensional line
- From: V Stokes <birdy00@bu.edu>
- N-Dimensional line