Category Archives: Mathematics

Properties of ellipses by matrix coefficients – III – coordinates of points with extremal radii

Posted on by

This post requires Javascript to display formulas!

A centered, rotated ellipse can be defined by matrices which operate on position-vectors for points on the ellipse. The topic of this post series is the relation of the coefficients of such matrices to some basic geometrical properties of an ellipse. In the previous posts

Properties of ellipses by matrix coefficients – I – Two defining matrices
Properties of ellipses by matrix coefficients – II – coordinates of points with extremal y-values

we have found that we can use (at least) two matrix based approaches:

  • One reflects a combination of two affine operations applied to a unit cycle. This approach lead us to a non-symmetric matrix, which we called AE. Its coefficients ((a, b), (c, d)) depend on the lengths of the ellipses’ principal axes and trigonometric functions of its rotation angle.
  • The second approach is based on coefficients of a quadratic form which describes an ellipse as a special type of a conic section. We got a symmetric matrix, which we called Aq.

We have shown how the coefficients α, β, γ of Aq and a further coefficient δ of the quadratic form can be expressed in terms of the coefficients of AE. Furthermore, we have derived equations for the lengths σ1, σ2 of the ellipse’s principal axes and the rotation angle by which the major axis is rotated against the x-axis of the Euclidean coordinate system [ECS] we work with. We have also found equations for the components of the position vectors to those points of the ellipse with maximum y-values. A major result was that the eigenvalues and eigenvectors of Aq completely control the ellipse’s properties.

In this post we determine the components of the vectors to the end-points of the ellipse’s principal axes in terms of the coefficients of Aq. Afterward we shall test our formulas by a Python program and plots for a specific example.

Reduced matrix equation for an ellipse

Our centered, but rotated ellipse is defined by a quadratic form, i.e. by a polynomial equation with quadratic terms in the components xe and ye of position vectors to points on the ellipse:

\[
\alpha\,x_e^2 \, + \, \beta \, x_e y_e \, + \, \gamma \, y_e^2 \:=\: \delta
\]

The quadratic polynomial can be formulated as a matrix operation applied to position vectors ve = (xe, ye)T. With the the quadratic and symmetric matrix Aq

\[ \pmb{\operatorname{A}}_q \:=\:
\begin{pmatrix} \alpha & \beta / 2 \\ \beta / 2 & \gamma \end{pmatrix}
\]

we can rewrite the polynomial equation for the centered ellipse as

\[
\pmb{v}_e^T \circ \pmb{\operatorname{A}}_q \circ \pmb{v}_e^T \:=\: \delta \,=\, \sigma_1^2\, \sigma_2^2, \quad \operatorname{with}\: \pmb{v_e} \,=\, \begin{pmatrix} x_e \\ y_e \end{pmatrix}.
\]

Just to cover another relation you may find in books: We could have included the δ-term in a somewhat artificial (3×3)-matrix

\[
\pmb{\operatorname{A}}_q^e \:=\:
\begin{pmatrix} \alpha & \beta / 2 & 0 \\ \beta / 2 & \gamma & 0 \\ 0 & 0 & – \delta \end{pmatrix},
\]

and applied this matrix to an artificially extended position vector to reproduce our definition equation:

\[
\pmb{\operatorname{A}}_q^e \, \circ \, \begin{pmatrix} x_e \\ y_e \\ 1 \end{pmatrix} \:=\: 0
\]

However, this formal aspect will not help much to solve the equations coming below. We want to describe the vectors to the principal axes’ end-points by mathematical expressions that depend on α, β, γ, δ – and both matrices will of course deliver the same results. But: There are still two different approaches to achieve our objective.

Method 1 to determine the vectors to the principal axes’ end points

My readers have certainly noticed that we have already gathered all required information to solve our task. In the first post of this series we have performed an eigendecomposition of our symmetric matrix Aq. We found that the two eigenvectors of Aq for respective eigenvalues λ1 and λ2 point along the principal axes of our rotated ellipse:

\[ \begin{align}
\lambda_1 \: &: \quad \pmb{\xi_1} \:=\: \left(\, {1 \over \beta} \left( (\alpha \,-\, \gamma) \,-\, \left[\, \beta^2 \,+\, \left(\gamma \,-\, \alpha \right)^2\,\right]^{1/2} \right), \: 1 \, \right)^T \\
\lambda_2 \: &: \quad \pmb{\xi_2} \:=\: \left(\, {1 \over \beta} \left( (\alpha \,-\, \gamma) \,+\, \left[\, \beta^2 \,+\, \left(\gamma \,-\, \alpha \right)^2\,\right]^{1/2} \right), \: 1 \, \right)^T
\end{align}
\]

The T symbolizes a transposition operation. The eigenvalues are related to the Aq-coefficients by the following equations:

\[ \begin{align}
\lambda_1 \:&=\: {1 \over 2} \left(\, \left( \alpha \,+\, \gamma \right) \,+\, \left[ \beta^2 \,+\, \left(\gamma \,-\, \alpha \right)^2 \,\right]^{1/2} \,\right) \\
\lambda_2 \:&=\: {1 \over 2} \left(\, \left( \alpha \,+\, \gamma \right) \,-\, \left[ \beta^2 \,+\, \left(\gamma \,-\, \alpha \right)^2 \,\right]^{1/2} \,\right)
\end{align}
\]

These eigenvalues correspond to the squares of the lengths of the ellipse’s axes.

\[ \begin{align}
\lambda_1 \:&=\: \sigma_1^2 \\
\lambda_2 \:&=\: \sigma_2^2 \\
\end{align}
\]

Therefore, we can simply take the components of the normalized vectors

\[ \begin{align}
\lambda_1 \: &: \quad \pmb{\xi_1^n} \:=\: {1 \over \|\pmb{\xi_1}\|}\, \pmb{\xi_1} \\
\lambda_2 \: &: \quad \pmb{\xi_2^n} \:=\: {1 \over \|\pmb{\xi_2}\|}\, \pmb{\xi_2}
\end{align}
\]

and multiply them with the square-root of the respective eigenvalues to get the vector components to the end-points of the ellipse’s axes:

\[ \begin{align}
\pmb{\xi_1}^{rmax} \:&=\: \sqrt{\lambda_1} * \pmb{\xi_1^n} \\
\pmb{\xi_2}^{rmax} \:&=\: \sqrt{\lambda_2} * \pmb{\xi_1^n}
\end{align}
\]

This is trivial regarding the algebraic operations, but results in lengthy (and boring) expressions in terms of the matrix coefficients. So, I skip to write down all the terms. (We do not need it for setting up ordered numerical programs.)

Remember that you could in addition replace (α, β, γ, δ) by coefficients (a, b, c, d) of matrix AE. See the first post of this series for the formulas. This would, however, produce even longer equation terms.

Equation for points with maximum radius values

We define again some convenience variables:

\[ \begin{align}
a_h \,& =\, {\alpha \over \gamma} \\
b_h \,& =\, {1 \over 2 } {\beta \over \gamma} \\
d_h \,& =\, {\delta \over \gamma} \\
g_h \,& =\, a_h \,-\, b_h^2 \\
f_h \,& =\, 1 \,+\, b_h^2 \,-\, g_h \phantom{\huge{(}}
\end{align}
\]

and

\[
\xi_h \,=\, { \left[\, 4\,d_h\, g_h\, b_h^2 \,+\, d_h\, f_h^2 \,\right] \over 2\, \left[\, 4\,b_h^2\,g_h^2 \,+\, g_h\,f_h^2\,\right] } \phantom{\Huge{(}} \\
\]
\[
\eta_h \,=\, { b_h^2 \, d_h^2 \over \left[\, 4\,b_h^2\,g_h^2 \,+\, g_h\,f_h^2\,\right] } \phantom{\Huge{(}}
\]

We find for ye:

\[ \begin{align}
y_e \:&=\: – b_h \, x_e \, \pm \, \left[\,d_h \,-\, \left( a_h \,-\, b_h^2 \right)\, x_e^2 \, \right]^{1/2} \\
\:&=\: – b_h \, x_e \, \pm \, \left[\,d_h \,-\, g_h\, x_e^2 \, \right]^{1/2} \phantom{\huge{(}}
\end {align}
\]

We pick the ye with the positive term in the following steps. (The way for the solution with the negative term in ye is analogous.) The square of ye is:

\[
y_e^2 \:=\: – d_h \,+\, \left(b^2 \,-\, g\right)\, x_e^2 \,-\, 2 \, b_h \, x_e \,\left[d \,-\, g\,x_e^2 \right]^{1/2}
\]

To find an extremal value of the radius we differentiate and set the derivative to zero:

\[
{\partial \, \left(y_e^2 \,+\, x_e^2\right) \over \partial \, x_e} \:=\: 0 \: \Rightarrow
\]
\[ \begin{align}
& \left(\, 1\,+\, b_h^2 \,-\, g_h\,\right) \, x_e \,-\, b_h \, \left[\, d_h \,-\, g_h\, x_e^2 \, \right]^{1/2} \\
&+\, b_h\,g_h\,x_e^2 \, {1 \over \left[\, d_h \,-\, g_h \, x_e^2 \, \right]^{1/2} } \:=\: 0
\end{align}
\]

This results in

\[
f_h\, x_e \, \left[ d_h \,-\, g_h\, x_e^2\right]^{1/2} \:=\: b_h\,d_h \,-\, 2\, b_h\,g_h\,x_e^2 \, .
\]

Solution for xe-values of the end-points of the principal axes

We take the square of both sides and reorder terms to get

\[
\left[\, 4\,b_h^2\,g_h^2 \,+\, g_h\,f_h^2\,\right]\, x_e^4 \,-\, \left[\, 4\,d_h\, g_h\, b_h^2 \,+\, d_h\, f_h^2 \,\right]\, x_e^2 \,+\, b_h^2 \, d_h^2 \:=\: 0
\]
\[
x_e^4 \,-\, { \left[\, 4\,d_h\, g_h\, b_h^2 \,+\, d_h\, f_h^2 \,\right] \over \left[\, 4\,b_h^2\,g_h^2 \,+\, g_h\,f_h^2\,\right] } \, x_e^2 \:=\:
-\, { b_h^2 \, d_h^2 \over \left[\, 4\,b_h^2\,g_h^2 \,+\, g_h\,f_h^2\,\right] }
\]

With

\[ \begin{align}
\xi_h \,&=\, { \left[\, 4\,d_h\, g_h\, b_h^2 \,+\, d_h\, f_h^2 \,\right] \over 2\, \left[\, 4\,b_h^2\,g_h^2 \,+\, g_h\,f_h^2\,\right] } \\
\eta_h \,&=\, { b_h^2 \, d_h^2 \over \left[\, 4\,b_h^2\,g_h^2 \,+\, g_h\,f_h^2\,\right] } \phantom{\Huge{)^A}}
\end{align}
\]

we have

\[
x_e^4 \, -\, 2 \,\xi_h \,x_e^2 \:=\: – \eta_h
\]

With the help of a quadratic supplement we get

\[
\left[\, x_e^2 \, -\, \xi_h \right]^2 \:=\: \xi_h^2 \:-\: \eta_h
\]

and find the solution

\[
x_e \:=\: \pm \, \sqrt{ \xi_h \,\pm\, \sqrt{\, \xi_h^2 \,-\,\eta_h \,} }
\]

A detailed analysis also for the other ye-expression (see above) leads to further solutions for the coordinates (=vector component values) of points with extremal values for the radii. These are the end-points of the principal axes of the ellipse:

\[ \begin{align}
x_{e1}^{rmax} \:&=\: -\, \sqrt{ \, \xi_h \,-\, \sqrt{\, \xi_h^2 \,-\,\eta \,} } \\
y_{e1}^{rmax} \:&=\: +\, \sqrt{ \, d_h \,-\, g_h \, \left(x_{e1}^{rmax}\right)^2 \,} \,-\, b_h \, x_{e1}^{rmax} \\
x_{e2}^{rmax} \:&=\: -\, x_{e1}^{rmax} \phantom{\huge{)}} \\
y_{e2}^{rmax} \:&=\: -\, y_{e1}^{rmax} \\
x_{e3}^{rmax} \:&=\: +\, \sqrt{ \, \xi_h \,+\, \sqrt{\, \xi_h^2 \,-\,\eta \,} } \phantom{\huge{)}} \\
y_{e3}^{rmax} \:&=\: +\, \sqrt{ \, d_h \,-\, g_h \, \left(x_{e3}^{rmax}\right)^2 \,} \,-\, b_h \, x_{e3}^{rmax} \\
x_{e4}^{rmax} \:&=\: -\, x_{e3}^{rmax} \phantom{\huge{)}} \\
y_{e4}^{rmax} \:&=\: -\, y_{e3}^{rmax}
\end {align}
\]

I leave it to the reader to expand the convenience variables into terms containing the original coefficients α, β, γ, δ.

Plots

It is easy to write a Python program, which calculates and plots the data of an ellipse and the special points with extremal values of the radii and extremal values of ye. The general steps which I followed were:

Step 0: Create 100 points a unit circle. Save the coordinates in Python lists (or Numpy arrays). Use Matplotlib’s plot(x,y)-function to plot the vectors.

Step 1: Create an axis-parallel ellipse with values for the axes ha = 2.0 and hb = 1.0 along the x- and the y-axis of the Euclidean coordinate system [ECS]. Do this by applying a diagonal scaling matrix Dσ1, σ2 (see the first post of this series).

Step 2: Rotate the ellipse bei π/3 (60 °). Do this by applying a rotation matrix Rπ/3 to the position vectors of your ellipse (with the help of Numpy). Alternatively, you can first create the matrices, perform a matrix multiplication and then apply the resulting matrix to the position vectors of your unit circle.

(The limiting lines have been calculated by the formulas given above.)

Step 3: Determine the coefficients of combined matrix AE = Rπ/3 ○ Dσ1, σ2

I got for the coefficients ( (a, b), (c, d) ) of AE : 

A_ell = 
[[ 1.         -0.8660254 ]
[ 1.73205081  0.5       ]]

Step 3: Determine the coefficients of the matrix Aq by the formulas given in the first post of this series. I got 

A_q = 
[[ 3.25       -1.29903811]
[-1.29903811  1.75      ]]

For δ I got: 

delta =  4.0

which is consistent with the length-values of the principal axes.

Step 4: Determine values for the eigenvalues λ1 and λ2 from the Aq-coefficients by the formulas given in the first post. Also calculate them by using Numpy’s
eigenvalues, eigenvectors = numpy.linalg.eig(A_q). Theory tells us that these values should be exactly λ1 = 4 and λ1 = 1. I got 

Eigenvalues from A_q:  lambda_1 = 4. :: lambda_2 = 1.

Step 5: Determine the components of the normalized eigenvectors with the help of numpy.linalg.eig(A_q). I got:

Components of normalized eigenvectors by theoretical formulas from A_q coefficients: 
ev_1_n :  -0.8660254037844386  :  0.5000000000000002
ev_2_n :  0.5000000000000001  :  0.8660254037844385

Eigenvectors from A_q via numpyy.linalg.eig():  
ev_1_num :  0.8660254037844387  :  -0.5000000000000001
ev_2_num :  0.5000000000000001  :  0.8660254037844387

The deviation between ev_1_n and ev_1_num is just due to a difference by -1. This is correct as the eigenvectors are unique only up to a minus-sign in all components.

Step 6: Calculate the sinus of the rotation angle of our ellipse from Aq– and Aq-coefficients. The theoretical value is sin(2 π/3) = sin(2 pi/3) = 0.8660254037844387. I got: 

sin(2. * rotation angle) of major axis of the ellipse against the ECS x-axis from A_E coefficients: 
sin_2phi-A_E  =  0.8660254037844388

sin(2. * rotation angle) of major axis of the ellipse against the ECS x-axis from from eigenvectors of A_q:
sin_2phi-ev_A_q =  0.8660254037844387 

sin(2. * rotation angle) of major axis of the ellipse against the ECS x-axis from A_q-coefficients:
sin_2phi-coeff-A_q =  0.8660254037844388

Perfect!

Step 7: Plot the end-points of the normalized eigenvectors of Aq:

Note that in our example case the end-point of the eigenvector along the minor axis must be located exactly on the elliptic curve as the ellipses minor axes has a length of b=1!

Step 8: Calculate the components of the vectors to data-points of the ellipse with maximal absolute ye-values from the Aq-coefficients given in the previous post. Plot these data-points (here in green color).

Step 9: Calculate the components of the vectors to data-points of the ellipse with maximal values of the radii with the help of the complex formulas presented in this post and plot these points in addition.

Conclusion

In this mini-series of posts we have performed some small mathematical exercises with respect to centered and rotated ellipses. We have calculated basic geometrical properties of such ellipses from the coefficients of matrices which define ellipses in algebraic form. Linear Algebra helped us to understand that the eigenvectors and eigenvalues of a symmetric matrix, whose coefficients stem from a quadratic equation (for a conic section), control both the orientation and the lengths of the ellipse’s axes completely.

