raw Math
RAW Math Linear Algebra Statistics

How to plot a covariance error ellipse

Robert Eisele

A typical way to visualize two-dimensional Gaussian-distributed data is to plot a confidence ellipse. Assume the data is drawn as \(D\sim\mathcal{N}(\mu, \Sigma)\), and the goal is to plot an ellipse representing a confidence level \(p\) by computing the radii of the ellipse, its center, and its rotation. The following plot shows randomly drawn data together with the ellipses for \(p\in\{0.9, 0.95, 0.99\}\):

Why the Level Sets Are Ellipses at All

The density of a two-dimensional Gaussian is largest at the mean \(\mu\) and decays in every direction from there, so a natural way to draw "the region containing a given fraction of the probability mass" is to trace a curve of constant density. For an uncorrelated Gaussian with standard deviations \(\sigma_x,\sigma_y\), the density factors into independent one-dimensional Gaussians, and setting it equal to a constant leaves exactly

\[\frac{(x-\mu_x)^2}{\sigma_x^2}+\frac{(y-\mu_y)^2}{\sigma_y^2}=\text{const},\]

which, after centering at the mean, is precisely the equation of an axis-aligned ellipse. So the ellipse is not an arbitrary drawing convention: it falls directly out of the shape of the Gaussian density itself, and the only remaining question is which constant on the right-hand side captures a given confidence level \(p\).

Derivation

Take the uncorrelated case first, with the data already centered at the origin. Then \(x/\sigma_x\) and \(y/\sigma_y\) are two independent standard normal variables, and the left-hand side above, divided by the constant, becomes a sum of two independent squared standard normals — by definition a \(\chi^2\) random variable with two degrees of freedom (one for each dimension). Calling this quantity \(s\), the constant-density curve for confidence level \(p\) is exactly the one enclosing a \(\chi^2_2\)-probability of \(p\):

\[\left(\frac{x}{\sigma_x}\right)^2 + \left(\frac{y}{\sigma_y}\right)^2 = s,\qquad P(\chi^2_2\le s) = p.\]

This value \(s\) is known as the Mahalanobis radius (squared) of the ellipse. In general, finding \(s\) requires the inverse CDF of a \(\chi^2\) distribution, available as a table in any statistics textbook or as the Matlab function s=chi2inv(p, k) for \(k\) degrees of freedom. Two dimensions are a fortunate special case, though: a \(\chi^2\) distribution with exactly two degrees of freedom is the same as an exponential distribution with rate \(\frac12\), whose CDF has the closed form \(P(\chi^2_2\le s)=1-e^{-s/2}\). Setting this equal to \(p\) and solving for \(s\) gives a formula with no table or numerical inversion needed at all:

\[s = -2\log(1-p)\]

For \(p=0.9\) this gives \(s=4.6057\), for \(p=0.95\) it gives \(s=5.99146\), and for \(p=0.99\) it gives \(s=9.21034\) — matching a \(\chi^2\) table exactly, but derived in closed form instead of looked up. The ellipse can then be drawn with radii \(\sigma_x\sqrt{s}\) and \(\sigma_y\sqrt{s}\).

Generalization for a Given Covariance

