3D Gaussian Splatting & Covariance — Visual Experiments¶

Goal: Build intuition for two tightly linked notes: why 3DGS represents a scene with explicit Gaussian primitives, and why the covariance matrix is the geometric heart of each primitive.

What we'll explore¶

#	Experiment	Question answered
1	A Gaussian is a soft ellipsoid	What do mean, anisotropy, and orientation actually change in 3D space?
2	Why direct covariance is fragile	Why doesn't 3DGS optimize the six symmetric entries of $\Sigma$ directly?
3	The $R S S^T R^T$ decomposition	What changes when we alter scale only, versus rotation only?
4	Projection to screen space	How does a 3D covariance become a 2D ellipse through $\Sigma' = J W \Sigma W^T J^T$?
5	Alpha compositing and failure modes	Why do covariance choice and depth ordering determine blur, streaks, and occlusion artifacts?

---

Linked theory: research/06-3d-gaussian-splatting.md · research/07-covariance-matrix-3dgs.md · Run: uv run jupyter lab from latent-anything-theory/.

In [1]:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401
from scipy.spatial.transform import Rotation as Rotation3D
import warnings
warnings.filterwarnings('ignore')

np.random.seed(42)
plt.rcParams.update({
    'figure.dpi': 110,
    'font.size': 11,
    'axes.titlesize': 12,
    'axes.labelsize': 11,
    'figure.facecolor': 'white',
    'axes.facecolor': '#f8f9fa',
    'axes.grid': True,
    'grid.alpha': 0.35,
})
print('Setup complete.')

import numpy as np import matplotlib.pyplot as plt from matplotlib.patches import Ellipse from mpl_toolkits.mplot3d import Axes3D # noqa: F401 from scipy.spatial.transform import Rotation as Rotation3D import warnings warnings.filterwarnings('ignore') np.random.seed(42) plt.rcParams.update({ 'figure.dpi': 110, 'font.size': 11, 'axes.titlesize': 12, 'axes.labelsize': 11, 'figure.facecolor': 'white', 'axes.facecolor': '#f8f9fa', 'axes.grid': True, 'grid.alpha': 0.35, }) print('Setup complete.')

Setup complete.

---

A small geometry toolkit¶

The helpers below do four things:

1. Build a covariance from scale + rotation.
2. Convert that covariance into a visible 3D ellipsoid.
3. Project a 3D covariance to screen space with a pinhole-camera Jacobian.
4. Render simple 2D Gaussian splats with front-to-back alpha compositing.

They are intentionally lightweight: the point is to make the math transparent, not to reproduce the full CUDA renderer used by 3DGS.

In [2]:

def rotation_matrix_3d(euler_deg):
    return Rotation3D.from_euler('xyz', euler_deg, degrees=True).as_matrix()


def covariance_from_scale_rotation(scales, euler_deg=(0.0, 0.0, 0.0)):
    scales = np.asarray(scales, dtype=float)
    S = np.diag(scales)
    R = rotation_matrix_3d(euler_deg)
    return R @ S @ S.T @ R.T


def unit_sphere_grid(nu=72, nv=36):
    u = np.linspace(0, 2 * np.pi, nu)
    v = np.linspace(0, np.pi, nv)
    x = np.outer(np.cos(u), np.sin(v))
    y = np.outer(np.sin(u), np.sin(v))
    z = np.outer(np.ones_like(u), np.cos(v))
    return x, y, z


def ellipsoid_surface(cov, mean=(0.0, 0.0, 0.0), radius_level=2.0):
    vals, vecs = np.linalg.eigh(cov)
    vals = np.maximum(vals, 1e-9)
    radii = radius_level * np.sqrt(vals)
    xs, ys, zs = unit_sphere_grid()
    sphere = np.stack([xs, ys, zs], axis=0).reshape(3, -1)
    ellipsoid = (vecs @ np.diag(radii) @ sphere).reshape(3, *xs.shape)
    ellipsoid[0] += mean[0]
    ellipsoid[1] += mean[1]
    ellipsoid[2] += mean[2]
    return ellipsoid


def principal_axes(cov, mean=(0.0, 0.0, 0.0), radius_level=2.0):
    vals, vecs = np.linalg.eigh(cov)
    order = np.argsort(vals)[::-1]
    vals = vals[order]
    vecs = vecs[:, order]
    mean = np.asarray(mean, dtype=float)
    segments = []
    for idx in range(3):
        direction = vecs[:, idx] * radius_level * np.sqrt(max(vals[idx], 1e-9))
        segments.append((mean - direction, mean + direction))
    return segments, vals, vecs


def plot_ellipsoid(ax, cov, mean=(0.0, 0.0, 0.0), title='', color='#4C78A8', alpha=0.28):
    X, Y, Z = ellipsoid_surface(cov, mean=mean)
    ax.plot_surface(X, Y, Z, color=color, alpha=alpha, linewidth=0, shade=True)
    segments, vals, _ = principal_axes(cov, mean=mean)
    axis_colors = ['#E45756', '#72B7B2', '#F58518']
    for (start, end), c in zip(segments, axis_colors):
        ax.plot([start[0], end[0]], [start[1], end[1]], [start[2], end[2]], color=c, lw=2.2)
    ax.scatter([mean[0]], [mean[1]], [mean[2]], color='black', s=28)
    ax.set_title(title)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_box_aspect([1, 1, 1])
    ax.view_init(elev=24, azim=36)
    lim = 3.3
    ax.set_xlim(-lim + mean[0], lim + mean[0])
    ax.set_ylim(-lim + mean[1], lim + mean[1])
    ax.set_zlim(-lim + mean[2], lim + mean[2])
    return vals


