3

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?

  • 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. – Ben Grossmann Dec 15 '19 at 12:01
  • 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 - https://stats.stackexchange.com/questions/314046/why-does-andrew-ng-prefer-to-use-svd-and-not-eig-of-covariance-matrix-to-do-pca –  Dec 15 '19 at 12:10
  • SVD and EIG are both matrix-computations that are used for PCA. Everything discussed on that post is related to computing the PCA. – Ben Grossmann Dec 15 '19 at 12:12

2 Answers2

3

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 $$
Ben Grossmann
  • 225,327
  • 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. –  Dec 15 '19 at 12:42
  • @allthemikeysaretaken If there's a question in your comment that you wanted answered, then I'm not sure what that question is. – Ben Grossmann Dec 15 '19 at 12:55
  • 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? –  Dec 15 '19 at 12:58
  • @allthemikeysaretaken In that case you can still use PCA to find the find the plane containing your curve, which will be the plane spanned by $u_1$ and $u_2$, which allows you to reduce the question to a $2D$ problem. That being said, I think that the directions $u_1$ and $u_2$ and the factors $\sigma_1,\sigma_2$ can no longer can be usefully interpreted. – Ben Grossmann Dec 15 '19 at 13:03
  • I am also curious about one other thing. Do you know why the mean is used (instead of the median) before computing the SVD? To my understanding, the median is supposed to be more robust to outliars. –  Dec 15 '19 at 13:08
  • 1
    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. – Ben Grossmann Dec 15 '19 at 13:11
  • 3
    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. – John Hughes Dec 15 '19 at 16:35
  • @JohnHughes that's a fair point. I unconsciously made that assumption based on the data set presented. – Ben Grossmann Dec 15 '19 at 17:22
  • So did I...until I started writing some code, and realized it. :) – John Hughes Dec 15 '19 at 18:46
2

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.

John Hughes
  • 93,729
  • 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)? –  Dec 15 '19 at 12:17
  • 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. – John Hughes Dec 15 '19 at 12:29
  • Do you think an entirely different approach is better suited to this problem? –  Dec 15 '19 at 12:31
  • 1
    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. – John Hughes Dec 15 '19 at 13:08