Lab
Visualizing Quadratic Forms in \(\mathbb{R}^2\) and \(\mathbb{R}^3\)
This lab explores quadratic forms and their geometric interpretations in 2D and 3D spaces. Quadratic forms are fundamental mathematical structures that appear in many contexts, including linear algebra, optimization, statistics, and computer graphics. In this lab, we will:
- Understand the connection between symmetric matrices and quadratic forms
- Visualize conic sections in \(\mathbb{R}^2\) as level sets of quadratic forms
- Apply the Spectral Theorem to analyze and transform quadratic forms
- Extend these concepts to quadratic surfaces in \(\mathbb{R}^3\)
This lab will help solidify concepts from our Linear Algebra course, particularly:
- Symmetric matrices and their properties
- Eigenvalues and eigenvectors
- The Spectral Theorem
- Linear transformations and change of basis
Required Libraries
import numpy as np
import matplotlib.pyplot as plt
import plotly.graph_objects as go
import sympy as sp
from IPython.display import display, MarkdownPart 1: Quadratic Forms in \(\mathbb{R}^2\)
A quadratic form in \(\mathbb{R}^2\) can be written as:
\(q(\mathbf{x}) = \mathbf{x}^T A \mathbf{x}\)
where \(\mathbf{x} = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}\) and \(A\) is a \(2 \times 2\) symmetric matrix.
When expanded, this gives:
\[q(x) = a_{11}x_1^2 + 2a_{12}x_1x_2 + a_{22}x_2^2\]
Make sure you understand the following helper functions:
- The first function creates a random rotation, which is easy in \(\mathbb{R}^2\) and harder in \(\mathbb{R}^3\), and it uses the eigenvalues and the spectral theorem to create a \(2\times 2\) symmetric matrix.
- The second function looks at the eigenvalues to determine the conic type
- The last function returns the quadratic forms in sympy for easy visualization
def create_symmetric_matrix_2d(eigenvalues, random_seed=42):
"""Create a 2×2 symmetric matrix with specified eigenvalues and random orientation"""
np.random.seed(random_seed)
# Create a random rotation matrix (for eigenvectors)
theta = np.random.uniform(0, 2*np.pi)
Q = np.array([
[np.cos(theta), -np.sin(theta)],
[np.sin(theta), np.cos(theta)]
])
# Create the diagonal matrix of eigenvalues
D = np.diag(eigenvalues)
# Apply the spectral theorem formula: A = QDQ^T
A = Q @ D @ Q.T
return A, Q, D, theta
def get_conic_type(eigenvalues):
"""Determine the type of conic section based on eigenvalue signs"""
pos_count = sum(1 for ev in eigenvalues if ev > 1e-10)
neg_count = sum(1 for ev in eigenvalues if ev < -1e-10)
zero_count = sum(1 for ev in eigenvalues if abs(ev) <= 1e-10)
if pos_count == 2:
return "ellipse"
elif pos_count == 1 and neg_count == 1:
return "hyperbola"
elif pos_count == 1 and zero_count == 1:
return "(degenerate as two parallel lines)"
elif neg_count == 2:
return "empty set (imaginary ellipse)"
else:
return "unknown"
def display_equation(A, c=1):
"""Display the equation of the conic in both original and canonical form"""
# Define symbolic variables
x, y = sp.symbols('x y')
# Extract matrix elements
a, b, d = A[0, 0], A[0, 1], A[1, 1]
# Create the quadratic form expression
quadratic_form = a*x**2 + 2*b*x*y + d*y**2
# Get eigenvalues for canonical form
eigenvalues = np.linalg.eigvalsh(A)
# Create canonical form expression
x_prime, y_prime = sp.symbols("x' y'")
canonical_form = eigenvalues[0] * x_prime**2 + eigenvalues[1] * y_prime**2
# Return forms for display
return (quadratic_form,canonical_form)Graph the conic sections
This function graphs the quadratic forms and you don’t need to understand what it does, but you can.
def plot_conic(A, c=1, title=None, grid_size=2, num_points=1000):
"""Plot the conic section \mathbf{x}^T A \mathbf{x} = c"""
# Extract matrix elements
a, b, d = A[0, 0], A[0, 1], A[1, 1]
# Get eigenvalues and eigenvectors to determine conic type
eigenvalues, eigenvectors = np.linalg.eigh(A)
conic_type = get_conic_type(eigenvalues)
# Create a grid of points
x = np.linspace(-grid_size, grid_size, num_points)
y = np.linspace(-grid_size, grid_size, num_points)
X, Y = np.meshgrid(x, y)
# Calculate the quadratic form
Z = a*X**2 + 2*b*X*Y + d*Y**2
# Create figure
fig, ax = plt.subplots(figsize=(10, 8))
# Plot the conic section
ax.contour(X, Y, Z, levels=[c], colors=['blue'], linewidths=3)
# Draw the eigenvectors (principal axes)
colors = ['red', 'green']
for i in range(2):
# Scale eigenvector for better visualization
scale = 1.0 / np.sqrt(abs(eigenvalues[i]) if abs(eigenvalues[i]) > 1e-10 else 1)
vector = eigenvectors[:, i] * scale
eigenvalue_sign = "+" if eigenvalues[i] > 0 else ("-" if eigenvalues[i] < 0 else "0")
ax.arrow(0, 0, vector[0], vector[1],
head_width=0.1, head_length=0.1,
fc=colors[i], ec=colors[i],
length_includes_head=True)
ax.text(vector[0]*1.1, vector[1]*1.1,
f"λ = {eigenvalues[i]:.2f} {eigenvalue_sign}",
color=colors[i])
# Add a marker at the origin
ax.plot(0, 0, 'ko', markersize=8)
# Set labels and title
ax.set_xlabel('X')
ax.set_ylabel('Y')
if not title:
title = f"{conic_type.replace('_', ' ').title()}"
ax.set_title(title)
# Set equal aspect ratio
ax.set_aspect('equal')
# Set grid and limits
ax.grid(True, alpha=0.3)
ax.set_xlim(-grid_size, grid_size)
ax.set_ylim(-grid_size, grid_size)
return fig, axExercise 1: Create Examples
The signs and values of the eigenvalues determine the type of conic graph:
| Eigenvalue Signs | Surface Type |
|---|---|
| All positive | Ellipse |
| One positive and one negative | Hyperbola |
| One positive and one zero | Two parallel lines |
| Two negatives | Empty set |
Change the eigenvalues to obtain different shapes
# Create examples by changing the eigenvalues
eigenvalues = [4, 1]
A, Q, D, theta = create_symmetric_matrix_2d(eigenvalues, random_seed=42)
print("Symmetric Matrix A:")
print(A)
print("\nEigenvectors Q (columns):")
print(Q)
print("\nEigenvalues D:")
print(D)
print(f"\nRotation angle: {np.degrees(theta):.2f} degrees")
# Verify spectral decomposition
print("\nVerify A = QDQ^T:")
print(np.round(Q @ D @ Q.T, 10))
# Plot
fig, ax = plot_conic(A)
# Display the equation
original, canonical = display_equation(A)
#
display(Markdown(f"""Original form: ${sp.latex(original)} = 1$"""))
display(Markdown(f"""Canonical form: ${sp.latex(canonical)} = 1$"""))Exercise 2:
- Use the functions to plot the quadratic function \(x^2+xy+y^2=1\). What is the rotation angle?
- Create a quadratic form in \(\mathbb{R}^2\) representing a hyperbola with axes aligned to the coordinate axes. How would you modify the matrix to rotate the hyperbola by 45 degrees counter clockwise?
Part 2: Quadratic Forms in \(\mathbb{R}^3\)
We can extend these concepts to 3D space, where quadratic forms create quadric surfaces.
A quadratic form in \(\mathbb{R}^3\) is: \(q(\mathbf{x}) = \mathbf{x}^T A \mathbf{x}\)
where \(\mathbf{x} = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix}\) and \(A\) is a \(3 \times 3\) symmetric matrix.
Let’s create similar helper functions for 3D.
- Notice how the random orthogonal matrix is created. First by creating a \(3\times 3\) random matrix with independent normally distributed values and then applying the QR decomposition.
- The second and third functions are similar to the ones in \(\mathbb{R}^2\).
- The fourth function simply applies the Spectral Theorem
# Helper functions
def create_symmetric_matrix_3d(eigenvalues, random_seed=42):
"""Create a 3×3 symmetric matrix with specified eigenvalues and random orientation"""
np.random.seed(random_seed)
# Create a random orthogonal matrix (for eigenvectors)
# Using QR decomposition of a random matrix
H = np.random.randn(3, 3)
Q, R = np.linalg.qr(H)
# Make sure Q is orthogonal (det(Q) = 1)
if np.linalg.det(Q) < 0:
Q[:, 0] = -Q[:, 0]
# Create the diagonal matrix of eigenvalues
D = np.diag(eigenvalues)
# Apply the spectral theorem formula: A = QDQ^T
A = Q @ D @ Q.T
return A, Q, D
def get_surface_type(eigenvalues):
"""Determine the type of quadratic surface based on eigenvalue signs"""
pos_count = sum(1 for ev in eigenvalues if ev > 1e-10)
neg_count = sum(1 for ev in eigenvalues if ev < -1e-10)
zero_count = sum(1 for ev in eigenvalues if abs(ev) <= 1e-10)
if pos_count == 3:
return "ellipsoid"
elif pos_count == 2 and neg_count == 1:
return "hyperboloid_one_sheet"
elif pos_count == 1 and neg_count == 2:
return "hyperboloid_two_sheets"
elif pos_count == 2 and zero_count == 1:
return "elliptic_cylinder"
elif pos_count == 1 and neg_count == 1 and zero_count == 1:
return "hyperbolic_cylinder"
elif pos_count == 1 and zero_count == 2:
return "parabolic_cylinder"
else:
return "unknown"
def display_3d_equation(A, c=1):
"""Display the equation of the quadric in both original and canonical form"""
# Define symbolic variables
x, y, z = sp.symbols('x y z')
# Extract matrix elements
a, b, c_val, d, e, f = A[0,0], A[0,1], A[0,2], A[1,1], A[1,2], A[2,2]
# Create the quadratic form expression
quadratic_form = a*x**2 + 2*b*x*y + 2*c_val*x*z + d*y**2 + 2*e*y*z + f*z**2
# Get eigenvalues for canonical form
eigenvalues = np.linalg.eigvalsh(A)
# Create canonical form expression
x_prime, y_prime, z_prime = sp.symbols("x' y' z'")
canonical_form = eigenvalues[0] * x_prime**2 + eigenvalues[1] * y_prime**2 + eigenvalues[2] * z_prime**2
# Display the equations
return (quadratic_form,canonical_form)
def apply_spectral_theorem(A):
"""Apply the spectral theorem to decompose the matrix"""
eigenvalues, eigenvectors = np.linalg.eigh(A)
return eigenvalues, eigenvectorsGraph the surfaces in plotly
These are complicated functions that we are just using.
def generate_surface_points(A, c=1, num_points=50, range_scale=2.0):
"""Generate points for the quadratic surface x^T A x = c"""
eigenvalues, eigenvectors = apply_spectral_theorem(A)
surface_type = get_surface_type(eigenvalues)
if surface_type == "ellipsoid":
# Parametrization in spherical coordinates
u = np.linspace(0, 2 * np.pi, num_points)
v = np.linspace(0, np.pi, num_points)
# Generate points on unit sphere
x_mesh = np.outer(np.cos(u), np.sin(v))
y_mesh = np.outer(np.sin(u), np.sin(v))
z_mesh = np.outer(np.ones_like(u), np.cos(v))
# Transform and scale
points = np.zeros((num_points, num_points, 3))
for i in range(num_points):
for j in range(num_points):
# Unit sphere point
point = np.array([x_mesh[i, j], y_mesh[i, j], z_mesh[i, j]])
# Scale by 1/sqrt(eigenvalues)
scaled_point = point / np.sqrt(np.abs(eigenvalues))
# Transform to original coordinates
points[i, j, :] = eigenvectors @ scaled_point * np.sqrt(c)
elif surface_type == "hyperboloid_one_sheet":
# For hyperboloid of one sheet: x²/a² + y²/b² - z²/c² = 1
u = np.linspace(0, 2 * np.pi, num_points)
v = np.linspace(-range_scale, range_scale, num_points)
U, V = np.meshgrid(u, v)
# Find which eigenvalue is negative
neg_index = np.argmin(eigenvalues)
pos_indices = [i for i in range(3) if i != neg_index]
# Create array mapping eigen-indices to parameter equations
param_array = np.zeros((num_points, num_points, 3))
# Set the two positive eigenvalue dimensions
param_array[:, :, pos_indices[0]] = np.sqrt(c / eigenvalues[pos_indices[0]]) * np.cos(U) * np.cosh(V)
param_array[:, :, pos_indices[1]] = np.sqrt(c / eigenvalues[pos_indices[1]]) * np.sin(U) * np.cosh(V)
# Set the negative eigenvalue dimension
param_array[:, :, neg_index] = np.sqrt(c / -eigenvalues[neg_index]) * np.sinh(V)
# Transform back to original coordinates
points = np.zeros((num_points, num_points, 3))
for i in range(num_points):
for j in range(num_points):
points[i, j, :] = eigenvectors @ param_array[i, j, :]
elif surface_type == "hyperboloid_two_sheets":
# For hyperboloid of two sheets: x²/a² - y²/b² - z²/c² = 1
u = np.linspace(0, 2 * np.pi, num_points)
v = np.linspace(0.1, range_scale, num_points) # Start slightly above 0 to avoid singularity
U, V = np.meshgrid(u, v)
# Find which eigenvalue is positive
pos_index = np.argmax(eigenvalues)
neg_indices = [i for i in range(3) if i != pos_index]
# Create points for both sheets
points = np.zeros((2, num_points, num_points, 3))
for sheet in range(2):
sign = 1 if sheet == 0 else -1
# Create array mapping eigen-indices to parameter equations
param_array = np.zeros((num_points, num_points, 3))
# Set the positive eigenvalue dimension (the axis of the hyperboloid)
param_array[:, :, pos_index] = sign * np.sqrt(c / eigenvalues[pos_index]) * np.cosh(V)
# Set the negative eigenvalue dimensions
param_array[:, :, neg_indices[0]] = np.sqrt(c / -eigenvalues[neg_indices[0]]) * np.cos(U) * np.sinh(V)
param_array[:, :, neg_indices[1]] = np.sqrt(c / -eigenvalues[neg_indices[1]]) * np.sin(U) * np.sinh(V)
# Transform back to original coordinates
for i in range(num_points):
for j in range(num_points):
points[sheet, i, j, :] = eigenvectors @ param_array[i, j, :]
elif surface_type == "elliptic_cylinder":
# Find which eigenvalue is zero (or close to zero)
zero_index = np.argmin(np.abs(eigenvalues))
nonzero_indices = [i for i in range(3) if i != zero_index]
# The cylinder extends along the eigenvector of the zero eigenvalue
# We'll parametrize as a cylinder with elliptical cross-section
u = np.linspace(0, 2 * np.pi, num_points)
v = np.linspace(-range_scale, range_scale, num_points)
U, V = np.meshgrid(u, v)
# Create array mapping eigen-indices to parameter equations
param_array = np.zeros((num_points, num_points, 3))
# Set the nonzero eigenvalue dimensions (the elliptical cross-section)
param_array[:, :, nonzero_indices[0]] = np.sqrt(c / eigenvalues[nonzero_indices[0]]) * np.cos(U)
param_array[:, :, nonzero_indices[1]] = np.sqrt(c / eigenvalues[nonzero_indices[1]]) * np.sin(U)
# Set the zero eigenvalue dimension (the axis of the cylinder)
param_array[:, :, zero_index] = V
# Transform back to original coordinates
points = np.zeros((num_points, num_points, 3))
for i in range(num_points):
for j in range(num_points):
points[i, j, :] = eigenvectors @ param_array[i, j, :]
elif surface_type == "hyperbolic_cylinder":
# Find indices of positive, negative, and zero eigenvalues
pos_index = np.argmax(eigenvalues)
neg_index = np.argmin(eigenvalues)
zero_index = 3 - pos_index - neg_index # The remaining index
if abs(eigenvalues[zero_index]) > 1e-10: # Check it's actually zero
zero_index = neg_index if abs(eigenvalues[neg_index]) < abs(eigenvalues[pos_index]) else pos_index
neg_index = 0 if zero_index != 0 and pos_index != 0 else 1 if zero_index != 1 and pos_index != 1 else 2
pos_index = 3 - zero_index - neg_index
# The cylinder extends along the eigenvector of the zero eigenvalue
u = np.linspace(-range_scale, range_scale, num_points)
v = np.linspace(-range_scale, range_scale, num_points)
U, V = np.meshgrid(u, v)
# Create points for both sheets of the hyperbola
points = np.zeros((2, num_points, num_points, 3))
for sheet in range(2):
sign = 1 if sheet == 0 else -1
# Create array mapping eigen-indices to parameter equations
param_array = np.zeros((num_points, num_points, 3))
# Set the positive eigenvalue dimension
param_array[:, :, pos_index] = sign * np.sqrt(c / eigenvalues[pos_index]) * np.cosh(U)
# Set the negative eigenvalue dimension
param_array[:, :, neg_index] = np.sqrt(c / -eigenvalues[neg_index]) * np.sinh(U)
# Set the zero eigenvalue dimension (the axis of the cylinder)
param_array[:, :, zero_index] = V
# Transform back to original coordinates
for i in range(num_points):
for j in range(num_points):
points[sheet, i, j, :] = eigenvectors @ param_array[i, j, :]
else:
# For other surfaces or unknown types, return a placeholder
points = None
return points, eigenvalues, eigenvectors, surface_type
def plot_quadratic_surface(A, c=1, title=None):
"""Plot the quadratic surface x^T A x = c"""
points, eigenvalues, eigenvectors, surface_type = generate_surface_points(A, c)
if points is None:
print(f"Surface type '{surface_type}' visualization not implemented")
return None
fig = go.Figure()
if surface_type == "ellipsoid":
x = points[:, :, 0]
y = points[:, :, 1]
z = points[:, :, 2]
fig.add_trace(go.Surface(
x=x, y=y, z=z,
colorscale='Viridis',
opacity=0.8,
showscale=False
))
elif surface_type == "hyperboloid_one_sheet":
x = points[:, :, 0]
y = points[:, :, 1]
z = points[:, :, 2]
fig.add_trace(go.Surface(
x=x, y=y, z=z,
colorscale='Plasma',
opacity=0.8,
showscale=False
))
elif surface_type == "hyperboloid_two_sheets":
# Add both sheets
for sheet in range(2):
x = points[sheet, :, :, 0]
y = points[sheet, :, :, 1]
z = points[sheet, :, :, 2]
fig.add_trace(go.Surface(
x=x, y=y, z=z,
colorscale='Cividis',
opacity=0.8,
showscale=False
))
elif surface_type == "elliptic_cylinder":
x = points[:, :, 0]
y = points[:, :, 1]
z = points[:, :, 2]
fig.add_trace(go.Surface(
x=x, y=y, z=z,
colorscale='Blues',
opacity=0.8,
showscale=False
))
elif surface_type == "hyperbolic_cylinder":
# Add both sheets
for sheet in range(2):
x = points[sheet, :, :, 0]
y = points[sheet, :, :, 1]
z = points[sheet, :, :, 2]
fig.add_trace(go.Surface(
x=x, y=y, z=z,
colorscale='Reds',
opacity=0.8,
showscale=False
))
# Add principal axes (eigenvectors)
colors = ['red', 'green', 'blue']
for i in range(3):
# Scale eigenvector for visualization
scale = 2.5 / max(abs(e) for e in eigenvalues if abs(e) > 1e-10)
vector = eigenvectors[:, i] * scale
eigenvalue_sign = "+" if eigenvalues[i] > 1e-10 else ("-" if eigenvalues[i] < -1e-10 else "0")
fig.add_trace(go.Scatter3d(
x=[0, vector[0]], y=[0, vector[1]], z=[0, vector[2]],
mode='lines+markers',
line=dict(color=colors[i], width=6),
marker=dict(size=5, color=colors[i]),
name=f"Axis {i+1} (λ = {eigenvalues[i]:.2f} {eigenvalue_sign})"
))
# Add a marker at the origin
fig.add_trace(go.Scatter3d(
x=[0], y=[0], z=[0],
mode='markers',
marker=dict(size=8, color='black'),
name="Origin"
))
if not title:
title = f"{surface_type.replace('_', ' ').title()}"
# Set plot layout
fig.update_layout(
title=title,
scene=dict(
xaxis_title='X',
yaxis_title='Y',
zaxis_title='Z',
aspectmode='data'
),
width=800,
height=800,
showlegend=True
)
return figExercise 3: Create Examples
Types of Quadratic Surfaces
The signs and values of the eigenvalues determine the type of quadratic surface:
| Eigenvalue Signs | Surface Type |
|---|---|
| All positive | Ellipsoid |
| Two positive, one negative | Hyperboloid of one sheet |
| One positive, two negative | Hyperboloid of two sheets |
| Two positive, one zero | Elliptic cylinder |
| One positive, one negative, one zero | Hyperbolic cylinder |
| One positive, two zero | Two parallel planes |
Change the eigenvalues to generate the different objects
# Create
eigenvalues = [1 , 2, 3]
A, Q, D = create_symmetric_matrix_3d(eigenvalues)
print("Symmetric Matrix A:")
print(A)
print("\nEigenvectors Q (columns):")
print(Q)
print("\nEigenvalues D:")
print(D)
# Verify spectral decomposition
print("\nVerify A = QDQ^T:")
print(np.round(Q @ D @ Q.T, 10))
# Display the equation
original, canonical = display_3d_equation(A)
#
display(Markdown(f"""Original form: ${sp.latex(original)} = 1$"""))
display(Markdown(f"""Canonical form: ${sp.latex(canonical)} = 1$"""))
# Plot the ellipsoid using Plotly for interactive visualization
#fig = plot_quadric_surface(A)
fig = plot_quadratic_surface(A)
fig.show() # In Jupyter, this displays an interactive 3D plotExercise 4
- Use the functions to plot the quadratic surface \[2x^2+xy+xz+3y^2-yz-z^2=1\]