This knowledge is useful in some Machine Learning [ML] context where elliptic data appear as projections of multivariate normal distributions. Multivariate Gaussian probability functions control properties of a lot of natural objects. Experience shows that certain types of neural networks may transform such data into multivariate normal distributions in latent spaces. An evaluation of the numerical data coming from such ML-experiments often delivers the coefficients of defining matrices for ellipses.

In my blog I now return to the study of with shearing operations applied to circles, spheres, ellipses and 3-dimensional ellipsoids. Later I will continue with the study of multivariate normal distributions in latent spaces of Autoencoders. For both of these topics the knowledge we have gathered regarding the matrices behind ellipses will help us a lot.

 

Properties of ellipses by matrix coefficients – II – coordinates of points with extremal y-values

Posted on by

This post requires Javascript to display formulas!

In the previous post of this mini-series

Properties of ellipses by matrix coefficients – I – Two defining matrices

I have discussed how the coefficients of two matrices which can be used to define a centered, rotated ellipse can be used to calculate geometrical properties of the ellipse:

The lengths σ1, σ2 of the ellipse’s principal axes and the rotation angle by which the major axis is rotated against the x-axis of the Euclidean coordinate system [ECS] we work with.

But there are other properties which are interesting, too. A centered, rotated ellipse has two points with extremal values in their y-coordinates. Can we express the coordinates – or equivalently the components of respective position vectors – in terms of the basic matrix coefficients?

The answer is, of course, yes. This post provides a derivation of respective formulas.

Matrix equation for an ellipse

In the last post we have shown that a centered ellipse is defined by a quadratic form, i.e. by a polynomial equation with quadratic terms in the components xe and ye of position vectors for points of the ellipse:

\[
\alpha\,x_e^2 \, + \, \beta \, x_e y_e \, + \, \gamma \, y_e^2 \:=\: \delta
\]

The quadratic polynomial can be formulated as a matrix operation applied to position vectors ve of points on an ellipse. With the the quadratic and symmetric matrix Aq

\[ \pmb{\operatorname{A}}_q \:=\:
\begin{pmatrix} \alpha & \beta / 2 \\ \beta / 2 & \gamma \end{pmatrix}
\]

we can rewrite the polynomial equation for the centered ellipse as

\[
\pmb{v}_e^T \circ \pmb{\operatorname{A}}_q \circ \pmb{v}_e^T \:=\: \delta \,=\, \sigma_1^2\, \sigma_2^2, \quad \operatorname{with}\: \pmb{v_e} \,=\, \begin{pmatrix} x_e \\ y_e \end{pmatrix}.
\]

We could have artificially included the δ-term in a (3×3)-matrix, but this formal aspect will not help much to solve the equations coming below.

We also know already:

\[
\operatorname{det}\left( \pmb{\operatorname{A}}_q \right) \:=\: \left(\alpha\gamma \,-\, {1\over 4}\beta^2\right) \, \ge \, 0
\]

Aq is an invertible matrix (as was to be expected).

Points of the ellipse with extremal ye-values

Our first step to get the xe and ye coordinates for extremal points is to evaluate the quadratic form with respect to ye. We define:

\[ \begin{align}
a_h \,& =\, {\alpha \over \gamma} \\
b_h \,& =\, {1 \over 2 } {\beta \over \gamma} \\
d_h \,& =\, {\delta \over \gamma}
\end{align}
\]

We reorder terms of the equation and add a supplemental quadratic term on both sides:

\[
y_e^2 \, +\, b_h\,x_e\,y_e \, +\, b_h^2 \, x_e^2 \:=\: d_h \,-\, a_h \, x_e^2 \,+\, b_h^2 \, x_e^2
\]

We evaluate the complete quadratic term ob the left side to get

\[
y_e \:=\: – b_h \, x_e \, \pm \, \left[\,d_h \,-\, \left( a_h \,-\, b_h^2 \right)\, x_e^2 \, \right]^{1/2}
\]

Let us first focus on the positive of the two alternative terms:

\[
y_e \:=\: – b_h \, x_e \,+\, \left[\,d_h \,-\, \left( a_h \,-\, b_h^2 \right)\, x_e^2 \, \right]^{1/2}
\]

We get an extremal ye by evaluating the derivative with respect to xe

\[
{ \partial \, y_e \over \partial \, x_e } \:=\: 0
\]

This means

\[
– b_h \, – \, { \left(a_h \,-\, b_h^2 \right)\,x_e \over \left[\,d_h \,-\, \left( a_h \,-\, b_h^2 \right)\, x_e^2 \, \right]^{1/2} } \:=\: 0
\]

Getting rid of the denominator gives

\[
– b_h \, \left[\,d_h \,-\, \left( a_h \,-\, b_h^2 \right)\, x_e^2 \, \right]^{1/2} \, – \, \left(a_h \,-\, b_h^2 \right)\,x_e \:=\: 0
\]

and

\[
x_e \:=\: { – b_h \, \left[\,d_h \,-\, \left( a_h \,-\, b_h^2 \right)\, x_e^2 \, \right]^{1/2} \over \left(a_h \,-\, b_h^2 \right) }
\]

Squaring both sides and reordering leads to

\[
x_e^2 \:=\: { – b_h^2 \, d_h \over a_h \, \left(a_h \,-\, b_h^2 \right) }.
\]

Using the old coefficients provides the final result for xe:

\[
x_e \:=\: \pm \, {1 \over 2} \, \beta \, \left[ { \delta \over \alpha \,\left( \alpha \, \gamma \,-\, {1 \over 4} \beta^2 \right) } \right]^{1/2}
\]

We now follow the alternative solution for ye (see above). After a calculation of the ye-values from the derived xe values, we get the components for the two position vectors to the points with extremal y-values on the ellipse:

\[ \begin{align}
x_{e1}^{max} \:&=\: – \, {1 \over 2} \, \beta \, \left[ { \delta \over \alpha \,\left( \alpha \, \gamma \,-\, {1 \over 4} \beta^2 \right) } \right]^{1/2} \\
y_{e1}^{max} \:&=\: – \, {1 \over 2} \, {\beta \over \gamma} \, x_{e1}^{max} \,+\,
\left[ {\delta \over \gamma} \,-\, \left( {\alpha \over \gamma} \,-\, {1 \over 4} {\beta^2 \over \gamma^2} \right) \, x_{e1}^{max} \right]^{1/2}
\end{align}
\]

and

\[ \begin{align}
x_{e2}^{max} \:&=\: + \, {1 \over 2} \, \beta \, \left[ { \delta \over \alpha \,\left( \alpha \, \gamma \,-\, {1 \over 4} \beta^2 \right) } \right]^{1/2} \\
y_{e2}^{max} \:&=\: + \, {1 \over 2} \, {\beta \over \gamma} \, x_{e2}^{max} \,-\,
\left[ {\delta \over \gamma} \,-\, \left( {\alpha \over \gamma} \,-\, {1 \over 4} {\beta^2 \over \gamma^2} \right) \, x_{e2}^{max} \right]^{1/2}
\end{align}
\]

Coefficients of the matrix Ae

In my previous post I have discussed yet another matrix AE which also can be used to define an ellipse. This matrix summarizes two affine transformations of a centered unit circle: AE = Rφ ○ Dσ1, σ2. Dσ1, σ2.

You find the relations between the coefficients (a, b, c, d) of matrix AE and the coefficients (α, β, γ) of matrix Aq and δ in my previous post. This will allow you to calculate the vectors to the extremal points of an ellipse in terms of the coefficients (a, b, c, d).

Conclusion

In this post we have again used the coefficients of a matrix which defines an ellipse via a quadratic form to get information about a geometrical property.
We can now calculatethe components of the position vectors to the two points of an ellipse with extremal y-values as functions of the matrix coefficients.

In the next post

Properties of ellipses by matrix coefficients – III – coordinates of points with extremal radii

I will show you how to calculate the components of the end points of the principal axes of the ellipse with the help of our matrix for a quadratic form. I will also use our theoretical results for plots of some ellipses’ axes and of their extremal points. We will also compare theoretical predictions with numerically evaluated values.

 

Properties of ellipses by matrix coefficients – I – Two defining matrices

Posted on by

This post requires Javascript to display formulas!

For my two current post series on multivariate normal distributions [MNDs] and on shear operations

Multivariate Normal Distributions – II – random vectors and their covariance matrix
Fun with shear operations and SVD

