Find the rotation matrix to rotate 3-d vectors such that the z-component of the rotated vectors is approximately constant

conic sectionslinear algebraoptimizationpythonrotations

I have an N-body simulation. Each body in the simulation has an array of positions as a function of time. For example, the body Earth has the following positional coordinates (in meters) over 10 years (using time-step of 2 days):

 .. body: Earth

[[ 1.50124878e+11 -8.10072107e+09  0.00000000e+00]
 [ 1.49423093e+11  5.14365190e+09  0.00000000e+00]
 [ 1.49069175e+11  1.02812108e+10  0.00000000e+00]
 ...
 [ 1.49035495e+11 -1.83159842e+10  0.00000000e+00]
 [ 1.49667650e+11 -1.32192204e+10  0.00000000e+00]
 [ 1.50124878e+11 -8.10072107e+09  0.00000000e+00]]

The shape of this array is (1826, 3); that is, 3 position-vector components (x, y, z) taken over 1826 different times. In position-space (for which each scattered point represents the position at a unique time), this looks like:

3d positions example

Since I know the shape is an ellipse in the xy-plane, I can fit the ellipse directly. By fit, I mean find the optimal parameters that minimize some error function (like least squares). For the general conic section in the xy-plane, the formula of the generalized conic (discussed in this post) is: $𝑎𝑥^2+𝑏𝑥𝑦+𝑐𝑦^2+𝑑𝑥+𝑒𝑦+𝑓=0$

But, what if I can find a body for which the z-components of the positional vector are not constant (at f=0)? In this case, fitting the conic section now becomes more difficult (conceptually for me, but also computationally). One solution I've seen briefly mentioned online is to use dimensionality reduction; ie, reduce the data from a 3-d dataset to a 2-d dataset.

I am not sure, but I think that the best way to go about reducing the last dimension of data would be to find the proper rotation matrix such that the z-components of the rotated position vector will be a constant (if the data will allow for it); then, I can use the general conic formula (linked above). If this idea makes sense, then how does one go about finding this rotation matrix? If this idea is nonsense, how does one go about fitting the conic section to 3-d points?

Best Answer

Suppose that we are given (possibly noisy) measurements $v_1,\dots,v_n \in \Bbb R^3$, and we know a priori that these points should trace out an ellipse, so we would like to find the ellipse in $3$-space of best fit.

In order to reduce dimensionality, we could begin by computing the projection of the data onto a plane of best fit. The approach to this via PCA can be described as follows.

  • Begin with the $3 \times n$ matrix $M_1$ whose rows are $v_1,v_2,\dots,v_n$.
  • Subtract the mean from each row. That is, compute the row-vector $\bar v = \frac 1n \sum_{j=1}^n v_j,$ then take$$ M_2 = M_1 - \vec 1v = \pmatrix{v_1 - \bar v\\ \vdots \\ v_n - \bar v}$$
  • Compute a thin SVD of $M$. That is, find $U,\Sigma,V$ such that $M = U \Sigma V^T$. Here, $U$ is a $3 \times n$ matrix whose columns are orthonormal, $\Sigma$ is a diagonal matrix with non-negative diagonal entries $\sigma_1 \geq \sigma_2 \geq \sigma_3 \geq 0$, and $U$ is a $3 \times 3$ orthogonal matrix.
  • The point of the SVD is to extract the characterization of our ellipse. In particular, the first two rows of $U$ (call them $u_1$ and $u_2$) will point in the directions of the major and minor axes (respectively) of our ellipse. The length of the major axis is $\sigma_1$ and the length of the minor axis is $\sigma_2$. We can now parameterize the desired ellipse by $$ p(t) = \bar v + \sigma_1 \cos t\,u_1 + \sigma_2 \sin t\,u_2. $$
  • If you instead want an implicit description of your ellipse, then note that the ellipse will satisfy the following equations: if we take $p = (x,y,z)$, then $$ (p - \bar v) \cdot u_3 = 0, \qquad \frac{[(p - \bar v) \cdot u_1]^2}{\sigma_1^2} + \frac{[(p - \bar v) \cdot u_2]^2}{\sigma_2^2} = 1 $$
Related Question