def quadratic_form_on_plane(cov, grid_lim=3.0, n=220):
    xs = np.linspace(-grid_lim, grid_lim, n)
    ys = np.linspace(-grid_lim, grid_lim, n)
    X, Y = np.meshgrid(xs, ys)
    points = np.stack([X, Y, np.zeros_like(X)], axis=-1)
    inv_cov = np.linalg.pinv(cov)
    q = np.einsum('...i,ij,...j->...', points, inv_cov, points)
    return X, Y, q


def pinhole_jacobian(point_cam, fx=1.7, fy=1.7):
    x, y, z = point_cam
    return np.array([
        [fx / z, 0.0, -fx * x / (z ** 2)],
        [0.0, fy / z, -fy * y / (z ** 2)],
    ], dtype=float)


def project_covariance(cov_world, mean_world, world_to_cam=None, fx=1.7, fy=1.7):
    if world_to_cam is None:
        world_to_cam = np.eye(3)
    mean_world = np.asarray(mean_world, dtype=float)
    mean_cam = world_to_cam @ mean_world
    cov_cam = world_to_cam @ cov_world @ world_to_cam.T
    J = pinhole_jacobian(mean_cam, fx=fx, fy=fy)
    cov_2d = J @ cov_cam @ J.T
    center = np.array([fx * mean_cam[0] / mean_cam[2], fy * mean_cam[1] / mean_cam[2]])
    return center, cov_2d, mean_cam, J


def ellipse_axes_and_angle(cov2d, radius_level=2.0):
    vals, vecs = np.linalg.eigh(cov2d)
    order = np.argsort(vals)[::-1]
    vals = vals[order]
    vecs = vecs[:, order]
    width = 2 * radius_level * np.sqrt(max(vals[0], 1e-9))
    height = 2 * radius_level * np.sqrt(max(vals[1], 1e-9))
    angle = np.degrees(np.arctan2(vecs[1, 0], vecs[0, 0]))
    return width, height, angle, vals


def add_covariance_ellipse(ax, center, cov2d, color, label, radius_level=2.0, alpha=0.24):
    width, height, angle, vals = ellipse_axes_and_angle(cov2d, radius_level=radius_level)
    ellipse = Ellipse(
        xy=center,
        width=width,
        height=height,
        angle=angle,
        facecolor=color,
        edgecolor=color,
        alpha=alpha,
        lw=2.2,
        label=label,
    )
    ax.add_patch(ellipse)
    ax.scatter([center[0]], [center[1]], color=color, s=32)
    return width, height, angle, vals


def rotation_matrix_2d(theta_deg):
    theta = np.deg2rad(theta_deg)
    c, s = np.cos(theta), np.sin(theta)
    return np.array([[c, -s], [s, c]], dtype=float)


def gaussian_2d_pdf(grid_x, grid_y, mean, cov):
    shifted = np.stack([grid_x - mean[0], grid_y - mean[1]], axis=-1)
    inv_cov = np.linalg.pinv(cov)
    q = np.einsum('...i,ij,...j->...', shifted, inv_cov, shifted)
    return np.exp(-0.5 * q)


def render_splats_2d(splats, n=220, extent=(-2.8, 2.8, -2.3, 2.3)):
    xs = np.linspace(extent[0], extent[1], n)
    ys = np.linspace(extent[2], extent[3], n)
    X, Y = np.meshgrid(xs, ys)
    rgb = np.zeros((n, n, 3), dtype=float)
    transmittance = np.ones((n, n), dtype=float)
    for splat in sorted(splats, key=lambda item: item['depth']):
        alpha = splat['opacity'] * gaussian_2d_pdf(X, Y, splat['mean'], splat['cov'])
        alpha = np.clip(alpha, 0.0, 0.98)
        color = np.asarray(splat['color'], dtype=float)[None, None, :]
        rgb += transmittance[..., None] * alpha[..., None] * color
        transmittance *= (1.0 - alpha)
    return X, Y, np.clip(rgb, 0.0, 1.0), 1.0 - transmittance


print('Geometry helpers ready.')

