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?

  • $\begingroup$ An easier to implement method would be to project the data onto a best-fit plane, rather than finding a rotation such that the ellipse is made parallel to the $xy$-plane. One way to do this would be using princpal-component analysis. $\endgroup$ – Omnomnomnom Dec 15 '19 at 12:01
  • $\begingroup$ I have read that SVD is preferred over PCA, but I am not that familiar with either. Will SVD also work for this purpose? Source - stats.stackexchange.com/questions/314046/… $\endgroup$ – allthemikeysaretaken Dec 15 '19 at 12:10
  • $\begingroup$ SVD and EIG are both matrix-computations that are used for PCA. Everything discussed on that post is related to computing the PCA. $\endgroup$ – Omnomnomnom Dec 15 '19 at 12:12

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 $$
  • $\begingroup$ My goal is to use the optimized parameters (semi- major and minor axes $a$ and $b$) to describe the ellipse, estimate period and eccentricity, and compare orbital paths over multiple orbits. Hopefully other things too, just trying to learn more. $\endgroup$ – allthemikeysaretaken Dec 15 '19 at 12:42
  • $\begingroup$ @allthemikeysaretaken If there's a question in your comment that you wanted answered, then I'm not sure what that question is. $\endgroup$ – Omnomnomnom Dec 15 '19 at 12:55
  • $\begingroup$ I was responding to "If you instead want an implicit description of your ellipse". That aside, I am also curious if this method can apply to other conic sections, like a parabola or hyperbola, by the principal axis theorem. Would these apply? $\endgroup$ – allthemikeysaretaken Dec 15 '19 at 12:58
  • 1
    $\begingroup$ The mean is more natural in this setting. Unlike the median, it is invariant under changes of basis. Also, it is only by subtracting the mean that we end up with a matrix whose singular-values correspond to variances. $\endgroup$ – Omnomnomnom Dec 15 '19 at 13:11
  • 3
    $\begingroup$ This answer is generally rock-solid, but I'd like to point out that we need to know more than that the point "trace out" an ellipse -- we need to know that they're distributed somewhat uniformly around the ellipse. As an extreme case, trace out the unit circle uniformly with 50 points; then trace out $-0.1 < \theta < 0.1$ about ten thousand times more. Your algorithm will fit a nice ellipse (of rather high eccentricity) to the short almost vertical segment, and treat the rest of the circle as noise. I mention this because OP's looking at orbits, and data from 1.5 orbits could screw things up. $\endgroup$ – John Hughes Dec 15 '19 at 16:35

I think the idea has some fundamental problems. The first is simple to deal with: you might make the z-coords all equal some constant $k$, and I might make them all equal $-k$ by multiplying each of your results by $$\pmatrix{-1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & -1} $$ so there are at least two equally good solutions. Fine. Restrict to $k > 0$ and you've solved that problem.

The second is more serious: if you find a solution, I can multiply it by $$\pmatrix{c & -s & 0 \\ s & c & 0 \\ 0 & 0 & -1} $$ where $s$ and $c$ are the sine and cosine of any angle, and I've just rotated your ellipse in the plane, so the ellipse formula you derive depends on a choice...that's probably bad.

So ... I think you need to spend a little more time thinking about what possible result you'd like to get from your input.

  • $\begingroup$ This makes sense! That said, my beginner-level understanding of PCA is that it used to find the principal axes, which are the orthogonal axes that are (somehow) determined from the variance between points. That said, is this not just another rotation-transformation? If so, I would expect this same problem to persist when using PCA. Do I understand correctly, or is the expected PCA output supposed to be 2 axes (major and minor)? $\endgroup$ – allthemikeysaretaken Dec 15 '19 at 12:17
  • $\begingroup$ The expected PCA output is a "frame" (3 coordinate unit vectors, $u, v, w$ mutually perpendicular) with the property that when you form a matrix $M = [u,v,w]$, and then compute $Q = BM$, you get a "transformed" collection of data that (for your problem) will lie almost entirely in one coordinate plane, and within that plane, look like an axis-aligned ellipse. (To be honest, you probably need to do PCA on $C$, where $C$ consists of the data in $B$ with the average subtracted out.) If the original data is near-circular, the axis-aligned-ness may be less good than you hoped. $\endgroup$ – John Hughes Dec 15 '19 at 12:29
  • $\begingroup$ Do you think an entirely different approach is better suited to this problem? $\endgroup$ – allthemikeysaretaken Dec 15 '19 at 12:31
  • 1
    $\begingroup$ No, I think PCA is fine. What PCA sould give you is the opportunity to compute the major and minor axes, hence the eccentricity, of an orbit. I'm just saying that if the orbit is circular, then the directions of these axes will be indeterminate, and you should be aware that the results PCA produces should be regarded as "just a pair of random perpendicular vectors in the right plane" rather than having any meaning. $\endgroup$ – John Hughes Dec 15 '19 at 13:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.