In the general case, the covariances \(\sigma_{xy}\) and \(\sigma_{yx}\) are not zero, so the natural coordinate system of the ellipse is not axis-aligned, and the density's exponent is the quadratic form \((\mathbf{x}-\mu)^T\Sigma^{-1}(\mathbf{x}-\mu)\) instead of the simple sum of squares above. The trick is to reduce this back to the already-solved uncorrelated case by a change of basis: since the covariance matrix \(\Sigma=\left(\begin{array}{cc}\sigma_x^2&\sigma_{xy}\\\sigma_{yx}&\sigma_y^2\end{array}\right)\) is symmetric, it has an orthonormal eigenbasis \(\mathbf{q}_1,\mathbf{q}_2\) with eigenvalues \(\lambda_1,\lambda_2\), and rotating the data into this basis via \(\mathbf{y}=\mathbf{Q}^T(\mathbf{x}-\mu)\) turns the covariance matrix diagonal, \(\mathrm{Cov}(\mathbf{y})=\mathrm{diag}(\lambda_1,\lambda_2)\) — exactly the uncorrelated case already solved, just expressed along the rotated axes \(\mathbf{q}_1,\mathbf{q}_2\) instead of the original x- and y-axes. By definition, a (non-degenerate) covariance matrix is symmetric positive definite, so both \(\lambda_1,\lambda_2>0\) act exactly like squared standard deviations along these new axes, giving radii \(\sqrt{\lambda_1 s}\) and \(\sqrt{\lambda_2 s}\) oriented along \(\mathbf{q}_1\) and \(\mathbf{q}_2\) rather than along x and y. Since \(\Sigma\) is symmetric positive definite, its eigendecomposition already coincides with its singular value decomposition, so \(\mathbf{q}_1,\mathbf{q}_2\) are equally the left and right singular vectors of \(\Sigma\), and \(\lambda_1,\lambda_2\) its singular values.

The actual radii of the ellipse are therefore \(\sqrt{\lambda_1}\) and \(\sqrt{\lambda_2}\), the two eigenvalues \(\lambda_1\) and \(\lambda_2\) of the scaled covariance matrix \(s\cdot\Sigma\). Since \(\text{eigenvalue}(s\cdot\Sigma)=s\cdot\text{eigenvalue}(\Sigma)\) (scaling a matrix scales its eigenvalues by the same factor, directly from the eigenvalue equation), the eigenvalues of \(\Sigma\) can be scaled by \(s\) after the fact rather than recomputing them from scratch for every confidence level.

Matlab Implementation

With all this in hand, the procedure can be implemented in Matlab directly:

function plotErrorEllipse(mu, Sigma, p)

    if nargin == 2
      p = 0.95; 
    end

    s = -2 * log(1 - p);

    [V, D] = eig(Sigma * s);

    t = linspace(0, 2 * pi);
    a = (V * sqrt(D)) * [cos(t(:))'; sin(t(:))'];

    plot(a(1, :) + mu(1), a(2, :) + mu(2));

JavaScript Implementation

In JavaScript there is no built-in linear algebra support, so the eigenvalues must be computed explicitly. Rather than drawing the ellipse with many short line segments, it is more efficient to use the canvas ellipse function or the SVG ellipse tag, for which the parameters need to be computed explicitly.

For a covariance matrix of the form \(\Sigma=\left(\begin{array}{cc}a&b\\b&d\end{array}\right)\), the two eigenvalues follow from the closed-form \(2\times 2\) eigenvalue formula:

\[\lambda_{1,2} = \frac{a+d}{2} \pm \frac{1}{2}\sqrt{(a-d)^2+4b^2},\qquad \lambda_1\ge\lambda_2.\]

Computing an eigenvector and normalizing it, as the general derivation above does, is more work than a \(2\times 2\) problem actually needs. Instead of finding an eigenvector first and reading its angle off afterward, the angle \(\theta\) that diagonalizes the quadratic form can be found directly: substitute the rotated coordinates \(x=x'\cos\theta-y'\sin\theta,\ y=x'\sin\theta+y'\cos\theta\) into \(ax^2+2bxy+dy^2\) and multiply everything out. Every one of the three terms turns into a combination of \(x'^2\), \(y'^2\) and \(x'y'\); collecting only the coefficient in front of the surviving \(x'y'\) term (the one that has to vanish for the rotated form to be a clean \(\lambda_1x'^2+\lambda_2y'^2\)) gives

\[\underbrace{-2a\sin\theta\cos\theta}_{\text{from }ax^2} + \underbrace{2b(\cos^2\theta-\sin^2\theta)}_{\text{from }2bxy} + \underbrace{2d\sin\theta\cos\theta}_{\text{from }dy^2}.\]