def rotation_matrix_3d(euler_deg): return Rotation3D.from_euler('xyz', euler_deg, degrees=True).as_matrix() def covariance_from_scale_rotation(scales, euler_deg=(0.0, 0.0, 0.0)): scales = np.asarray(scales, dtype=float) S = np.diag(scales) R = rotation_matrix_3d(euler_deg) return R @ S @ S.T @ R.T def unit_sphere_grid(nu=72, nv=36): u = np.linspace(0, 2 * np.pi, nu) v = np.linspace(0, np.pi, nv) x = np.outer(np.cos(u), np.sin(v)) y = np.outer(np.sin(u), np.sin(v)) z = np.outer(np.ones_like(u), np.cos(v)) return x, y, z def ellipsoid_surface(cov, mean=(0.0, 0.0, 0.0), radius_level=2.0): vals, vecs = np.linalg.eigh(cov) vals = np.maximum(vals, 1e-9) radii = radius_level * np.sqrt(vals) xs, ys, zs = unit_sphere_grid() sphere = np.stack([xs, ys, zs], axis=0).reshape(3, -1) ellipsoid = (vecs @ np.diag(radii) @ sphere).reshape(3, *xs.shape) ellipsoid[0] += mean[0] ellipsoid[1] += mean[1] ellipsoid[2] += mean[2] return ellipsoid def principal_axes(cov, mean=(0.0, 0.0, 0.0), radius_level=2.0): vals, vecs = np.linalg.eigh(cov) order = np.argsort(vals)[::-1] vals = vals[order] vecs = vecs[:, order] mean = np.asarray(mean, dtype=float) segments = [] for idx in range(3): direction = vecs[:, idx] * radius_level * np.sqrt(max(vals[idx], 1e-9)) segments.append((mean - direction, mean + direction)) return segments, vals, vecs def plot_ellipsoid(ax, cov, mean=(0.0, 0.0, 0.0), title='', color='#4C78A8', alpha=0.28): X, Y, Z = ellipsoid_surface(cov, mean=mean) ax.plot_surface(X, Y, Z, color=color, alpha=alpha, linewidth=0, shade=True) segments, vals, _ = principal_axes(cov, mean=mean) axis_colors = ['#E45756', '#72B7B2', '#F58518'] for (start, end), c in zip(segments, axis_colors): ax.plot([start[0], end[0]], [start[1], end[1]], [start[2], end[2]], color=c, lw=2.2) ax.scatter([mean[0]], [mean[1]], [mean[2]], color='black', s=28) ax.set_title(title) ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('z') ax.set_box_aspect([1, 1, 1]) ax.view_init(elev=24, azim=36) lim = 3.3 ax.set_xlim(-lim + mean[0], lim + mean[0]) ax.set_ylim(-lim + mean[1], lim + mean[1]) ax.set_zlim(-lim + mean[2], lim + mean[2]) return vals def quadratic_form_on_plane(cov, grid_lim=3.0, n=220): xs = np.linspace(-grid_lim, grid_lim, n) ys = np.linspace(-grid_lim, grid_lim, n) X, Y = np.meshgrid(xs, ys) points = np.stack([X, Y, np.zeros_like(X)], axis=-1) inv_cov = np.linalg.pinv(cov) q = np.einsum('...i,ij,...j->...', points, inv_cov, points) return X, Y, q def pinhole_jacobian(point_cam, fx=1.7, fy=1.7): x, y, z = point_cam return np.array([ [fx / z, 0.0, -fx * x / (z ** 2)], [0.0, fy / z, -fy * y / (z ** 2)], ], dtype=float) def project_covariance(cov_world, mean_world, world_to_cam=None, fx=1.7, fy=1.7): if world_to_cam is None: world_to_cam = np.eye(3) mean_world = np.asarray(mean_world, dtype=float) mean_cam = world_to_cam @ mean_world cov_cam = world_to_cam @ cov_world @ world_to_cam.T J = pinhole_jacobian(mean_cam, fx=fx, fy=fy) cov_2d = J @ cov_cam @ J.T center = np.array([fx * mean_cam[0] / mean_cam[2], fy * mean_cam[1] / mean_cam[2]]) return center, cov_2d, mean_cam, J def ellipse_axes_and_angle(cov2d, radius_level=2.0): vals, vecs = np.linalg.eigh(cov2d) order = np.argsort(vals)[::-1] vals = vals[order] vecs = vecs[:, order] width = 2 * radius_level * np.sqrt(max(vals[0], 1e-9)) height = 2 * radius_level * np.sqrt(max(vals[1], 1e-9)) angle = np.degrees(np.arctan2(vecs[1, 0], vecs[0, 0])) return width, height, angle, vals def add_covariance_ellipse(ax, center, cov2d, color, label, radius_level=2.0, alpha=0.24): width, height, angle, vals = ellipse_axes_and_angle(cov2d, radius_level=radius_level) ellipse = Ellipse( xy=center, width=width, height=height, angle=angle, facecolor=color, edgecolor=color, alpha=alpha, lw=2.2, label=label, ) ax.add_patch(ellipse) ax.scatter([center[0]], [center[1]], color=color, s=32) return width, height, angle, vals def rotation_matrix_2d(theta_deg): theta = np.deg2rad(theta_deg) c, s = np.cos(theta), np.sin(theta) return np.array([[c, -s], [s, c]], dtype=float) def gaussian_2d_pdf(grid_x, grid_y, mean, cov): shifted = np.stack([grid_x - mean[0], grid_y - mean[1]], axis=-1) inv_cov = np.linalg.pinv(cov) q = np.einsum('...i,ij,...j->...', shifted, inv_cov, shifted) return np.exp(-0.5 * q) def render_splats_2d(splats, n=220, extent=(-2.8, 2.8, -2.3, 2.3)): xs = np.linspace(extent[0], extent[1], n) ys = np.linspace(extent[2], extent[3], n) X, Y = np.meshgrid(xs, ys) rgb = np.zeros((n, n, 3), dtype=float) transmittance = np.ones((n, n), dtype=float) for splat in sorted(splats, key=lambda item: item['depth']): alpha = splat['opacity'] * gaussian_2d_pdf(X, Y, splat['mean'], splat['cov']) alpha = np.clip(alpha, 0.0, 0.98) color = np.asarray(splat['color'], dtype=float)[None, None, :] rgb += transmittance[..., None] * alpha[..., None] * color transmittance *= (1.0 - alpha) return X, Y, np.clip(rgb, 0.0, 1.0), 1.0 - transmittance print('Geometry helpers ready.')

Geometry helpers ready.

What you see¶

1. What the output shows: Only a short confirmation message. The real value is in the helper functions now loaded into memory.
2. Why it matters: Each later experiment reuses the same small mathematical pipeline: build $\Sigma$ from scale + rotation, inspect it in 3D, project it to 2D, then render splats by alpha compositing.
3. Key takeaway: We are studying one coherent object from multiple angles, not five disconnected plotting tricks.

---

Experiment 1: A Gaussian is a soft ellipsoid, not a fuzzy sphere¶

Question: What do mean, anisotropy, and orientation actually change in 3D space?

The 3DGS note says a primitive is not just a point with color: it is a Gaussian with a center $\mu$ and a covariance $\Sigma$. The center decides where the primitive lives; the covariance decides how it occupies space.

