Theory Notebook

Converted from theory.ipynb for web reading.

Optimization on Manifolds

Optimization on manifolds turns curved constraints into geometry-aware algorithms for spheres, Stiefel and Grassmann manifolds, SPD matrices, low-rank models, and natural gradients.

This notebook is the executable companion to notes.md. It uses spheres, tangent projections, geodesic interpolation, SPD matrices, and orthogonality constraints as concrete geometry laboratories.

Code cell 2

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

try:
    import seaborn as sns
    sns.set_theme(style="whitegrid", palette="colorblind")
    HAS_SNS = True
except ImportError:
    plt.style.use("seaborn-v0_8-whitegrid")
    HAS_SNS = False

mpl.rcParams.update({
    "figure.figsize":    (10, 6),
    "figure.dpi":         120,
    "font.size":           13,
    "axes.titlesize":      15,
    "axes.labelsize":      13,
    "xtick.labelsize":     11,
    "ytick.labelsize":     11,
    "legend.fontsize":     11,
    "legend.framealpha":   0.85,
    "lines.linewidth":      2.0,
    "axes.spines.top":     False,
    "axes.spines.right":   False,
    "savefig.bbox":       "tight",
    "savefig.dpi":         150,
})
np.random.seed(42)
print("Plot setup complete.")

Code cell 3

COLORS = {
    "primary":   "#0077BB",
    "secondary": "#EE7733",
    "tertiary":  "#009988",
    "error":     "#CC3311",
    "neutral":   "#555555",
    "highlight": "#EE3377",
}

def header(title):
    print("\n" + "=" * 72)
    print(title)
    print("=" * 72)

def check_true(condition, name):
    ok = bool(condition)
    print(f"{'PASS' if ok else 'FAIL'} - {name}")
    assert ok, name

def check_close(value, target, tol=1e-8, name="value"):
    ok = abs(float(value) - float(target)) <= tol
    print(f"{'PASS' if ok else 'FAIL'} - {name}: got {float(value):.6f}, expected {float(target):.6f}")
    assert ok, name

def normalize(x):
    x = np.asarray(x, dtype=float)
    n = np.linalg.norm(x)
    if n == 0:
        raise ValueError("cannot normalize zero vector")
    return x / n

def tangent_projection_sphere(x, v):
    x = normalize(x)
    v = np.asarray(v, dtype=float)
    return v - np.dot(x, v) * x

def exp_sphere(x, v):
    x = normalize(x)
    v = tangent_projection_sphere(x, v)
    n = np.linalg.norm(v)
    if n < 1e-12:
        return x.copy()
    return np.cos(n) * x + np.sin(n) * v / n

def retract_sphere(x, v):
    return normalize(np.asarray(x, dtype=float) + np.asarray(v, dtype=float))

def slerp(x, y, t):
    x = normalize(x)
    y = normalize(y)
    dot = np.clip(np.dot(x, y), -1.0, 1.0)
    theta = np.arccos(dot)
    if theta < 1e-12:
        return x.copy()
    return (np.sin((1 - t) * theta) * x + np.sin(t * theta) * y) / np.sin(theta)

def riemannian_gradient_sphere(x, euclidean_grad):
    return tangent_projection_sphere(x, euclidean_grad)

def sphere_descent(target, steps=20, eta=0.3):
    target = normalize(target)
    x = normalize(np.array([1.0, 0.2, 0.1]))
    values = []
    for _ in range(steps):
        values.append(float(-np.dot(target, x)))
        egrad = -target
        rgrad = riemannian_gradient_sphere(x, egrad)
        x = retract_sphere(x, -eta * rgrad)
    values.append(float(-np.dot(target, x)))
    return x, np.array(values)

def stiefel_tangent_projection(Q, Z):
    Q = np.asarray(Q, dtype=float)
    Z = np.asarray(Z, dtype=float)
    sym = 0.5 * (Q.T @ Z + Z.T @ Q)
    return Z - Q @ sym

def qr_retraction(Y):
    Q, R = np.linalg.qr(Y)
    signs = np.sign(np.diag(R))
    signs[signs == 0] = 1.0
    return Q @ np.diag(signs)