With the double-angle identities \(2\sin\theta\cos\theta=\sin(2\theta)\) and \(\cos^2\theta-\sin^2\theta=\cos(2\theta)\), this collapses to \((d-a)\sin(2\theta)+2b\cos(2\theta)\) — the same rotate-away-the-cross-term move used to bring any conic section into standard position. Setting it to zero is exactly the condition for the cross term to vanish:

\[\tan(2\theta) = \frac{2b}{a-d} \quad\Longrightarrow\quad \theta=\frac12\operatorname{atan2}(2b,\,a-d).\]

Using \(\operatorname{atan2}\) instead of a plain \(\arctan\) avoids a division by zero when \(a=d\) (a circular error region, where every rotation is equally valid) and picks out one specific rotation among the two that are \(90^\circ\) apart. Which of the two gets picked doesn't actually matter: swapping \(\lambda_1\) and \(\lambda_2\) while adding \(90^\circ\) to \(\theta\) describes the exact same ellipse, so there is no need for a case distinction between \(\sigma_x>\sigma_y\) and \(\sigma_x<\sigma_y\) as long as the two radii stay paired with the \(\lambda\) they came from. All that remains is computing the radii of the ellipse.

function plotErrorEllipse(ctx, mu, Sigma, p) {

  p = p || 0.95;

  var s = -2 * Math.log(1 - p);

  var a = Sigma[0][0];
  var b = Sigma[0][1];
  var d = Sigma[1][1];

  var diff = a - d;
  var tmp = Math.sqrt(diff * diff + 4 * b * b);

  var lambda1 = (a + d + tmp) / 2;
  var lambda2 = (a + d - tmp) / 2;
  var theta = 0.5 * Math.atan2(2 * b, diff);

  ctx.ellipse(
          mu[0], mu[1],
          Math.sqrt(s * lambda1),
          Math.sqrt(s * lambda2),
          theta,
          0, Math.PI * 2,
          false);
}

One practical wrinkle when handing \(\theta\) straight to ctx.ellipse: it is only the correct on-screen angle as long as one data unit spans the same number of pixels along both axes. If the plot stretches x and y by different amounts, a direction that sits at \(\theta\) in data space no longer sits at \(\theta\) on screen, because the anisotropic scaling is applied after the rotation is computed. Writing \(k\) for the aspect ratio (pixels per y-unit divided by pixels per x-unit), the direction \((\cos\theta,\sin\theta)\) is stretched to \((\cos\theta, k\sin\theta)\) on screen, so the angle that actually needs to be passed to the drawing function is \(\theta_{\text{screen}}=\operatorname{atan2}(k\sin\theta,\cos\theta)\) rather than \(\theta\) itself.

Instead of eigendecomposing the covariance matrix directly, the same axes and radii can also be read off the singular value decomposition of the mean-centered data matrix, which is the numerically preferred route for principal component analysis on larger datasets.

Related Stuff

raw / Linear Algebra
02 Aug
·
Robert Eisele

Introduction to Singular Value Decomposition

This derives the singular value decomposition (SVD) of a real matrix of any shape directly from the eigendecomposition of the symmetric positive semidefinite matrix A^T A, motivated visually by how any matrix maps the unit circle to an ellipse and by the general rotate-stretch-rotate geometric picture this gives for every matrix. It works through a full numerical 2x2 example, shows how the right and left singular vectors are eigenvectors of A^T A and A A^T and how this completes the four fundamental subspaces picture, illustrated with a diagram of the Strang diagram itself, distinguishes the full and compact SVD, answers precisely when the eigendecomposition of a symmetric matrix coincides with its SVD (positive semidefiniteness) and how singular values differ from eigenvalues in general (including for defective and rotation matrices), derives the uniqueness of the SVD up to sign flips and rotations within repeated singular values with a worked cautionary example of what goes wrong when both sides are eigendecomposed independently, explains the numerically robust way to derive one side from the other, proves the Eckart-Young theorem for optimal low-rank approximation, and lists its properties (including the determinant as a product of singular values) and applications in least squares (including the explicit pseudo-inverse sum formula, exact-solvability criterion, and affine solution space), homogeneous least squares, data compression, and principal component analysis.
raw / Linear Algebra
08 May
·
Robert Eisele

