# Transforming noisy multivariate data

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

-
Physically, I know there are rotations involved. The data is essentially gravity along x,y,z axes… –  Pam Aug 15 '13 at 21:41
So in general case you may want to adapt this method and with model of ellipsoid you have the transformation –  Kuba Aug 15 '13 at 21:47
Don't forget to upvote any helpful answers...keeps everyone motivated... :) –  cormullion Aug 16 '13 at 10:31
Also (saw this too late to answer) you may want to look at SingularValueDecomposition. –  gpap Aug 16 '13 at 10:33

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]


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]


-
brilliant. The key was nudging it in direction of the offset… –  Pam Aug 16 '13 at 10:34
Still getting used to mathematica, any way of figuring out what the best fit transform (i.e. affine and offset) values are. Would help me tune out my experimental procedure… –  Pam Aug 16 '13 at 11:02
Take a look at the value of sol2; it is a pair containing {minimal value of objective function, rules for values of optimization variables at minimum}. So you can get the value of sxx out by evaluating sxx /. sol2[[2]], and so on. Or simply do affine /. sol2[[2]] and offset /. sol2[[2]]. –  Rahul Narain Aug 16 '13 at 11:08
tnx… that was quick.. ! –  Pam Aug 16 '13 at 11:09

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.

-