In [3]:

cov_sphere = covariance_from_scale_rotation([1.0, 1.0, 1.0], [0.0, 0.0, 0.0])
cov_flat = covariance_from_scale_rotation([1.8, 0.9, 0.28], [0.0, 0.0, 0.0])
cov_rot = covariance_from_scale_rotation([1.8, 0.9, 0.28], [30.0, 20.0, 55.0])

fig = plt.figure(figsize=(15.5, 4.8))
fig.suptitle('Experiment 1: mean + covariance determine the shape of a 3DGS primitive', fontsize=13, fontweight='bold')

ax1 = fig.add_subplot(1, 3, 1, projection='3d')
plot_ellipsoid(ax1, cov_sphere, mean=(0.0, 0.0, 0.0), title='Isotropic covariance\nall directions equally spread', color='#4C78A8')

ax2 = fig.add_subplot(1, 3, 2, projection='3d')
plot_ellipsoid(ax2, cov_flat, mean=(0.0, 0.0, 0.0), title='Anisotropic covariance\none thin axis, two broad axes', color='#54A24B')

ax3 = fig.add_subplot(1, 3, 3, projection='3d')
plot_ellipsoid(ax3, cov_rot, mean=(1.1, -0.8, 0.5), title='Same eigenvalues, new orientation + mean', color='#F58518')

plt.tight_layout()
plt.show()

cov_sphere = covariance_from_scale_rotation([1.0, 1.0, 1.0], [0.0, 0.0, 0.0]) cov_flat = covariance_from_scale_rotation([1.8, 0.9, 0.28], [0.0, 0.0, 0.0]) cov_rot = covariance_from_scale_rotation([1.8, 0.9, 0.28], [30.0, 20.0, 55.0]) fig = plt.figure(figsize=(15.5, 4.8)) fig.suptitle('Experiment 1: mean + covariance determine the shape of a 3DGS primitive', fontsize=13, fontweight='bold') ax1 = fig.add_subplot(1, 3, 1, projection='3d') plot_ellipsoid(ax1, cov_sphere, mean=(0.0, 0.0, 0.0), title='Isotropic covariance\nall directions equally spread', color='#4C78A8') ax2 = fig.add_subplot(1, 3, 2, projection='3d') plot_ellipsoid(ax2, cov_flat, mean=(0.0, 0.0, 0.0), title='Anisotropic covariance\none thin axis, two broad axes', color='#54A24B') ax3 = fig.add_subplot(1, 3, 3, projection='3d') plot_ellipsoid(ax3, cov_rot, mean=(1.1, -0.8, 0.5), title='Same eigenvalues, new orientation + mean', color='#F58518') plt.tight_layout() plt.show()

What you see¶

1. What the plot shows: The left primitive is a sphere-like Gaussian because all three directions have the same variance. The middle primitive becomes a flattened ellipsoid because one axis is much thinner than the others. The right primitive keeps the same axis lengths as the middle one, but rotating the covariance changes the ellipsoid's orientation, and changing the mean shifts the whole primitive in space.
2. Why it looks this way: A covariance matrix is not a blur amount attached to a point. Its eigenvalues set the spread along principal axes, and its eigenvectors set the orientation of those axes. That is why 3DGS can model local surface patches much better than isotropic point clouds.
3. Key takeaway: In 3DGS, the primitive is really a soft oriented ellipsoid; the covariance carries most of the local geometry.

---

Experiment 2: Why direct covariance is fragile¶

Question: Why doesn't 3DGS optimize the six symmetric entries of $\Sigma$ directly?

A true covariance must stay symmetric positive semidefinite. If gradient updates push a symmetric matrix outside that cone, the level sets stop being valid ellipsoids and the Gaussian interpretation breaks.

In [4]:

valid_cov = covariance_from_scale_rotation([1.7, 0.8, 0.35], [0.0, 0.0, 25.0])
bad_rot = rotation_matrix_3d([0.0, 0.0, 25.0])
invalid_cov = bad_rot @ np.diag([1.7**2, 0.8**2, -0.20]) @ bad_rot.T

Xv, Yv, q_valid = quadratic_form_on_plane(valid_cov)
Xi, Yi, q_invalid = quadratic_form_on_plane(invalid_cov)
eig_valid = np.linalg.eigvalsh(valid_cov)
eig_invalid = np.linalg.eigvalsh(invalid_cov)

fig, axes = plt.subplots(1, 3, figsize=(16.0, 4.7))
fig.suptitle('Experiment 2: one negative eigenvalue is enough to leave the Gaussian world', fontsize=13, fontweight='bold')

levels_valid = [0.5, 1.0, 2.0, 4.0, 7.0]
axes[0].contour(Xv, Yv, q_valid, levels=levels_valid, cmap='viridis')
axes[0].set_title('Valid SPD covariance\nclosed ellipses on the z=0 slice')
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].set_aspect('equal')

levels_invalid = [-6.0, -3.0, -1.0, 0.5, 1.0, 2.0, 4.0]
axes[1].contour(Xi, Yi, q_invalid, levels=levels_invalid, cmap='RdBu_r')
axes[1].set_title('Symmetric but invalid matrix\nopen saddle-like level sets')
axes[1].set_xlabel('x')
axes[1].set_ylabel('y')
axes[1].set_aspect('equal')