some geometrical and algebraic properties of ellipses are of interest. Geometrically, we think of an ellipse in terms of their perpendicular two principal axes, focal points, ellipticity and rotation angles with respect to a coordinate system. As elliptic data appear in many contexts of physics, chaos theory, engineering, optics … ellipses are well studied mathematical objects. So, why a post about ellipses in the Machine Learning section of a blog?

In my present working context ellipses appear as a side result of statistical multivariate normal distributions [MNDs]. The projections of multidimensional contour hyper-surfaces of a MND within the ℝn onto coordinate planes of an Euclidean Coordinate System [ECS] result in 2-dimensional ellipses. These ellipses are typically rotated against the axes of the ECS – and their rotation angles reflect data correlations. The general relations of statistical vector data with projections of multidimensional MNDs is somewhat intricate.

Data produced in numerical experiments, e.g. in a Machine Learning context, most often do not give you the geometrical properties of ellipses, which some theory may have predicted, directly. Instead you may get averaged values of statistical vector distributions which correspond to a kind of algebraic coefficients. These coefficients can often be regarded as elements of a matrix. The mathematical reason is that ellipses can in general be defined by matrices operating on position vectors. In particular: Coefficients of quadratic polynomial expressions used to describe ellipses as conic sections correspond to the coefficients of a matrix operating on position vectors.

So, when I became confronted with multidimensional MNDs and their projections on coordinate planes the following questions became interesting:

  • How can one derive the lengths σ1, σ2 of the perpendicular principal axes of an ellipse from data for the coefficients of a matrix which defines the ellipse by a polynomial expression?
  • By which formula do the matrix coefficients provide the inclination angle of the ellipse’s primary axes with the x-axis of a chosen coordinate system?

You may have to dig a bit to find correct and reproducible answers in your math books. Regarding the resulting mathematical expressions I have had some bad experiences with ChatGPT. But as a former physicist I take the above questions as a welcome exercise in solving quadratic equations and doing some linear algebra. So, for those of my readers who are a bit interested in elementary math I want to answer the posed questions step by step and indicate how one can derive the respective formulas. The level is moderate – you need some knowledge in trigonometry and/or linear algebra.

Centered ellipses and two related matrices

Below I regard ellipses whose centers coincide with the origin of a chosen ECS. For our present purpose we thus get rid of some boring linear terms in the equations we have to solve. We do not loose much of general validity by this step: Results of an off-center ellipse follow from applying a simple translation operation to the resulting vector data. But I admit: Your (statistical) data must give you some information about the center of your ellipse. We assume that this is the case.

Our ellipses can, however, be rotated with respect to a chosen ECS. I.e., their longer principal axes may be inclined by some angle φ towards the x-axis of our ECS.

There are actually two different ways to define a centered ellipse by a matrix:

  • Alternative 1: We define the (rotated) ellipse by a matrix AE which results from the (matrix) product of two simpler matrices: AE = Rφ ○ Dσ1, σ2. Dσ1, σ2 corresponds to a scaling operation applied to position vectors for points located on a centered unit circle. Aφ describes a subsequent rotation. AE summarizes these geometrical operations in a compact form.
  • Alternative 2: We define the (rotated) ellipse by a matrix Aq which combines the x- and y-elements of position vectors in a polynomial equation with quadratic terms in the components (see below). The matrix defines a so called quadratic form. Geometrically interpreted, a quadratic form describes an ellipse as a special case of a conic section. The coefficients of the polynomial and the matrix must, of course, fulfill some particular properties.

While it is relatively simple to derive the matrix elements from known values for σ1, σ2 and φ it is a bit harder to derive the ellipse’s properties from the elements of either of the two defining matrices. I will cover both matrices in this post. For many practical purposes the derivation of central elliptic properties from given elements of Aq is more relevant and thus of special interest in the following discussion.

Matrix AE of a centered and rotated ellipse: Scaling of a unit circle followed by a rotation

Our starting point is a unit circle C whose center coincides with our ECS’s origin. The components of vectors vc to points on the circle C fulfill the following conditions:

\[
\pmb{C} \::\: \left\{ \pmb{v}_c \:=\: \begin{pmatrix} x_c \\ y_c \end{pmatrix} \:=\: \begin{pmatrix} \operatorname{cos}(\psi) \\ \operatorname{sin}(\psi) \end{pmatrix}, \quad 0\,\leq\,\psi\, \le 2\pi \right\}
\]

and

\[
x_c^2 \, +\, y_c^2 \,=\; 1
\]

We define an ellipse E(σ1, σ2) by the application of two linear operations to the vectors of the unit circle:

\[
\pmb{E}_{\sigma_1, \, \sigma_2} \::\: \left\{ \, \pmb{v}_e \:=\: \begin{pmatrix} x_e \\ y_e \end{pmatrix} \:=\: \pmb{\operatorname{R}}_{\phi} \circ \pmb{\operatorname{D}}_E \circ \pmb{v_c} , \quad \pmb{v_c} \in \pmb{C} \, \right\}
\]

DE is a diagonal matrix which describes a stretching of the circle along the ECS-axes, and Rφ is an orthogonal rotation matrix. The stretching (or scaling) of the vector-components is done by

\[
\pmb{\operatorname{D}}_E \:=\: \begin{pmatrix} \sigma_1 & 0 \\ 0 & \sigma_2 \end{pmatrix},
\]
\[
\pmb{\operatorname{D}}_E^{-1} \:=\: \begin{pmatrix} {1 \over \sigma_1} & 0 \\ 0 & {1 \over \sigma_2} \end{pmatrix},
\]

The coefficients σ1, σ2 obviously define the lengths of the principal axes of the yet unrotated ellipse. To be more precise: σ1 is half of the diameter in x-direction, σ1 is half of the diameter in y-direction.

The subsequent rotation by an angle φ against the x-axis of the ECS is done by

\[
\pmb{\operatorname{R}}_{\phi} \:=\:
\begin{pmatrix} \operatorname{cos}(\phi) & – \,\operatorname{sin}(\phi) \\ \operatorname{sin}(\phi) & \operatorname{cos}(\phi)\end{pmatrix}
\:=\: \begin{pmatrix} u_1 & -\,u_2 \\ u_2 & u_1 \end{pmatrix}
\]
\[
\pmb{\operatorname{R}}_{\phi}^T \:=\: \pmb{\operatorname{R}}_{\phi}^{-1} \:=\: \pmb{\operatorname{R}}_{-\,\phi}
\]

The combined linear transformation results in a matrix AE with coefficients ((a, b), (c, d)):

\[ \begin{align}
\pmb{\operatorname{A}}_E \:=\: \pmb{\operatorname{R}}_{\phi} \circ \pmb{\operatorname{D}}_E \:=\:
\begin{pmatrix} \sigma_1\,u_1 & -\,\sigma_2\,u_2 \\ \sigma_1\,u_2 & \sigma_2\,u_1 \end{pmatrix} \:=\:: \begin{pmatrix} a & b \\ c & d \end{pmatrix}
\end{align}
\]

These is the first set of matrix coefficients we are interested in.

Note:

\[ \pmb{\operatorname{A}}_E^{-1} \:=\:
\begin{pmatrix} {1 \over \sigma_1} \,u_1 & {1 \over \sigma_1}\,u_2 \\ -{1 \over \sigma_2}\,u_2 & {1 \over \sigma_2}\,u_1 \end{pmatrix}
\]
\[
\pmb{v}_e \:=\: \begin{pmatrix} x_e \\ y_e \end{pmatrix} \:=\: \pmb{\operatorname{A}}_E \circ \begin{pmatrix} x_c \\ y_c \end{pmatrix}
\]
\[
\pmb{v}_k \:=\: \begin{pmatrix} x_c \\ y_c \end{pmatrix} \:=\: \pmb{\operatorname{A}}_E^{-1} \circ \begin{pmatrix} x_e \\ y_e \end{pmatrix}
\]

We use

\[ \begin{align}
u_1 \,&=\, \operatorname{cos}(\phi),\quad u_2 \,=\,\operatorname{sin}(\phi), \quad u_1^2 \,+\, u_2^2 \,=\, 1 \\
\lambda_1 \,&: =\, \sigma_1^2, \quad\quad \lambda_2 \,: =\, \sigma_2^2
\end{align}
\]

and find

