Generate and Visualize Ellipsoids (Ellipses) Nov 17, 2024 visited：

Generate high-dimensional ellipsoids with bounded condition number, and visualize 2D ellipses using Matplotlib.

Generating an Ellipsoid

A ellipsoid in -dimensional Euclidean space can be characterized as:

where is a symmetric positive-definite matrix, and is the center of the ellipsoid. The vector is easy to construct, so the main problem becomes: how do we generate a symmetric p.d. matrix?

Let be a random matrix, then it is almost surely full-rank. Consequently, for any nonzero , we have , and . Therefore, is a p.d. matrix. To generate a p.d. matrix, we can generate the random matrix and compute . The code (using Eigen) is as follows:

typedef Eigen::MatrixXd Matrix;

std::random_device rd;
std::mt19937 engine(rd());

// Generate a real symmetric p.d. matrix
Matrix GeneratePDMatrix(int d) {
	std::uniform_real_distribution<double> distribution(-1.0, 1.0);
	Matrix A = Matrix::NullaryExpr(d, d, [&](){return distribution(engine);});
	return A.transpose() * A;
}

However, in practice this method encountered some numerical problems in high-dimensions (when ). My suspicion is that we lost control on the eigenvalues of . They can be either very small or very large, and the condition number is unbounded. Below is a figure shows the distribution of the eigenvalues for 100 p.d. matrices generated using the approach.

Another approach is to generate an orthonormal matrix , and a positive diagonal matrix . Then the matrix is a symmetric p.d. matrix. The rows of are the eigenvectors, and the entries of are the eigenvalues. In such a way we can control the range of the eigenvalues of .

Let be a randomly generated full-rank matrix. We can make it orthonormal by the Gram-Schmidt process. Let be the projection from to , namely . Let





Then are orthogonal, and the matrix is orthonormal. Now we can generate the eigenvalues as we want, and let . Then is a real symmetric p.d. matrix with “good” eigenvalues (within our desired range).

The time complexity of Gram-Schmidt is , and the matrix multiplication is also . The total process takes time. Code implements this method (using Eigen) is as follows:

typedef Eigen::MatrixXd Matrix;

std::random_device rd;
std::mt19937 engine(rd());

// Project u to v
Matrix Proj(Matrix u, Matrix v) {
	auto a = u.transpose() * v;
	auto b = v.transpose() * v;
	double length = a(0, 0) / b(0, 0);
	return v * length;
}

// The Gram-Schmidt process
Matrix GramSchmidt(const Matrix& A) {
	int d = A.rows();
	Matrix a, P(d, d);
	for (int i = 0; i < d; ++i) {
		a = A(Eigen::all, i);
		for (int j = 0; j < i; ++j) {
			a -= Proj(A(Eigen::all, i), P(Eigen::all, j));
		}
		P(Eigen::all, i) = a;
	}
	for (int i = 0; i < d; ++i) {
		P(Eigen::all, i).normalize();
	}
	return P;
}

// Generate a real symmetric p.d. matrix with eigenvalues in [l, r]
Matrix GeneratePDMatrixWithinRange(int d, double l, double r) {
	std::uniform_real_distribution<double> entry(-1.0, 1.0), eigenval(l, r);
	Matrix A = Matrix::NullaryExpr(d, d, [&]() {return entry(engine);});
	Matrix P = GramSchmidt(A);
	Matrix D(d, d);
	D.setZero();
	for (int i = 0; i < d; ++i) {
		D(i, i) = eigenval(engine);
	}
	return P * D * P.transpose();
}

The figure below shows the distribution of the eigenvalues for the matrices generated using the above method with parameter .

Visualization of Ellipses

Given a 2D ellipse characterized by

where and , how can we visualize it using Matplotlib? Though is an Ellipse patch in the library, it specifies the ellipse by its width, height, and rotation angle. How do we visualize an ellipse using the matrix form?

Let be the square root of . Recall that

The ellipse is equivalent to

Let , then we have

Here is the unit disk. Rearranging the terms gives

In other words, let , the ellipse can be constructed by transforming the unit disk to .

(Alternatively, one can also get and . Both and produces the same ellipse, but the transformation of the unit disk is different.)

To visualize the ellipse, we can sample points on the unit circle, and map them to the contour of the ellipse via .

Code implements this method is as follows:

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import sqrtm

def draw_ellipse(ax, Q, c):
    A = np.linalg.inv(sqrtm(Q))
    ## Alternatively:
    # eigenvals, eigenvecs = np.linalg.eig(Q)
    # A = (1 / np.sqrt(eigenvals)) * eigenvecs
    theta = np.linspace(0, 2 * np.pi, 100)
    xs = (A @ [np.sin(theta), np.cos(theta)]).transpose() + c
    ax.plot(xs[:,0], xs[:,1])

A sample result is given below, where the green and red arrows illustrates the directions of the eigenvectors.

Q = np.array(
    [[1.2360807e+00, -3.1592775e-01], 
    [-3.1592775e-01, 1.3596055e+00]])
c = np.array([2, 1])

fig, ax = plt.subplots()
ax.set_aspect(1)

draw_ellipse(ax, Q, c)

The different between the transformations and are visualized as below, where the green and red arrows illustrate the directions of the eigenvectors.