def spd_from_eigs(eigs):
    Q = np.array([[np.cos(0.4), -np.sin(0.4)], [np.sin(0.4), np.cos(0.4)]])
    return Q @ np.diag(eigs) @ Q.T

def mat_log_spd(A):
    vals, vecs = np.linalg.eigh(A)
    return vecs @ np.diag(np.log(vals)) @ vecs.T

def mat_invsqrt_spd(A):
    vals, vecs = np.linalg.eigh(A)
    return vecs @ np.diag(1.0 / np.sqrt(vals)) @ vecs.T

def spd_distance(A, B):
    C = mat_invsqrt_spd(A) @ B @ mat_invsqrt_spd(A)
    return float(np.linalg.norm(mat_log_spd(C), ord="fro"))

print("Differential-geometry helpers ready.")

Demo 1: Constrained optimization as unconstrained manifold optimization

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 5

header("Demo 1 - Constrained optimization as unconstrained manifold optimization: sphere tangent projection")
x = normalize(np.array([1.0, 1.0, 1.0]))
v = np.array([2.0, -1.0, 0.5])
tangent = tangent_projection_sphere(x, v)
print("Point on sphere:", np.round(x, 3).tolist())
print("Projected tangent:", np.round(tangent, 3).tolist())
check_close(np.dot(x, tangent), 0.0, tol=1e-10, name="orthogonal to base point")

Demo 2: Why projection alone is not always geometry-aware

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 7

header("Demo 2 - Why projection alone is not always geometry-aware: local chart for circle")
theta = np.linspace(-1.2, 1.2, 200)
circle = np.column_stack([np.cos(theta), np.sin(theta)])
fig, ax = plt.subplots(figsize=(6, 6))
ax.plot(circle[:, 0], circle[:, 1], color=COLORS["primary"], label="chart image")
ax.scatter([1], [0], color=COLORS["highlight"], label="base point")
ax.set_title("Local coordinate patch on the circle")
ax.set_xlabel("$x_1$")
ax.set_ylabel("$x_2$")
ax.axis("equal")
ax.legend()
fig.tight_layout()
plt.show()
plt.close(fig)
print("Chart parameter range:", (float(theta.min()), float(theta.max())))
check_true(circle.shape == (200, 2), "chart maps one coordinate to ambient R^2")

Demo 3: Tangent-space updates and retractions

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 9