\[ \begin{align}
a \,&=\, \sigma_1\,u_1, \quad b \,=\, -\, \sigma_2\,u_2, \\
c \,&=\, \sigma_1\,u_2, \quad d \,=\, \sigma_2\,u_1
\end{align}
\]
\[
\operatorname{det}\left( \pmb{\operatorname{A}}_E \right) \:=\: a\,d \,-\, b\,c \:=\: \sigma_1\, \sigma_2
\]

σ1 and σ2 are factors which give us the lengths of the principal axes of the ellipse. σ1 and σ2 have positive values. We therefore demand:

\[
\operatorname{det}\left( \pmb{\operatorname{A}}_E \right) \:=\: a\,d \,-\, b\,c \:\ge\: 0
\]

Ok, we have defined an ellipse via a matrix AE, whose coefficients are directly based on geometrical properties. But as said: Often an ellipse is described by a an equation with quadratic terms in x and y coordinates of data points. The quadratic form has its background in algebraic properties of conic sections. As a next step we derive such a quadratic equation and relate the coefficients of the quadratic polynomial with the elements of our matrix AE. The result will in turn define another very useful matrix Aq.

Quadratic forms – Case 1: Centered ellipse, principal axes aligned with ECS-axes

We start with a simple case. We take a so called axis-parallel ellipse which results from a scaling matrix DE, only. I.e., in this case, the rotation matrix is assumed to be just the identity matrix. We can omit it from further calculations:

\[
\pmb{\operatorname{R}}_{\phi} \:=\:
\begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \,=\, \pmb{\operatorname{I}}, \quad u_1 \,=\, 1,\: u_2 \,=\, 0, \: \phi = 0
\]

We need an expression in terms of (xe, ye). To get quadratic terms of vector components it often helps to invoke a scalar product. The scalar product of a vector with itself gives us the squared norm or length of a vector. In our case the norms of the inversely re-scaled vectors obviously have to fulfill:

\[
\left[\, \pmb{\operatorname{D}}_E^{-1} \circ \begin{pmatrix} x_e \\ y_e \end{pmatrix} \, \right]^T \,\bullet \, \left[\, \pmb{\operatorname{D}}_E^{-1} \circ \begin{pmatrix} x_e \\ y_e \end{pmatrix} \,\right] \:=\: 1
\]

(The bullet represents the scalar product of the vectors.) This directly results in:

\[
{1 \over \sigma_1^2} x_e^2 \, + \, {1 \over \sigma_2^2} y_e^2 \:=\: 1
\]

We eliminate the denominator to get a convenient quadratic form:

\[
\lambda_2\,x_e^2 \,+\, \lambda_1\, y_e^2 \:=\: \lambda_1 \lambda_2 \quad \left( = \, \operatorname{det}\left(\pmb{\operatorname{D}}_E\right) \right) \phantom{\huge{(}}
\]

If we were given the quadratic form more generally by coefficients α, β and γ

\[
\alpha \,x_e^2 \,+\, \beta\, x_e y_e \,+\, \gamma\, y_e^2 \:=\: \delta
\]

we could directly relate these coefficients with the geometrical properties of our ellipse:

Axis-parallel ellipse:

\[ \begin{align}
\alpha \,&=\, c^2 \,=\, \sigma_2^2 \,=\, \lambda_2 \\
\gamma \,&=\, a^2 \,=\, \sigma_1^2 \,=\, \lambda_1 \\
\beta \,&=\, b \,=\, c \,=\, 0 \\
\delta \,&=\, a^2\, d^2 \,=\, \sigma_1^2 \, \sigma_2^2 \,=\, \lambda_1 \lambda_2 \\
\phi &= 0
\end{align}
\]

I.e., we can directly derive σ1, σ2 and φ from the coefficients of the quadratic form. But an axis-parallel ellipse is a very simple ellipse. Things get more difficult for a rotated ellipse.

Quadratic forms – Case 2: General centered and rotated ellipse

We perform the same trick to get a quadratic polynomial with the vectors ve of a rotated ellipse:

\[
\left[ \,\pmb{\operatorname{A}}_E^{-1} \circ \begin{pmatrix} x_e \\ y_e \end{pmatrix} \, \right]^T \, \bullet \,
\left[ \, \pmb{\operatorname{A}}_E^{-1} \circ \begin{pmatrix} x_e \\ y_e \end{pmatrix} \,\right] \:=\: 1
\]

I skip the lengthy, but simple algebraic calculation. We get (with our matrix elements a, b, c, d):

\[
\left( c^2 \,+\, d^2 \right)\,x_e^2 \,\, – \,\, 2\left( a\,c\, +\, b\,d \right)\,x_e y_e \,\, + \,\, \left(a^2 \,+\, b^2\right)\,y_e^2
\:=\: \sigma_1^2 \, \sigma_2^2
\]

The rotation has obviously lead to mixing of components in the polynomial. The coefficient for xeye is > 0 for the non-trivial case.

Quadratic form: A matrix equation to define an ellipse

We rewrite our equation again with general coefficients α, β and γ

\[
\alpha\,x_e^2 \, + \, \beta \, x_e y_e \, + \, \gamma \, y_e^2 \:=\: \delta
\]

These are coefficients which may come from some theory or from averages of numerical data. The quadratic polynomial can in turn be reformulated as a matrix operation with a symmetric matrix Aq:

\[
\pmb{v}_e^T \circ \pmb{\operatorname{A}}_q \circ \pmb{v}_e^T \:=\: \delta
\]

with

\[ \pmb{\operatorname{A}}_q \:=\:
\begin{pmatrix} \alpha & \beta / 2 \\ \beta / 2 & \gamma \end{pmatrix}
\]
\[ \pmb{\operatorname{A}}_q \:=\:
\begin{pmatrix} c^2 \,+\, d^2 & a\,c \, +\, b\,d \\ a\,c \, +\, b\,d & a^2 \,+\, b^2 \end{pmatrix}
\]
\[ \begin{align}
\alpha \:&=\: c^2 \,+\, d^2 \:=\: \sigma_1^2 u_2^2 \, + \, \sigma_2^2 u_1^2 \\
\gamma \:&=\: a^2 \,+\, b^2 \:=\: \sigma_1^2 u_1^2 \, + \, \sigma_2^2 u_2^2 \\
\beta \:&=\: – 2\left(a c \,+\, b d \right) \:=\: -2 \left( \sigma_1^2 u_1 u_2 \, + \, \sigma_2^2 u_1 u_2 \right) \\
\delta \:&=\: \left(ad \,-\, bc\right)^2 \:=\: \sigma_1^2 \, \sigma_2^2
\end{align}
\]

Note also:

\[ \begin{align}
\alpha \,+\, \gamma \:&=\: a^2 \,+\, b^2 \,+\, c^2 \,+\, d^2 \:=\: \sigma_1^2 \, + \, \sigma_2^2 \\
\alpha \,-\, \gamma \:&=\: a^2 \,+\, b^2 \,-\, c^2 \,-\, d^2 \:=\: \sigma_1^2 \, \operatorname{cos}\left(2\phi\right) \, + \, \sigma_2^2 \, \operatorname{cos}\left(2\phi\right) \end{align}
\]

These terms are intimately related to the geometrical data; expect them to play a major role in further considerations.

With the help of the coefficients of AE we can also show that det(Aq) > 0:

\[
\operatorname{det}\left( \pmb{\operatorname{A}}_q \right) \:=\: \left(\alpha\gamma \,-\, {1\over 4}\beta^2\right) \:=\: \left(bc \,-\, ad\right)^2 \, \ge \, 0
\]

Thus Aq is an invertible matrix (as was to be expected).

Above we have got α, β, γ, δ as some relatively simple functions of a, b, c, d. The inversion is not so trivial and we do not even try it here. Instead we focus on how we can express σ1, σ2 and φ as functions of either (a, b, c, d) or (α, β, γ, δ).

How to derive σ1, σ2 and φ from the coefficients of AE or Aq in the general case?

Let us assume we have (numerical) data for the coefficients of the quadratic form. Then we may want to calculate values for the length of the principal axes and the rotation angle φ of the corresponding ellipse. There are two ways to derive respective formulas:

  • Approach 1: Use trigonometric relations to directly solve the equation system.
  • Approach 2: Use an eigenvector decomposition of Aq.

Both ways are fun!

Direct derivation of σ1, σ2 and φ by using trigonometric relations

We start with the hard tour, namely by solving equations for λ1, λ2 and φ directly. This requires some knowledge in trigonometry. So far, we know the following:

\[ \begin{align}
\gamma \:&=\: a^2 \,+\, b^2 \:=\: \lambda_1 \operatorname{cos}^2\phi \, + \, \lambda_2 \operatorname{sin}^2\phi \\
\alpha \:&=\: c^2 \,+\, d^2 \:=\: \lambda_2 \operatorname{sin}^2\phi \, + \, \lambda_2 \operatorname{cos}^2\phi \\
\beta \:&=\: – 2\left(a c \,+\, b d \right) \:=\: -2 \left( \lambda_1 \,-\, \lambda_2 \right) \operatorname{cos}\phi \, \operatorname{sin}\phi \\
\delta \:&=\: \lambda_1 \, \lambda_2
\end{align}
\]

Trigonometric relations which we can use are:

\[ \begin{align}
\operatorname{sin}(2 \phi) \:&=\: 2 \,\operatorname{cos}(\phi)\, \operatorname{sin}(\phi) \\
\operatorname{cos}(2 \phi) \:&=\: 2 \,\operatorname{cos}^2(\phi)\, -\, 1 \\
\:&=\: 1 \,-\, 2\,\operatorname{sin}^2(\phi) \\
\:&=\: \operatorname{cos}^2(\phi)\, -\, \operatorname{sin}^2(\phi)
\end{align}
\]

Without loosing generality we assume

\[
\lambda_1 \:\ge \lambda_2
\]

In the end results would only differ by a rotation of π/2, if we had chosen otherwise. This leads to

\[ \begin{align}
2 \gamma \:&=\: 2\left( a^2 \,+\, b^2 \right) \:=\: \lambda_1 \left( 1 \,+\, \operatorname{cos}(2\phi) \right) \, + \,
\lambda_2 \left( 1 \,-\, \operatorname{cos}(2\phi) \right) \\
2 \alpha \:&=\: 2 \left( c^2 \,+\, d^2 \right) \:=\: \lambda_1 \left( 1 \,-\, \operatorname{cos}(2\phi) \right) \, + \,
\lambda_2 \left( 1 \,+\, \operatorname{cos}(2\phi) \right) \\
– \beta \:&=\: 2 \left(a c \,+\, b d \right) \:=\: \left( \lambda_1 \,-\, \lambda_2 \right) \operatorname{sin}(2\phi)
\end{align}
\]

We rearrange terms and get:

\[ \begin{align}
\left( \lambda_1 \,-\, \lambda_2 \right) \operatorname{cos}(2\phi) \,+\, \lambda_1 \,+\, \lambda_2 \:&=\: 2 \left( a^2 \,+\, b^2 \right) \\
\left( \lambda_1 \,-\, \lambda_2 \right) \operatorname{cos}(2\phi) \,-\, \lambda_1 \,-\, \lambda_2 \:&=\: – 2 \left( c^2 \,+\, d^2 \right) \\
\left( \lambda_1 \,-\, \lambda_2 \right) \operatorname{sin}(2\phi) \:&=\: 2 \left( a\,c \,+\, b\,d \right)
\end{align}
\]

Let us define some further variables before we add and subtract the first two of the above equations:

\[ \begin{align}
r \:&=\: {1 \over 2} \left( a^2 \,+\, b^2 \,+\, c^2 \,+\, d^2 \right) \:=\: {1 \over 2} \left( \gamma \,+\, \alpha \right) \\
s_1 \:&=\: {1 \over 2} \left( a^2 \,+\, b^2 \,-\, c^2 \,-\, d^2 \right) \:=\: {1 \over 2} \left( \gamma \,-\, \alpha \right) \\
s_2 \:&=\: \left( a\,c \,+\, b\,d \right) \:=\: – {1 \over 2} \, \beta \\
\pmb{s} \:&=\: \begin{pmatrix} s_1 \\ s_2 \end{pmatrix} \:=\: {1 \over 2} \begin{pmatrix} \gamma \,-\, \alpha \\ – \beta \end{pmatrix} \phantom{\Huge{(}} \\
s \:&=\: \sqrt{ s_1^2 \,+\, s_2^2 } \:=\: {1 \over 2} \left[ \, \beta^2 \,+\, \left( \gamma \,-\, \alpha \right)^2 \,\right]^{1/2} \phantom{\huge{(}}
\end{align}
\]

Then adding two of the equations with the sin2φ and cos2φ above and using the third one results in:

\[
{1 \over 2} \left( \lambda_1 \,-\, \lambda_2 \right) \begin{pmatrix} \operatorname{cos}(2\phi) \\ \operatorname{sin}(2\phi) \end{pmatrix} \:=\: \begin{pmatrix} s_1 \\ s_2 \end{pmatrix}
\]

Taking the vector norm on both sides (with λ1 ≥ λ2) and adding two of the equations above results in:

\[ \begin{align}
\lambda_1 \,\, – \,\, \lambda_2 \:&=\: 2 s \:\:=\: \left[ \, \beta^2 \,+\, \left( \gamma \,-\, \alpha \right) \,\right]^{1/2} \\
\lambda_1 \, + \, \lambda_2 \:&=\: 2 r \:\:=\: \gamma \,+\, \alpha
\end{align}
\]

This gives us:

\[ \begin{align}
\sigma_1^2 \,=\, \lambda_1 \:&=\: r \,+\, s \\
\sigma_2^2 \,=\, \lambda_2 \:&=\: r \,-\, s
\end{align}
\]

In terms of a, b, c, d:

\[ \begin{align}
\sigma_1^2 \,=\, \lambda_1 \:&=\: {1 \over 2} \left[ \, a^2+b^2+c^2 +d^2 \,+\, \left[ 4 (ac + bd)^2 \, +\, \left( c^2+d^2 -a^2 -b^2\right)^2 \, \right]^{1/2} \right] \\
\sigma_2^2 \,=\, \lambda_2 \:&=\: {1 \over 2} \left[ \, a^2+b^2+c^2 +d^2 \,-\, \left[ 4 (ac + bd)^2 \, +\, \left( c^2+d^2 -a^2 -b^2\right)^2 \, \right]^{1/2} \right]
\end{align}
\]

Who said that life has to be easy? In terms of α, β, γ, δ it looks a bit better:

\[ \begin{align}
\sigma_1^2 \,=\, \lambda_1 \:&=\: {1 \over 2} \left[ \, \left( \gamma \,+\, \alpha \right) \,+\, \left[ \beta^2 \, +\, \left( \gamma \,-\, \alpha \right)^2 \, \right]^{1/2} \right] \\
\sigma_2^2 \,=\, \lambda_2 \:&=\: {1 \over 2} \left[ \, \left( \gamma \,+\, \alpha \right) \,-\, \left[ \beta^2 \, +\, \left( \gamma \,-\, \alpha \right)^2 \, \right]^{1/2} \right]
\end{align}
\]

The reader can convince himself that with the definitions above we do indeed reproduce

\[
\lambda_1 \, \lambda_2 \:=\: r^2 \,-\, s^2 \:=\: \left(a\, d\,+\, b\,c \right)^2 \,=\, \operatorname{det}\left(\pmb{\operatorname{A}}_E\right)^2
\]

Determination of the inclination angle φ

For the determination of the angle φ we use:

\[
\begin{pmatrix} \operatorname{cos}(2\phi) \\ \operatorname{sin}(2\phi) \end{pmatrix} \:=\: {1 \over s} \begin{pmatrix} s_1 \\ s_2 \end{pmatrix}
\]

If we choose

\[
-\pi/2 \,\lt\, \phi \le \pi/2
\]

we get:

\[
\phi \:=\: {1 \over 2} \operatorname{arctan}\left({s_2 \over s_1}\right) \:=\: – {1 \over 2} \operatorname{arctan}\left( { – 2\left(ac \,+\, bd\right) \over (a^2 \,+\, b^2 \,-\, c^2 \,-\, d^2) } \right)
\]
\[
\phi \:=\: – {1 \over 2} \operatorname{arctan}\left( {\beta \over \gamma \,- \alpha } \right)
\]

Or equivalently with respect to α, β, γ:

\[
\operatorname{sin}\left(2\phi\right) \:=\: – \, {\beta \over \left[ \beta^2 \, +\, \left( \gamma \,-\, \alpha \right)^2 \, \right]^{1/2} }
\]
\[
\phi \:=\: {1 \over 2} \operatorname{arcsin}\left({s_2 \over s}\right) \:=\: – {1 \over 2} \operatorname{arcsin}\left( {\beta \over \left[ \beta^2 \, +\, \left( \gamma \,-\, \alpha \right)^2 \, \right]^{1/2} } \right)
\]

Note: All in all there are four different solutions. The reason is that we alternatively could have requested λ2 ≥ λ1 and also chosen the angle π + φ. So, the ambiguity is due to a selection of the considered principal axis and rotational symmetries.

In the special case that we have a circle

\[
\lambda_1 \,=\, \lambda_2 \,=\, r
\]

and then, of course, any angle φ will be allowed.

2nd way to a solution for σ1, σ2 and φ via eigendecomposition

For our second way of deriving formulas for σ1, σ2 and φ we use some linear algebra. This way is interesting for two reasons: It indicates how we can use the Python “linalg”-package together with Numpy to get results numerically. In addition we get familiar with a representation of the ellipse in a properly rotated ECS.

Above we have written down a symmetric matrix Aq describing an operation on the position vectors of points on our rotated ellipse:

\[
\pmb{v}_e^T \circ \pmb{\operatorname{A}}_q \circ \pmb{v}_e^T \:=\: \delta
\]

We know from linear algebra that every symmetric matrix can be decomposed into a product of orthogonal matrices O, OT and a diagonal matrix. This reflects the so called eigendecomposition of a symmetric matrix. It is a unique decomposition in the sense that it has a uniquely defined solution in terms of the coefficients of the following matrices:

\[
\pmb{\operatorname{A}}_q \:=\: \pmb{\operatorname{O}} \circ \pmb{\operatorname{D}}_q \circ \pmb{\operatorname{O}}^T
\]

with

\[
\pmb{\operatorname{D}}_{q} \:=\: \begin{pmatrix} \lambda_{u} & 0 \\ 0 & \lambda_{d} \end{pmatrix}
\]

The coefficients λu and λd are eigenvalues of both Dq and Aq. Reason: Orthogonal matrices do not change eigenvalues of a transformed matrix. So, the diagonal elements of Dq are the eigenvalues of Aq. Linear algebra also tells us that the columns of the matrix O are given by the components of the normalized eigenvectors of Aq.

We can interpret O as a rotation matrix Rψ for some angle ψ:

\[
\pmb{v}_e^T \circ \pmb{\operatorname{A}}_q \circ \pmb{v}_e^T \:=\: \pmb{v}_e^T \circ
\pmb{\operatorname{R}}_{\psi} \circ \pmb{\operatorname{D}}_q \circ \pmb{\operatorname{R}}_{\psi}^T \circ \pmb{v}_e \:=\: \delta
\]

This means

\[
\left[ \pmb{\operatorname{R}}_{-\psi} \circ \pmb{v}_e \right]^T \circ
\pmb{\operatorname{D}}_q \circ \left[ \pmb{\operatorname{R}}_{-\psi} \circ \pmb{v}_e \right] \:=\: \delta \:=\: \sigma_1^2\, \sigma_2^2
\]
\[
\pmb{v}_{-\psi}^T \circ \pmb{\operatorname{D}}_q \circ \pmb{v}_{-\psi} \:=\: \sigma_1^2\, \sigma_2^2
\]

The whole operation tells us a simple truth, which we are already familiar with:

By our construction procedure for a rotated ellipse we know that a rotated ECS exists, in which the ellipse can be described as the result of a scaling operation (along the coordinate axes of the rotated ECS) applied to a unit circle. (This ECS is, of course, rotated by an angle φ against our working ECS in which the ellipse appears rotated.)

Indeed:

\[
\left(\, x_{-\psi}, \,y_{-\psi}\,\right) \circ \begin{pmatrix} \lambda_u & 0 \\ 0 & \lambda_d \end{pmatrix} \circ \begin{pmatrix} x_{-\psi} \\ y_{-\psi} \end{pmatrix} \:=\: \sigma_1^2\, \sigma_2^2
\]
\[
{\lambda_u \over \sigma_1^2\, \sigma_2^2} \, x_{-\psi}^2 \, + \, {\lambda_d \over \sigma_1^2\, \sigma_2^2} \, y_{-\psi}^2 \:=\: 1
\]

We know exactly what angle ψ by which we have to rotate our ECS to get this result: ψ = φ. Therefore:

\[ \begin{align}
x_{-\psi} \:&=\: x_c, \\
y_{-\psi} \:&=\: y_c, \\
\lambda_u \:&=\: \lambda_2 \,=\, \sigma_2^2, \\
\lambda_d \:&=\: \lambda_1 \,=\, \sigma_1^2
\end{align}
\]

This already makes it plausible that the eigenvalues of our symmetric matrix Aq are just λ1 and λ2.

Mathematically, a lengthy calculation will indeed reveal that the eigenvalues of a symmetric matrix Aq with coefficients α, 1/2*β and γ have the following form:

\[
\lambda_{u/d} \:=\: {1 \over 2} \left[\, \left(\alpha \,+\, \gamma \right) \,\pm\, \left[ \beta^2 + \left(\gamma \,-\, \alpha \right)^2 \,\right]^{1/2} \, \right]
\]

This is, of course, exactly what we have found some minutes ago by directly solving the equations with the trigonometric terms.

We will prove the fact that these indeed are valid eigenvalues in a minute. Let us first look at respective eigenvectors ξ1/2. To get them we must solve the equations resulting from

\[
\left( \begin{pmatrix} \alpha & \beta / 2 \\ \beta / 2 & \gamma \end{pmatrix} \,-\, \begin{pmatrix} \lambda_{1/2} & 0 \\ 0 & \lambda_{1/2} \end{pmatrix} \right) \,\circ \, \pmb{\xi_{1/2}} \:=\: \pmb{0},
\]

with

\[
\pmb{\xi_1} \,=\, \begin{pmatrix} \xi_{1,x} \\ \xi_{1,y} \end{pmatrix}, \quad \pmb{\xi_2} \,=\, \begin{pmatrix} \xi_{2,x} \\ \xi_{2,y} \end{pmatrix}
\]

Again a lengthy calculation shows that the following vectors fulfill the conditions (up to a common factor in the components):

\[ \begin{align}
\lambda_1 \: &: \quad \pmb{\xi_1} \:=\: \left(\, {1 \over \beta} \left( (\alpha \,-\, \gamma) \,-\, \left[\, \beta^2 \,+\, \left(\gamma \,-\, \alpha \right)^2\,\right]^{1/2} \right), \: 1 \, \right)^T \\
\lambda_2 \: &: \quad \pmb{\xi_2} \:=\: \left(\, {1 \over \beta} \left( (\alpha \,-\, \gamma) \,+\, \left[\, \beta^2 \,+\, \left(\gamma \,-\, \alpha \right)^2\,\right]^{1/2} \right), \: 1 \, \right)^T
\end{align}
\]

for the eigenvalues

\[ \begin{align}
\lambda_1 \:&=\: {1 \over 2} \left(\, \left(\alpha \,+\, \gamma \right) \,+\, \left[ \beta^2 \,+\, \left(\gamma \,-\, \alpha \right)^2 \,\right]^{1/2} \, \right) \\
\lambda_2 \:&=\: {1 \over 2} \left(\, \left(\alpha \,+\, \gamma \right) \,-\, \left[ \beta^2 \,+\, \left(\gamma \,-\, \alpha \right)^2 \,\right]^{1/2} \, \right) \\
\end{align}
\]

The T at the formulas for the vectors symbolizes a transposition operation.

Note that the vector components given above are not normalized. This is important for performing numerical checks as Numpy and linear algebra programs would typically give you normalized eigenvectors with a length = 1. But you can easily compensate for this by working with

\[ \begin{align}
\lambda_1 \: &: \quad \pmb{\xi_1^n} \:=\: {1 \over \|\pmb{\xi_1}\|}\, \pmb{\xi_1} \\
\lambda_2 \: &: \quad \pmb{\xi_2^n} \:=\: {1 \over \|\pmb{\xi_2}\|}\, \pmb{\xi_2}
\end{align}
\]

Proof for the eigenvalues and eigenvector components

We just prove that the eigenvector conditions are e.g. fulfilled for the components of the second eigenvector ξ2 and λ2 = λd.

\[ \begin{align}
\left(\alpha \,-\, \lambda_2 \right) * \xi_{2,x} \,+\, {1 \over 2} \beta * \xi_{2,y} \,&=\, 0 \\
{1 \over 2} \beta * \xi_{2,x} \,+\, \left( \gamma \,-\, \lambda_2 \right) * \xi_{2,y} \,&=\, 0
\end{align}
\]

(The steps for the first eigenvector are completely analogous).

We start with the condition for the first component

\[ \begin{align}
&\left( \alpha \,-\,
{1\over 2}\left[\,\left(\alpha \, + \, \gamma\right) \,+\, \left[ \beta^2 \,+\, \left( \alpha \,-\, \gamma \right)^2 \right]^{1/2} \right] \right) * \\
& {1 \over \beta}\,
\left[\, \left(\alpha \,-\, \gamma\right) \,+\, \left[ \beta^2 \,+\, \left(\alpha \,-\, \gamma \right)^2 \right]^{1/2} \,\right] \,+\, {\beta \over 2 }
\,=\, 0
\end{align}
\]
\[ \begin{align}
& {1 \over 2 } \left[ \left(\alpha \,-\,\gamma\right) \,+\, \left[ \beta^2 \,+\, \left( \alpha \,-\, \gamma \right)^2 \right]^{1/2} \right] * \\
& {1 \over \beta}\,
\left[\, \left(\alpha \,-\, \gamma\right) \,+\, \left[ \beta^2 \,+\, \left(\alpha \,-\, \gamma \right)^2 \right]^{1/2} \,\right] \,+\, {\beta \over 2 }
\,=\, 0
\end{align}
\]
\[
{1 \over 2 \beta} \left[ (\alpha \,-\,\gamma)^2 \,-\, \beta^2 \,-\, (\alpha \,-\,\gamma)^2 \right] \,+\, {\beta \over 2 } \,=\, 0
\]

The last relation is obviously true. You can perform a similar calculation for the other eigenvector component:

\[ \begin{align}
{1 \over 2} \, \beta & {1 \over \beta}\,
\left[\, \left(\alpha \,-\, \gamma\right) \,+\, \left[ \beta^2 \,+\, \left(\alpha \,-\, \gamma \right)^2 \right]^{1/2} \,\right] \,+\, \\
&
\left( \gamma \, -\,
{1\over 2}\left[\,\left(\alpha \, + \, \gamma\right) \,+\, \left[ \beta^2 \,+\, \left( \alpha \,-\, \gamma \right)^2 \right]^{1/2} \right] \right) * 1 \,=\, 0
\end{align}
\]

Thus:

\[ \begin{align}
&{1 \over 2} \, \left(\alpha \,-\, \gamma\right) \,+\, {1 \over 2} \left[ \beta^2 \,+\, \left(\alpha \,-\, \gamma \right)^2 \right]^{1/2} \,-\, \\
&{1\over 2}\left(\alpha \,-\, \gamma\right) \,-\, {1\over 2}\left[ \beta^2 \,+\, \left( \alpha \,-\, \gamma \right)^2 \right]^{1/2} \,=\, 0
\end{align}
\]

True, again. In a very similar exercise one can show that the scalar product of the eigenvectors is equal to zero:

\[ \begin{align}
& {1 \over \beta}\,
\left[\, \left(\alpha \,-\, \gamma\right) \,+\, \left[ \beta^2 \,+\, \left(\alpha \,-\, \gamma \right)^2 \right]^{1/2} \,\right] \,+\, {\beta \over 2 } \,*\, \\
& {1 \over \beta}\,
\left[\, \left(\alpha \,-\, \gamma\right) \,-\, \left[ \beta^2 \,+\, \left(\alpha \,-\, \gamma \right)^2 \right]^{1/2} \,\right] \,+\, {\beta \over 2 } \,+\, 1\,*\,1 \\
&\, {1 \over \beta^2}\, *\, (-\beta^2) \,+\, 1 \,=\, 0
\end{align}
\]

I.e.:

\[
\pmb{\xi_1} \bullet \pmb{\xi_1} \,=\, \left( \xi_{1,x}, \, \xi_{1,y} \right) \circ \begin{pmatrix} \xi_{2,x} \\ \xi_{2,y} \end{pmatrix} \,= \, 0,
\]

which means that the eigenvectors are perpendicular to each other. Exactly, what we expect for the orientations of the principal axes of an ellipse against each other.

Rotation angle from coefficients of Aq

We still need a formula for the rotation angle(s). From linear algebra results related to an eigendecomposition we know that the orthogonal (rotation) matrices consist of columns of the normalized eigenvectors. With the components given in terms of our un-rotated ECS, in which we basically work. These vectors point along the principal axes of our ellipse. Thus the components of these eigenvectors define our aspired rotation angles of the ellipse’s principal axes against the x-axis of our ECS.

Let us prove this. By assuming

\[ \begin{align}
\operatorname{cos}(\phi_1) \,&=\, \xi_{1,x}^n \\
\operatorname{sin}(\phi_1) \,&=\, \xi_{1,y}^n
\end{align}
\]

and using

\[
\operatorname{cos}(2\phi_1) \,=\, 2\, \operatorname{sin}(\phi_1) \, \operatorname{cos}(\phi_1)
\]

we get:

\[ \begin{align}
\operatorname{sin}(2 \phi_1) \,=\,
2 * { \xi_{1,x} * \xi_{1,y} \over \left[\, \xi_{1,x}^2 \, + \, \xi_{1,y}^2 \,\right]^{1/2} }
\end{align}
\]

and thus

\[ \begin{align}
\operatorname{sin}(2 \phi_1) \,&=\,
2 \, { {1 \over \large{\beta}} \left( (\alpha \,-\, \gamma) \,-\, \left[\, \beta^2 \,+\, \left(\gamma \,-\, \alpha \right)^2\,\right]^{1/2} \right) \,*\, 1
\over
\left[\, \left( {1 \over \large{\beta}} \left( (\alpha \,-\, \gamma) \,-\, \left[\, \beta^2 \,+\, \left(\gamma \,-\, \alpha \right)^2\,\right]^{1/2} \right) \right)^2
\,+\, 1^2 \right]^{1/2} } \\
&=\, 2\, { {1 \over \large{\beta}} \left( t \,-\, z \right)
\over
{1 \over \large{\beta}^2 \phantom{\large{]}} } \left[\, \beta^2 \,+\, \left(\, t \,-\, z \,\right)^2 \right]^2 }
\end{align}
\]

with

\[ \begin{align}
t \,&=\, (\alpha \,-\, \gamma) \\
z \,&=\, \left[\, \beta^2 \,+\, \left(\gamma \,-\, \alpha \right)^2\,\right]^{1/2}
\end{align}
\]

This looks very differently from the simple expression we got above. And a direct approach is cumbersome. The trick is multiply nominator and denominator by a convenience factor

\[
\left( t \,+\, z \right),
\]

and exploit

\[ \begin{align}
\left( t \,-\, z \right) \, \left( t \,+\, z \right) \,&=\, t^2 \,-\, z^2 \\
\left( t \,-\, z \right) \, \left( t \,+\, z \right) \,&=\, – \beta^2
\end{align}
\]

to get

\[ \begin{align}
&2 * \beta { (t\,-\, z) * ( t\,+\, z) \over \left[ \beta^2 \, + \,( t \,-\, z )^2 \right] * (t \,+\, z) } \\
=\, &2 * \beta { – \beta^2 \over \beta^2 (t\,+\,z) \,-\, \beta^2 (t\,-\,z) } \\
=\, & – {\beta \over \left[\, \beta^2 \,+\, (\alpha \,-\, \gamma)^2 \,\right]^{1/2} }
\end{align}
\]

This means that our 2nd solution approach provides the result

\[
\operatorname{sin}(2 \phi_1) \, =\, – \, { \beta \over \left[\, \beta^2 \,+\, (\alpha \,-\, \gamma)^2 \,\right]^{1/2} }\,,
\]

which is of course identical to the result we got with our first solution approach. It is clear that the second axis has an inclination by φ +- π / 2:

\[
\phi_2\, =\, \phi_1 \,\pm\, \pi/2.
\]

In general the angles have a natural ambiguity of π.

Conclusion

In this post I have shown how one can derive essential properties of centered, but rotated ellipses from matrix-based representations. Such calculations become relevant when e.g. experimental or numerical data only deliver the coefficients of a quadratic form for the ellipse.

We have first established the relation of the coefficients of a matrix that defines an ellipse by a combined scaling and rotation operation with the coefficients of a matrix which defines an ellipse as a quadratic form of the components of position vectors. In addition we have shown how the coefficients of both matrices are related to quantities like the lengths of the principal axes of the ellipse and the inclination of these axes against the x-axis of the Euclidean coordinate system in which the ellipse is described via position vectors. So, if one of the matrices is given we can numerically calculate the ellipse’s main properties.

In the next post of this mini-series

Properties of ellipses by matrix coefficients – II – coordinates of points with extremal y-values

we have a look at the x- and y-coordinates of points on an ellipse with extremal y-values. All in terms of the matrix coefficients we are now familiar with.