Introduction to Eigenvalues and Eigenvectors

This introduces eigenvalues and eigenvectors as invariant directions of linear transformations. It derives the characteristic polynomial (including the alternate det(lambda I - A) sign convention), works through a full numerical example of a symmetric matrix with a repeated eigenvalue to show the characteristic polynomial, row reduction, and eigenspace computation concretely, explains eigenspaces and multiplicity including a defective, non-diagonalisable counter-example and why real symmetric matrices can never be defective, proves the trace/determinant formulas and the linear independence of eigenvectors for distinct eigenvalues, states Cauchy's interlacing theorem for the eigenvalues of a principal submatrix and the Perron-Frobenius theorem for nonnegative matrices, studies orthogonal matrices (proving every eigenvalue has absolute value 1), rotation matrices with complex eigenvalues and the degrees-of-freedom argument behind them, reflection matrices with real +-1 eigenvalues and a diagram of the mirror line and an eigenvector's invariant direction, and three-dimensional rotations via Euler's rotation theorem and the trace-angle relation, presents the eigenvalue decomposition, extends it to arbitrary real matrix powers and the symmetric square root, derives Gram factors, whitening transforms, and their use in generating correlated random vectors, derives the Rayleigh quotient and its variational characterization of the extreme eigenvalues, shows how repeated application of a matrix is governed by its dominant eigenvalue via power iteration, illustrated with an interactive diagram of successive normalized directions spiraling onto the dominant eigenvector and with the Fibonacci sequence and Binet's formula, proves the Cayley-Hamilton theorem (with a common false-proof pitfall explained) and uses it to derive the determinant of a matrix directly from traces of its powers via Newton's identities, and lists typical applications in linear algebra and engineering, including spectral graph theory and PageRank.
raw / Linear Algebra
02 Aug
·
Robert Eisele

Introduction to Principal Component Analysis

This derives principal component analysis (PCA) from first principles, starting with a visual scatter-plot illustration of why correlated data has a natural long axis, then deriving the first principal component as the direction maximizing projected variance via both the Rayleigh quotient and a Lagrange-multiplier argument, showing both lead to the covariance matrix's dominant eigenvector. It extends this to the full eigendecomposition for all remaining orthogonal components, proves the equivalent reconstruction-error-minimization view via a Pythagorean argument and shows it yields the identical optimal subspace, explains precisely what 'percentage of variance explained' means as a share of the trace of the covariance matrix (with a worked five-component example and an explained-variance bar chart), connects PCA to the singular value decomposition and the Eckart-Young theorem for optimal low-rank reconstruction, contrasts the first principal component with an ordinary least-squares regression line via their different (perpendicular vs. vertical) error criteria, illustrated with a diagram overlaying both fitted lines and both residual directions on the same data, discusses assumptions and limitations (linearity, sensitivity to scale and outliers, decorrelation versus true independence) with kernel PCA and independent component analysis mentioned as extensions, and closes with applications in dimensionality reduction and visualization (illustrated with a schematic MNIST digit-clustering diagram), noise reduction and compression (including the older Hotelling and Karhunen-Loeve transform names for the same construction), and Eigenfaces (Sirovich and Kirby 1987; Turk and Pentland 1991; Lanitis, Taylor and Cootes 1995 for shape eigenfaces) with nearest-neighbor recognition in score space. It also notes the sign- and degeneracy-related non-uniqueness of the principal components, mirroring the same ambiguity already discussed for the singular value decomposition, defines loadings alongside scores, derives the dual-PCA Gram-matrix trick for computing components when there are far fewer observations than variables (cross-linked to a dedicated chapter on EigenGame, a game-theoretic reformulation for computing components at even larger scale), derives whitening (sphering) as a rescaling of the scores to unit variance, and cautions that variance-based component selection is only a proxy for downstream task performance.
All Math