indices = np.arange(3)
bar_width = 0.36
axes[2].bar(indices - bar_width / 2, eig_valid, width=bar_width, color='#54A24B', label='valid SPD')
axes[2].bar(indices + bar_width / 2, eig_invalid, width=bar_width, color='#E45756', label='invalid symmetric')
axes[2].axhline(0.0, color='black', lw=1.2)
axes[2].set_xticks(indices)
axes[2].set_xticklabels(['$\\lambda_1$', '$\\lambda_2$', '$\\lambda_3$'])
axes[2].set_title('Eigenvalues of the two matrices')
axes[2].set_xlabel('eigenvalue index')
axes[2].set_ylabel('value')
axes[2].legend(frameon=False)

plt.tight_layout()
plt.show()

valid_cov = covariance_from_scale_rotation([1.7, 0.8, 0.35], [0.0, 0.0, 25.0]) bad_rot = rotation_matrix_3d([0.0, 0.0, 25.0]) invalid_cov = bad_rot @ np.diag([1.7**2, 0.8**2, -0.20]) @ bad_rot.T Xv, Yv, q_valid = quadratic_form_on_plane(valid_cov) Xi, Yi, q_invalid = quadratic_form_on_plane(invalid_cov) eig_valid = np.linalg.eigvalsh(valid_cov) eig_invalid = np.linalg.eigvalsh(invalid_cov) fig, axes = plt.subplots(1, 3, figsize=(16.0, 4.7)) fig.suptitle('Experiment 2: one negative eigenvalue is enough to leave the Gaussian world', fontsize=13, fontweight='bold') levels_valid = [0.5, 1.0, 2.0, 4.0, 7.0] axes[0].contour(Xv, Yv, q_valid, levels=levels_valid, cmap='viridis') axes[0].set_title('Valid SPD covariance\nclosed ellipses on the z=0 slice') axes[0].set_xlabel('x') axes[0].set_ylabel('y') axes[0].set_aspect('equal') levels_invalid = [-6.0, -3.0, -1.0, 0.5, 1.0, 2.0, 4.0] axes[1].contour(Xi, Yi, q_invalid, levels=levels_invalid, cmap='RdBu_r') axes[1].set_title('Symmetric but invalid matrix\nopen saddle-like level sets') axes[1].set_xlabel('x') axes[1].set_ylabel('y') axes[1].set_aspect('equal') indices = np.arange(3) bar_width = 0.36 axes[2].bar(indices - bar_width / 2, eig_valid, width=bar_width, color='#54A24B', label='valid SPD') axes[2].bar(indices + bar_width / 2, eig_invalid, width=bar_width, color='#E45756', label='invalid symmetric') axes[2].axhline(0.0, color='black', lw=1.2) axes[2].set_xticks(indices) axes[2].set_xticklabels(['$\\lambda_1$', '$\\lambda_2$', '$\\lambda_3$']) axes[2].set_title('Eigenvalues of the two matrices') axes[2].set_xlabel('eigenvalue index') axes[2].set_ylabel('value') axes[2].legend(frameon=False) plt.tight_layout() plt.show()

What you see¶

1. What the plot shows: The valid covariance produces neat closed ellipses on the plane, exactly what we expect from a Gaussian cross-section. The invalid symmetric matrix produces open, saddle-like contours because one eigenvalue has gone negative. The bar chart makes the failure explicit: symmetry alone is not enough.
2. Why it looks this way: A Gaussian needs $\mathbf{v}^T \Sigma \mathbf{v} \ge 0$ in every direction. Once an eigenvalue becomes negative, some directions get negative quadratic form values, so the object no longer behaves like an ellipsoid with non-negative variance. Directly optimizing the matrix entries makes this failure easy to trigger.
3. Key takeaway: 3DGS avoids optimizing raw covariance entries because it needs a parameterization that stays inside the valid SPD cone throughout training.

---

Experiment 3: The $R S S^T R^T$ decomposition separates size from orientation¶

Question: What changes when we alter scale only, versus rotation only?

This is the practical reason for the standard 3DGS parameterization. Scale controls the axis lengths; rotation controls only the coordinate frame of those axes.

In [5]:

cov_unit = covariance_from_scale_rotation([1.0, 1.0, 1.0], [0.0, 0.0, 0.0])
cov_scaled = covariance_from_scale_rotation([2.0, 1.0, 0.32], [0.0, 0.0, 0.0])
cov_scaled_rot = covariance_from_scale_rotation([2.0, 1.0, 0.32], [35.0, 30.0, 55.0])

fig = plt.figure(figsize=(15.5, 4.8))
fig.suptitle('Experiment 3: scale changes the radii; rotation changes the frame', fontsize=13, fontweight='bold')

ax1 = fig.add_subplot(1, 3, 1, projection='3d')
plot_ellipsoid(ax1, cov_unit, title='Step 1: unit sphere\nall scales equal', color='#4C78A8')

ax2 = fig.add_subplot(1, 3, 2, projection='3d')
plot_ellipsoid(ax2, cov_scaled, title='Step 2: apply scale only\nshape changes, frame stays aligned', color='#72B7B2')

ax3 = fig.add_subplot(1, 3, 3, projection='3d')
plot_ellipsoid(ax3, cov_scaled_rot, title='Step 3: rotate the scaled ellipsoid\nsame lengths, new orientation', color='#F58518')

plt.tight_layout()
plt.show()