header("Demo 3 - Tangent-space updates and retractions: spherical geodesic interpolation")
x = normalize(np.array([1.0, 0.0, 0.0]))
y = normalize(np.array([0.0, 1.0, 0.0]))
pts = np.array([slerp(x, y, t) for t in np.linspace(0, 1, 25)])
norms = np.linalg.norm(pts, axis=1)
print("First point:", pts[0].round(3).tolist())
print("Middle point:", pts[len(pts)//2].round(3).tolist())
check_true(np.all(np.abs(norms - 1.0) < 1e-10), "slerp stays on sphere")
fig, ax = plt.subplots(figsize=(6, 6))
ax.plot(pts[:, 0], pts[:, 1], color=COLORS["secondary"], marker="o", label="geodesic arc")
ax.set_title("Great-circle interpolation")
ax.set_xlabel("$x_1$")
ax.set_ylabel("$x_2$")
ax.axis("equal")
ax.legend()
fig.tight_layout()
plt.show()
plt.close(fig)

Demo 4: Matrix manifolds in ML

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 11

header("Demo 4 - Matrix manifolds in ML: Riemannian gradient on sphere")
x = normalize(np.array([0.7, 0.2, 0.6]))
target = normalize(np.array([0.0, 1.0, 0.0]))
egrad = -target
rgrad = riemannian_gradient_sphere(x, egrad)
print("Euclidean gradient:", np.round(egrad, 3).tolist())
print("Riemannian gradient:", np.round(rgrad, 3).tolist())
check_close(np.dot(x, rgrad), 0.0, tol=1e-10, name="Riemannian gradient is tangent")
print("Interpretation: steepest descent must live in the tangent space.")

Demo 5: Optimization pipeline: gradient step retract transport

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 13

header("Demo 5 - Optimization pipeline: gradient step retract transport: sphere gradient descent")
target = np.array([0.0, 1.0, 0.0])
x_final, values = sphere_descent(target, steps=25, eta=0.25)
print("Final point:", np.round(x_final, 3).tolist())
print("Initial objective:", round(values[0], 4))
print("Final objective:", round(values[-1], 4))
check_true(values[-1] < values[0], "objective decreases on sphere")
fig, ax = plt.subplots()
ax.plot(values, color=COLORS["primary"], label="objective")
ax.set_title("Riemannian gradient descent on the sphere")
ax.set_xlabel("Step")
ax.set_ylabel("$-a^T x$")
ax.legend()
fig.tight_layout()
plt.show()
plt.close(fig)

Demo 6: Riemannian gradient descent update

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 15

header("Demo 6 - Riemannian gradient descent update: Stiefel tangent projection")
Q = np.eye(3, 2)
Z = np.array([[0.2, 1.0], [1.5, -0.3], [0.7, 0.4]])
Xi = stiefel_tangent_projection(Q, Z)
constraint = Q.T @ Xi + Xi.T @ Q
print("Tangent constraint matrix:", np.round(constraint, 10))
check_true(np.linalg.norm(constraint) < 1e-10, "Stiefel tangent condition holds")
Y = qr_retraction(Q + 0.2 * Xi)
check_true(np.linalg.norm(Y.T @ Y - np.eye(2)) < 1e-10, "QR retraction returns orthonormal columns")
print("Retracted columns are orthonormal.")

Demo 7: Retraction $R_p(\mathbf{v})$

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 17

header("Demo 7 - Retraction $R_p(\\mathbf{v})$: SPD manifold distance")
A = spd_from_eigs([1.0, 3.0])
B = spd_from_eigs([2.0, 5.0])
d_ab = spd_distance(A, B)
d_ba = spd_distance(B, A)
print("Distance A to B:", round(d_ab, 6))
print("Distance B to A:", round(d_ba, 6))
check_close(d_ab, d_ba, tol=1e-10, name="SPD distance symmetry")
check_true(d_ab > 0, "distinct SPD matrices have positive distance")

Demo 8: Vector transport preview

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 19

header("Demo 8 - Vector transport preview: exponential vs retraction")
x = normalize(np.array([1.0, 0.0, 0.0]))
v = tangent_projection_sphere(x, np.array([0.0, 0.2, 0.1]))
exp_point = exp_sphere(x, v)
ret_point = retract_sphere(x, v)
print("Exponential point:", np.round(exp_point, 5).tolist())
print("Retraction point:", np.round(ret_point, 5).tolist())
check_true(abs(np.linalg.norm(exp_point) - 1.0) < 1e-10, "exponential stays on sphere")
check_true(abs(np.linalg.norm(ret_point) - 1.0) < 1e-10, "retraction stays on sphere")
print("Distance between exp and retraction:", round(float(np.linalg.norm(exp_point - ret_point)), 6))

Demo 9: Riemannian Hessian preview

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 21

header("Demo 9 - Riemannian Hessian preview: sphere tangent projection")
x = normalize(np.array([1.0, 1.0, 1.0]))
v = np.array([2.0, -1.0, 0.5])
tangent = tangent_projection_sphere(x, v)
print("Point on sphere:", np.round(x, 3).tolist())
print("Projected tangent:", np.round(tangent, 3).tolist())
check_close(np.dot(x, tangent), 0.0, tol=1e-10, name="orthogonal to base point")

Demo 10: First-order optimality on manifolds

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 23

header("Demo 10 - First-order optimality on manifolds: local chart for circle")
theta = np.linspace(-1.2, 1.2, 200)
circle = np.column_stack([np.cos(theta), np.sin(theta)])
fig, ax = plt.subplots(figsize=(6, 6))
ax.plot(circle[:, 0], circle[:, 1], color=COLORS["primary"], label="chart image")
ax.scatter([1], [0], color=COLORS["highlight"], label="base point")
ax.set_title("Local coordinate patch on the circle")
ax.set_xlabel("$x_1$")
ax.set_ylabel("$x_2$")
ax.axis("equal")
ax.legend()
fig.tight_layout()
plt.show()
plt.close(fig)
print("Chart parameter range:", (float(theta.min()), float(theta.max())))
check_true(circle.shape == (200, 2), "chart maps one coordinate to ambient R^2")

Demo 11: Sphere optimization

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 25

header("Demo 11 - Sphere optimization: spherical geodesic interpolation")
x = normalize(np.array([1.0, 0.0, 0.0]))
y = normalize(np.array([0.0, 1.0, 0.0]))
pts = np.array([slerp(x, y, t) for t in np.linspace(0, 1, 25)])
norms = np.linalg.norm(pts, axis=1)
print("First point:", pts[0].round(3).tolist())
print("Middle point:", pts[len(pts)//2].round(3).tolist())
check_true(np.all(np.abs(norms - 1.0) < 1e-10), "slerp stays on sphere")
fig, ax = plt.subplots(figsize=(6, 6))
ax.plot(pts[:, 0], pts[:, 1], color=COLORS["secondary"], marker="o", label="geodesic arc")
ax.set_title("Great-circle interpolation")
ax.set_xlabel("$x_1$")
ax.set_ylabel("$x_2$")
ax.axis("equal")
ax.legend()
fig.tight_layout()
plt.show()
plt.close(fig)

Demo 12: Stiefel manifold and orthogonality constraints

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 27

header("Demo 12 - Stiefel manifold and orthogonality constraints: Riemannian gradient on sphere")
x = normalize(np.array([0.7, 0.2, 0.6]))
target = normalize(np.array([0.0, 1.0, 0.0]))
egrad = -target
rgrad = riemannian_gradient_sphere(x, egrad)
print("Euclidean gradient:", np.round(egrad, 3).tolist())
print("Riemannian gradient:", np.round(rgrad, 3).tolist())
check_close(np.dot(x, rgrad), 0.0, tol=1e-10, name="Riemannian gradient is tangent")
print("Interpretation: steepest descent must live in the tangent space.")

Demo 13: Grassmann manifold and subspace learning

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 29

header("Demo 13 - Grassmann manifold and subspace learning: sphere gradient descent")
target = np.array([0.0, 1.0, 0.0])
x_final, values = sphere_descent(target, steps=25, eta=0.25)
print("Final point:", np.round(x_final, 3).tolist())
print("Initial objective:", round(values[0], 4))
print("Final objective:", round(values[-1], 4))
check_true(values[-1] < values[0], "objective decreases on sphere")
fig, ax = plt.subplots()
ax.plot(values, color=COLORS["primary"], label="objective")
ax.set_title("Riemannian gradient descent on the sphere")
ax.set_xlabel("Step")
ax.set_ylabel("$-a^T x$")
ax.legend()
fig.tight_layout()
plt.show()
plt.close(fig)

Demo 14: SPD manifold and covariance learning

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 31

header("Demo 14 - SPD manifold and covariance learning: Stiefel tangent projection")
Q = np.eye(3, 2)
Z = np.array([[0.2, 1.0], [1.5, -0.3], [0.7, 0.4]])
Xi = stiefel_tangent_projection(Q, Z)
constraint = Q.T @ Xi + Xi.T @ Q
print("Tangent constraint matrix:", np.round(constraint, 10))
check_true(np.linalg.norm(constraint) < 1e-10, "Stiefel tangent condition holds")
Y = qr_retraction(Q + 0.2 * Xi)
check_true(np.linalg.norm(Y.T @ Y - np.eye(2)) < 1e-10, "QR retraction returns orthonormal columns")
print("Retracted columns are orthonormal.")

Demo 15: Riemannian line search and trust-region preview

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 33

header("Demo 15 - Riemannian line search and trust-region preview: SPD manifold distance")
A = spd_from_eigs([1.0, 3.0])
B = spd_from_eigs([2.0, 5.0])
d_ab = spd_distance(A, B)
d_ba = spd_distance(B, A)
print("Distance A to B:", round(d_ab, 6))
print("Distance B to A:", round(d_ba, 6))
check_close(d_ab, d_ba, tol=1e-10, name="SPD distance symmetry")
check_true(d_ab > 0, "distinct SPD matrices have positive distance")

Demo 16: Orthogonal weight constraints

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 35

header("Demo 16 - Orthogonal weight constraints: exponential vs retraction")
x = normalize(np.array([1.0, 0.0, 0.0]))
v = tangent_projection_sphere(x, np.array([0.0, 0.2, 0.1]))
exp_point = exp_sphere(x, v)
ret_point = retract_sphere(x, v)
print("Exponential point:", np.round(exp_point, 5).tolist())
print("Retraction point:", np.round(ret_point, 5).tolist())
check_true(abs(np.linalg.norm(exp_point) - 1.0) < 1e-10, "exponential stays on sphere")
check_true(abs(np.linalg.norm(ret_point) - 1.0) < 1e-10, "retraction stays on sphere")
print("Distance between exp and retraction:", round(float(np.linalg.norm(exp_point - ret_point)), 6))

Demo 17: Low-rank matrix factorization

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 37

header("Demo 17 - Low-rank matrix factorization: sphere tangent projection")
x = normalize(np.array([1.0, 1.0, 1.0]))
v = np.array([2.0, -1.0, 0.5])
tangent = tangent_projection_sphere(x, v)
print("Point on sphere:", np.round(x, 3).tolist())
print("Projected tangent:", np.round(tangent, 3).tolist())
check_close(np.dot(x, tangent), 0.0, tol=1e-10, name="orthogonal to base point")

Demo 18: PCA as Grassmann optimization

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 39

header("Demo 18 - PCA as Grassmann optimization: local chart for circle")
theta = np.linspace(-1.2, 1.2, 200)
circle = np.column_stack([np.cos(theta), np.sin(theta)])
fig, ax = plt.subplots(figsize=(6, 6))
ax.plot(circle[:, 0], circle[:, 1], color=COLORS["primary"], label="chart image")
ax.scatter([1], [0], color=COLORS["highlight"], label="base point")
ax.set_title("Local coordinate patch on the circle")
ax.set_xlabel("$x_1$")
ax.set_ylabel("$x_2$")
ax.axis("equal")
ax.legend()
fig.tight_layout()
plt.show()
plt.close(fig)
print("Chart parameter range:", (float(theta.min()), float(theta.max())))
check_true(circle.shape == (200, 2), "chart maps one coordinate to ambient R^2")

Demo 19: Natural-gradient learning and Fisher geometry

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 41

header("Demo 19 - Natural-gradient learning and Fisher geometry: spherical geodesic interpolation")
x = normalize(np.array([1.0, 0.0, 0.0]))
y = normalize(np.array([0.0, 1.0, 0.0]))
pts = np.array([slerp(x, y, t) for t in np.linspace(0, 1, 25)])
norms = np.linalg.norm(pts, axis=1)
print("First point:", pts[0].round(3).tolist())
print("Middle point:", pts[len(pts)//2].round(3).tolist())
check_true(np.all(np.abs(norms - 1.0) < 1e-10), "slerp stays on sphere")
fig, ax = plt.subplots(figsize=(6, 6))
ax.plot(pts[:, 0], pts[:, 1], color=COLORS["secondary"], marker="o", label="geodesic arc")
ax.set_title("Great-circle interpolation")
ax.set_xlabel("$x_1$")
ax.set_ylabel("$x_2$")
ax.axis("equal")
ax.legend()
fig.tight_layout()
plt.show()
plt.close(fig)

Demo 20: Hyperbolic and SPD representation learning

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 43

header("Demo 20 - Hyperbolic and SPD representation learning: Riemannian gradient on sphere")
x = normalize(np.array([0.7, 0.2, 0.6]))
target = normalize(np.array([0.0, 1.0, 0.0]))
egrad = -target
rgrad = riemannian_gradient_sphere(x, egrad)
print("Euclidean gradient:", np.round(egrad, 3).tolist())
print("Riemannian gradient:", np.round(rgrad, 3).tolist())
check_close(np.dot(x, rgrad), 0.0, tol=1e-10, name="Riemannian gradient is tangent")
print("Interpretation: steepest descent must live in the tangent space.")

Demo 21: Constrained optimization as unconstrained manifold optimization

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 45

header("Demo 21 - Constrained optimization as unconstrained manifold optimization: sphere gradient descent")
target = np.array([0.0, 1.0, 0.0])
x_final, values = sphere_descent(target, steps=25, eta=0.25)
print("Final point:", np.round(x_final, 3).tolist())
print("Initial objective:", round(values[0], 4))
print("Final objective:", round(values[-1], 4))
check_true(values[-1] < values[0], "objective decreases on sphere")
fig, ax = plt.subplots()
ax.plot(values, color=COLORS["primary"], label="objective")
ax.set_title("Riemannian gradient descent on the sphere")
ax.set_xlabel("Step")
ax.set_ylabel("$-a^T x$")
ax.legend()
fig.tight_layout()
plt.show()
plt.close(fig)

Demo 22: Why projection alone is not always geometry-aware

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 47

header("Demo 22 - Why projection alone is not always geometry-aware: Stiefel tangent projection")
Q = np.eye(3, 2)
Z = np.array([[0.2, 1.0], [1.5, -0.3], [0.7, 0.4]])
Xi = stiefel_tangent_projection(Q, Z)
constraint = Q.T @ Xi + Xi.T @ Q
print("Tangent constraint matrix:", np.round(constraint, 10))
check_true(np.linalg.norm(constraint) < 1e-10, "Stiefel tangent condition holds")
Y = qr_retraction(Q + 0.2 * Xi)
check_true(np.linalg.norm(Y.T @ Y - np.eye(2)) < 1e-10, "QR retraction returns orthonormal columns")
print("Retracted columns are orthonormal.")

Demo 23: Tangent-space updates and retractions

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 49

header("Demo 23 - Tangent-space updates and retractions: SPD manifold distance")
A = spd_from_eigs([1.0, 3.0])
B = spd_from_eigs([2.0, 5.0])
d_ab = spd_distance(A, B)
d_ba = spd_distance(B, A)
print("Distance A to B:", round(d_ab, 6))
print("Distance B to A:", round(d_ba, 6))
check_close(d_ab, d_ba, tol=1e-10, name="SPD distance symmetry")
check_true(d_ab > 0, "distinct SPD matrices have positive distance")

Demo 24: Matrix manifolds in ML

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 51

header("Demo 24 - Matrix manifolds in ML: exponential vs retraction")
x = normalize(np.array([1.0, 0.0, 0.0]))
v = tangent_projection_sphere(x, np.array([0.0, 0.2, 0.1]))
exp_point = exp_sphere(x, v)
ret_point = retract_sphere(x, v)
print("Exponential point:", np.round(exp_point, 5).tolist())
print("Retraction point:", np.round(ret_point, 5).tolist())
check_true(abs(np.linalg.norm(exp_point) - 1.0) < 1e-10, "exponential stays on sphere")
check_true(abs(np.linalg.norm(ret_point) - 1.0) < 1e-10, "retraction stays on sphere")
print("Distance between exp and retraction:", round(float(np.linalg.norm(exp_point - ret_point)), 6))

Demo 25: Optimization pipeline: gradient step retract transport

This demo turns the approved TOC concept into a small executable geometric calculation.

Code cell 53

header("Demo 25 - Optimization pipeline: gradient step retract transport: sphere tangent projection")
x = normalize(np.array([1.0, 1.0, 1.0]))
v = np.array([2.0, -1.0, 0.5])
tangent = tangent_projection_sphere(x, v)
print("Point on sphere:", np.round(x, 3).tolist())
print("Projected tangent:", np.round(tangent, 3).tolist())
check_close(np.dot(x, tangent), 0.0, tol=1e-10, name="orthogonal to base point")