Transforming noisy multivariate data

Question

I’ve been racking my brain all afternoon on trying to solve this…

I have a set of noisy co-ordinate data

data = {{x1,y1,z1}…..{xn,yn,zn}}

that forms an ellipsoid . I’d like to transform the co-ordinate data into a unit sphere. As an illustration, my original data is in blue, and I’d like to transform it to the red data.

What is the best way to find the transformation coefficients

affine = {{sxx,sxy,sxz}, {syx, syy,syz}, {szx,szy,szz}} 
offset = Table[{ox,oy,oz}, {i,1,Length[data]}]  

such that:

transform = affine.data + offset 

subject to the constraint that

x^2+y^2+z^2 = 1

I assume this requires some kind of least square fit. I can’t figure out how I can get *fit * or NMinimize to make this work. Maybe I am using the wrong functions….

Ideas greatly appreciated….

Thanks, P

Answer 1

The obvious approach of forming a least-squares problem and throwing FindMinimum at it totally works; it just needs a little bit of a nudge in the right direction.

First, let me generate some points distributed on an ellipsoid, with some noise added to make it interesting:

sphere = Normalize /@ RandomVariate[NormalDistribution[], {1000, 3}];
a = RandomVariate[NormalDistribution[], {3, 3}];
b = RandomVariate[NormalDistribution[], 3];
data = a.# + b + 0.1 RandomVariate[NormalDistribution[], 3] & /@ sphere;
ListPointPlot3D[data, BoxRatios -> Automatic]

Here's your transformation, lightly modified:

affine = {{sxx, sxy, sxz}, {syx, syy, syz}, {szx, szy, szz}};
offset = {ox, oy, oz};
transform = affine.# + offset & /@ data;

Now let's just take the criterion $\|\mathbf x_i\|^2=1$, define its total squared error $\sum_i \big(\|\mathbf x_i\|^2-1\big)^2$, and see how well FindMinimum can minimize it:

sol = FindMinimum[Total[(Norm[#]^2 - 1)^2 & /@ transform], Flatten[{affine, offset}]];
ListPointPlot3D[transform /. sol[[2]], BoxRatios -> Automatic]

Oops. What happened? The minimization has shrunk all the points down to a tiny patch on the surface of the unit sphere. Well, that does make all the points have $\|\mathbf x_i\|^2\approx1$, but that's not really what we want. What this shows is that our objective has multiple local minima, and if we want to get the right one we should pick a better starting point. A good first guess for the offset is probably the one that makes the data set centered around zero:

min = Min /@ Transpose[data];
max = Max /@ Transpose[data];
mid = (min + max)/2;

Then we try again knowing a good starting point for the variables, and voilà:

sol2 = FindMinimum[Total[(Norm[#]^2 - 1)^2 & /@ transform], 
   Transpose[{Flatten[{affine, offset}], Flatten[{IdentityMatrix[3], -mid}]}]];
ListPointPlot3D[transform /. sol2[[2]], BoxRatios -> Automatic]

Answer 2

I'm not completely sure if this is the type of transform that you are looking for, but the following might help at least with fitting a set of data to an ellipsoid. First, generate some data. I'm using the parametric equations for an ellipsoid just because I just learned them ;-).

data = With[{a = 1, b = 2, c = 1.5}, 
    Table[{a Cos[u] Sin[v], b Sin[u] Sin[v], c Cos[v]}, {u, 0, 2 Pi, 
      2 Pi/20}, {v, -2 Pi, Pi, Pi/30}]]~Flatten~{1, 2};
ListPointPlot3D[data, SphericalRegion -> True, PlotStyle -> Blue]

Fitting the data to the cartesian coordinate equation for an ellipsoid is straightforward, but we first need to append the value of the function (which in this case is 1) to each set of points. I do this to the data within the nonlinearmodelfit function:

nlm = NonlinearModelFit[
   Partition[Flatten@Riffle[data, 1], 4],
   x^2/a^2 + y^2/b^2 + z^2/c^2,
   {a, b, c},
   {x, y, z}];
nlm["BestFitParameters"]
(* {a -> 1., b -> 2., c -> 1.5} *)

Here's were I don't know what you are looking for. I presume that if you want to transform the data, you want to normalize the three coefficients a, b, c such that they are equal. (Note the first line of code below is very slow.)

newdata = # /({a, b, c} /. nlm["BestFitParameters"]) & /@ data;
nlm2 = NonlinearModelFit[
   Partition[Flatten@Riffle[newdata, 1], 4],
   x^2/a^2 + y^2/b^2 + z^2/c^2,
   {a, b, c},
   {x, y, z}];
nlm2["BestFitParameters"]
(* {a -> 1., b -> 1., c -> 1.} *)
ListPointPlot3D[{data, newdata}, SphericalRegion -> True, 
 PlotStyle -> {Blue, Red}]

I haven't played with the AspectRatio in this image, so my sphere (red) looks squished.