cov_unit = covariance_from_scale_rotation([1.0, 1.0, 1.0], [0.0, 0.0, 0.0]) cov_scaled = covariance_from_scale_rotation([2.0, 1.0, 0.32], [0.0, 0.0, 0.0]) cov_scaled_rot = covariance_from_scale_rotation([2.0, 1.0, 0.32], [35.0, 30.0, 55.0]) fig = plt.figure(figsize=(15.5, 4.8)) fig.suptitle('Experiment 3: scale changes the radii; rotation changes the frame', fontsize=13, fontweight='bold') ax1 = fig.add_subplot(1, 3, 1, projection='3d') plot_ellipsoid(ax1, cov_unit, title='Step 1: unit sphere\nall scales equal', color='#4C78A8') ax2 = fig.add_subplot(1, 3, 2, projection='3d') plot_ellipsoid(ax2, cov_scaled, title='Step 2: apply scale only\nshape changes, frame stays aligned', color='#72B7B2') ax3 = fig.add_subplot(1, 3, 3, projection='3d') plot_ellipsoid(ax3, cov_scaled_rot, title='Step 3: rotate the scaled ellipsoid\nsame lengths, new orientation', color='#F58518') plt.tight_layout() plt.show()

What you see¶

1. What the plot shows: The middle ellipsoid differs from the sphere only by axis lengths: one axis is stretched, one stays moderate, and one becomes thin. The right ellipsoid has exactly the same radii as the middle one, but the whole frame has rotated in 3D.
2. Why it looks this way: In $\Sigma = R S S^T R^T$, the diagonal matrix $S$ sets the principal-axis magnitudes before any rotation happens. The rotation matrix $R$ then reorients that already-built ellipsoid. This division is much easier to optimize and interpret than six free covariance entries.
3. Key takeaway: The decomposition is not cosmetic. It is a clean geometric factorization: build the ellipsoid with scales first, then place it in the world with a rotation.

---

Experiment 4: Projection turns a 3D covariance into a 2D ellipse¶

Question: How does a 3D Gaussian become a screen-space footprint through $\Sigma' = J W \Sigma W^T J^T$?

In 3DGS, rendering never uses the raw 3D covariance directly. The renderer first projects the primitive into image space, where it becomes a 2D Gaussian ellipse.

In [6]:

front_cov = covariance_from_scale_rotation([1.15, 0.70, 0.24], [0.0, 0.0, 0.0])
tilted_cov = covariance_from_scale_rotation([1.15, 0.70, 0.24], [28.0, 18.0, 40.0])

center_near, cov2d_near, _, _ = project_covariance(front_cov, mean_world=[0.0, 0.0, 2.2])
center_far, cov2d_far, _, _ = project_covariance(front_cov, mean_world=[0.0, 0.0, 5.0])
center_front, cov2d_front, _, _ = project_covariance(front_cov, mean_world=[0.0, 0.0, 3.2])
center_tilt, cov2d_tilt, _, _ = project_covariance(tilted_cov, mean_world=[0.0, 0.0, 3.2])

near_w, near_h, _, _ = ellipse_axes_and_angle(cov2d_near)
far_w, far_h, _, _ = ellipse_axes_and_angle(cov2d_far)
front_w, front_h, _, _ = ellipse_axes_and_angle(cov2d_front)
tilt_w, tilt_h, _, _ = ellipse_axes_and_angle(cov2d_tilt)

fig, axes = plt.subplots(1, 3, figsize=(16.0, 4.9))
fig.suptitle('Experiment 4: the same 3D primitive can have very different 2D footprints', fontsize=13, fontweight='bold')

ax = axes[0]
add_covariance_ellipse(ax, center_near, cov2d_near, color='#4C78A8', label='near: z=2.2')
add_covariance_ellipse(ax, center_far, cov2d_far, color='#E45756', label='far: z=5.0')
ax.set_title('Depth effect on the projected ellipse')
ax.set_xlabel('screen x')
ax.set_ylabel('screen y')
ax.set_aspect('equal')
ax.set_xlim(-2.2, 2.2)
ax.set_ylim(-2.2, 2.2)
ax.legend(frameon=False, loc='upper right')

ax = axes[1]
add_covariance_ellipse(ax, center_front, cov2d_front, color='#72B7B2', label='front-facing')
add_covariance_ellipse(ax, center_tilt + np.array([0.18, -0.04]), cov2d_tilt, color='#F58518', label='tilted covariance')
ax.set_title('Orientation effect at the same depth')
ax.set_xlabel('screen x')
ax.set_ylabel('screen y')
ax.set_aspect('equal')
ax.set_xlim(-2.2, 2.2)
ax.set_ylim(-2.2, 2.2)
ax.legend(frameon=False, loc='upper right')

labels = ['near major', 'near minor', 'far major', 'far minor', 'front major', 'tilted major']
values = [near_w, near_h, far_w, far_h, front_w, tilt_w]
colors = ['#4C78A8', '#4C78A8', '#E45756', '#E45756', '#72B7B2', '#F58518']
axes[2].bar(np.arange(len(values)), values, color=colors)
axes[2].set_title('Screen-space axis lengths (2-sigma footprint)')
axes[2].set_xlabel('measured axis')
axes[2].set_ylabel('length in projected coordinates')
axes[2].set_xticks(np.arange(len(values)))
axes[2].set_xticklabels(labels, rotation=30, ha='right')

plt.tight_layout()
plt.show()

front_cov = covariance_from_scale_rotation([1.15, 0.70, 0.24], [0.0, 0.0, 0.0]) tilted_cov = covariance_from_scale_rotation([1.15, 0.70, 0.24], [28.0, 18.0, 40.0]) center_near, cov2d_near, _, _ = project_covariance(front_cov, mean_world=[0.0, 0.0, 2.2]) center_far, cov2d_far, _, _ = project_covariance(front_cov, mean_world=[0.0, 0.0, 5.0]) center_front, cov2d_front, _, _ = project_covariance(front_cov, mean_world=[0.0, 0.0, 3.2]) center_tilt, cov2d_tilt, _, _ = project_covariance(tilted_cov, mean_world=[0.0, 0.0, 3.2]) near_w, near_h, _, _ = ellipse_axes_and_angle(cov2d_near) far_w, far_h, _, _ = ellipse_axes_and_angle(cov2d_far) front_w, front_h, _, _ = ellipse_axes_and_angle(cov2d_front) tilt_w, tilt_h, _, _ = ellipse_axes_and_angle(cov2d_tilt) fig, axes = plt.subplots(1, 3, figsize=(16.0, 4.9)) fig.suptitle('Experiment 4: the same 3D primitive can have very different 2D footprints', fontsize=13, fontweight='bold') ax = axes[0] add_covariance_ellipse(ax, center_near, cov2d_near, color='#4C78A8', label='near: z=2.2') add_covariance_ellipse(ax, center_far, cov2d_far, color='#E45756', label='far: z=5.0') ax.set_title('Depth effect on the projected ellipse') ax.set_xlabel('screen x') ax.set_ylabel('screen y') ax.set_aspect('equal') ax.set_xlim(-2.2, 2.2) ax.set_ylim(-2.2, 2.2) ax.legend(frameon=False, loc='upper right') ax = axes[1] add_covariance_ellipse(ax, center_front, cov2d_front, color='#72B7B2', label='front-facing') add_covariance_ellipse(ax, center_tilt + np.array([0.18, -0.04]), cov2d_tilt, color='#F58518', label='tilted covariance') ax.set_title('Orientation effect at the same depth') ax.set_xlabel('screen x') ax.set_ylabel('screen y') ax.set_aspect('equal') ax.set_xlim(-2.2, 2.2) ax.set_ylim(-2.2, 2.2) ax.legend(frameon=False, loc='upper right') labels = ['near major', 'near minor', 'far major', 'far minor', 'front major', 'tilted major'] values = [near_w, near_h, far_w, far_h, front_w, tilt_w] colors = ['#4C78A8', '#4C78A8', '#E45756', '#E45756', '#72B7B2', '#F58518'] axes[2].bar(np.arange(len(values)), values, color=colors) axes[2].set_title('Screen-space axis lengths (2-sigma footprint)') axes[2].set_xlabel('measured axis') axes[2].set_ylabel('length in projected coordinates') axes[2].set_xticks(np.arange(len(values))) axes[2].set_xticklabels(labels, rotation=30, ha='right') plt.tight_layout() plt.show()

What you see¶

1. What the plot shows: The same 3D covariance produces a larger screen-space ellipse when it is closer to the camera and a smaller one when it is farther away. At a fixed depth, tilting the covariance changes the orientation and aspect ratio of the 2D footprint. The bar chart makes those differences quantitative.
2. Why it looks this way: The Jacobian $J$ of the pinhole camera rescales and mixes spatial directions depending on depth and viewing geometry. That is exactly what $\Sigma' = J W \Sigma W^T J^T$ encodes: the image-plane uncertainty is the camera-transformed 3D covariance seen through a local linearization of perspective projection.
3. Key takeaway: 3DGS renders projected ellipses, not raw 3D ellipsoids. Depth and orientation strongly change the final footprint even when the original 3D Gaussian stays the same.

---

Experiment 5: Alpha compositing and the main failure modes¶

Question: Why do covariance choice and depth ordering determine blur, streaks, and occlusion artifacts?

The screen-space ellipses are only half the story. The renderer still has to blend them front-to-back. This is where covariance mistakes become visible as image artifacts.

In [7]:

base_blue_cov = np.array([[0.18, 0.04], [0.04, 0.10]])
base_orange_cov = np.array([[0.22, -0.06], [-0.06, 0.12]])
needle_cov = rotation_matrix_2d(42.0) @ np.diag([0.42, 0.010]) @ rotation_matrix_2d(42.0).T

good_splats = [
    {'mean': (-0.85, 0.15), 'cov': base_blue_cov, 'opacity': 0.85, 'color': (0.20, 0.55, 0.95), 'depth': 1.8},
    {'mean': (0.55, -0.10), 'cov': base_orange_cov, 'opacity': 0.80, 'color': (0.95, 0.52, 0.18), 'depth': 2.4},
]
blur_splats = [
    {'mean': (-0.85, 0.15), 'cov': 2.8 * base_blue_cov, 'opacity': 0.85, 'color': (0.20, 0.55, 0.95), 'depth': 1.8},
    {'mean': (0.55, -0.10), 'cov': base_orange_cov, 'opacity': 0.80, 'color': (0.95, 0.52, 0.18), 'depth': 2.4},
]
needle_splats = [
    {'mean': (-0.85, 0.15), 'cov': needle_cov, 'opacity': 0.85, 'color': (0.20, 0.55, 0.95), 'depth': 1.8},
    {'mean': (0.55, -0.10), 'cov': base_orange_cov, 'opacity': 0.80, 'color': (0.95, 0.52, 0.18), 'depth': 2.4},
]
wrong_order_splats = [
    {'mean': (-0.85, 0.15), 'cov': base_blue_cov, 'opacity': 0.85, 'color': (0.20, 0.55, 0.95), 'depth': 2.6},
    {'mean': (0.55, -0.10), 'cov': base_orange_cov, 'opacity': 0.80, 'color': (0.95, 0.52, 0.18), 'depth': 1.7},
]

_, _, img_good, _ = render_splats_2d(good_splats)
_, _, img_blur, _ = render_splats_2d(blur_splats)
_, _, img_needle, _ = render_splats_2d(needle_splats)
_, _, img_wrong, _ = render_splats_2d(wrong_order_splats)

fig, axes = plt.subplots(2, 2, figsize=(11.6, 9.0))
fig.suptitle('Experiment 5: how covariance and depth errors become visible image artifacts', fontsize=13, fontweight='bold')

cases = [
    ('Reference: moderate covariances\ncorrect front-to-back order', img_good),
    ('Oversized covariance\nblur and over-coverage', img_blur),
    ('Needle-like covariance\nstreaky footprint', img_needle),
    ('Wrong depth order\nincorrect occlusion', img_wrong),
]

extent = (-2.8, 2.8, -2.3, 2.3)
for ax, (title, image) in zip(axes.ravel(), cases):
    ax.imshow(image, extent=extent, origin='lower')
    ax.set_title(title)
    ax.set_xlabel('screen x')
    ax.set_ylabel('screen y')
    ax.set_aspect('equal')
    ax.grid(False)

plt.tight_layout()
plt.show()

base_blue_cov = np.array([[0.18, 0.04], [0.04, 0.10]]) base_orange_cov = np.array([[0.22, -0.06], [-0.06, 0.12]]) needle_cov = rotation_matrix_2d(42.0) @ np.diag([0.42, 0.010]) @ rotation_matrix_2d(42.0).T good_splats = [ {'mean': (-0.85, 0.15), 'cov': base_blue_cov, 'opacity': 0.85, 'color': (0.20, 0.55, 0.95), 'depth': 1.8}, {'mean': (0.55, -0.10), 'cov': base_orange_cov, 'opacity': 0.80, 'color': (0.95, 0.52, 0.18), 'depth': 2.4}, ] blur_splats = [ {'mean': (-0.85, 0.15), 'cov': 2.8 * base_blue_cov, 'opacity': 0.85, 'color': (0.20, 0.55, 0.95), 'depth': 1.8}, {'mean': (0.55, -0.10), 'cov': base_orange_cov, 'opacity': 0.80, 'color': (0.95, 0.52, 0.18), 'depth': 2.4}, ] needle_splats = [ {'mean': (-0.85, 0.15), 'cov': needle_cov, 'opacity': 0.85, 'color': (0.20, 0.55, 0.95), 'depth': 1.8}, {'mean': (0.55, -0.10), 'cov': base_orange_cov, 'opacity': 0.80, 'color': (0.95, 0.52, 0.18), 'depth': 2.4}, ] wrong_order_splats = [ {'mean': (-0.85, 0.15), 'cov': base_blue_cov, 'opacity': 0.85, 'color': (0.20, 0.55, 0.95), 'depth': 2.6}, {'mean': (0.55, -0.10), 'cov': base_orange_cov, 'opacity': 0.80, 'color': (0.95, 0.52, 0.18), 'depth': 1.7}, ] _, _, img_good, _ = render_splats_2d(good_splats) _, _, img_blur, _ = render_splats_2d(blur_splats) _, _, img_needle, _ = render_splats_2d(needle_splats) _, _, img_wrong, _ = render_splats_2d(wrong_order_splats) fig, axes = plt.subplots(2, 2, figsize=(11.6, 9.0)) fig.suptitle('Experiment 5: how covariance and depth errors become visible image artifacts', fontsize=13, fontweight='bold') cases = [ ('Reference: moderate covariances\ncorrect front-to-back order', img_good), ('Oversized covariance\nblur and over-coverage', img_blur), ('Needle-like covariance\nstreaky footprint', img_needle), ('Wrong depth order\nincorrect occlusion', img_wrong), ] extent = (-2.8, 2.8, -2.3, 2.3) for ax, (title, image) in zip(axes.ravel(), cases): ax.imshow(image, extent=extent, origin='lower') ax.set_title(title) ax.set_xlabel('screen x') ax.set_ylabel('screen y') ax.set_aspect('equal') ax.grid(False) plt.tight_layout() plt.show()

What you see¶

1. What the plot shows: The reference image looks compact and plausible. Enlarging one covariance makes its splat wash over too much of the screen, producing blur. Making it extremely anisotropic turns it into a streak. Swapping the depth order changes which primitive should appear in front, so the overlap region looks wrong even though the covariances themselves are reasonable.
2. Why it looks this way: The renderer multiplies each footprint by opacity and blends it through front-to-back transmittance. A large covariance spreads alpha over too many pixels, a needle covariance concentrates alpha along a thin direction, and wrong ordering corrupts occlusion. These are direct image-space consequences of geometry choices made earlier in the pipeline.
3. Key takeaway: Covariance is not just an internal parameter. It directly controls the visible image footprint, while depth sorting decides whether the resulting splats compose into a believable scene.

---

Summary¶

Experiment	Theme	Key insight
1. A Gaussian is a soft ellipsoid	Primitive intuition	Mean places the primitive; covariance gives it shape and orientation
2. Why direct covariance is fragile	Validity	Symmetry is not enough; a Gaussian needs positive-semidefinite covariance
3. The $R S S^T R^T$ decomposition	Parameterization	Scale builds the ellipsoid, rotation places its frame
4. Projection to screen space	Rendering geometry	The 3D covariance becomes a 2D ellipse through the camera Jacobian
5. Alpha compositing and failure modes	Image artifacts	Blur, streaks, and wrong occlusion are screen-space consequences of covariance + depth decisions

Connection to the broader theory¶

The two source notes are really one story. 3D Gaussian Splatting introduces the big idea: replace a global implicit field with explicit Gaussian primitives that can be rasterized in real time. Covariance Matrix trong 3DGS explains why this works geometrically: each primitive is only useful because its covariance encodes an oriented local surface patch, stays valid through the $R S S^T R^T$ parameterization, and projects cleanly into a 2D ellipse.

For Latent-Anything, this notebook highlights the key design lesson: a Gaussian primitive is already a structured latent token. It has position, size, orientation, opacity, and appearance-facing consequences. That makes it far easier to inspect and manipulate than a single opaque vector or a monolithic MLP.