# GEODESIC SEARCH on S^7 — Algorithm Specification ## Component: Radial Self-Finding Experiment — Geodesic Search Engine **Author**: Geometric Physicist (Information Geometry) **Status**: DESIGN_PHASE **Dependencies**: `EXPERIMENT_RADIAL_SELF_FIND.md` --- ## 1. Mathematical Preliminaries ### 1.1 The Manifold: S^7 (Fisher Sphere) The state space is the 7-sphere embedded in R^8: ``` S^7 = {x ∈ R^8 : ||x||² = 1, x_i ≥ 0} ``` In √p-coordinates, each point x corresponds to a probability distribution p ∈ Δ_7 (the 7-simplex) via: ``` x = (√p_1, √p_2, ..., √p_8), where Σ_i p_i = 1, p_i ≥ 0 ``` The **Fisher-Rao metric** on Δ_7: ``` g_ij = δ_ij / p_i + 1/p_8 (for i,j = 1,...,7) ``` In √p-coordinates, this becomes the **round metric** on S^7 (the metric induced from the Euclidean metric on R^8). This is the Chentsov theorem: the Fisher-Rao metric on the simplex is isometric to the round sphere. ### 1.2 Tangent Space At a point x ∈ S^7, the tangent space is: ``` T_x S^7 = {v ∈ R^8 : ⟨v, x⟩ = 0} ``` A **direction** d at x is a unit vector in T_x S^7: ``` d ∈ T_x S^7, ||d|| = 1 ``` The set of all directions at x forms a 6-sphere S^6 ⊂ T_x S^7. **Projection to tangent space** (for any vector v ∈ R^8): ``` Proj_{T_x}(v) = v - ⟨v, x⟩ · x ``` ### 1.3 Geodesics on S^7 Geodesics on the round sphere are **great circles**. The geodesic starting at x in direction d with arc-length parameter t is: ``` γ_d(t) = cos(t) · x + sin(t) · d (G1) ``` **Properties:** - γ_d(0) = x - γ_d'(0) = d - ||γ_d(t)|| = 1 for all t (stays on S^7) - Arc-length parameter: ||γ_d'(t)|| = 1 **Range restriction**: Since x_i ≥ 0 (probability simplex), geodesics may leave the valid region. The maximal valid step is: ``` t_max(d) = min_i { arctan(-x_i / d_i) : d_i < 0 } (G2) ``` If all d_i ≥ 0, then t_max = +∞ (stays in positive orthant). ### 1.4 Exponential Map The exponential map at x ∈ S^7 is the geodesic flow: ``` exp_x(v) = cos(||v||) · x + sin(||v||) · (v / ||v||) (EM) ``` for v ∈ T_x S^7, v ≠ 0. For v = 0, exp_x(0) = x. **Inverse (logarithm map)**: For y ∈ S^7 not antipodal to x: ``` log_x(y) = arccos(⟨x, y⟩) · (y - ⟨x, y⟩ x) / ||y - ⟨x, y⟩ x|| (LM) ``` ### 1.5 Parallel Transport To transport a tangent vector v ∈ T_x S^7 along the geodesic from x to y = γ_d(t), use the standard formula for parallel transport on the sphere: ``` P_{x→y}(v) = v - (⟨v, y⟩ / (1 + ⟨x, y⟩)) · (x + y) (PT) ``` **Properties:** - P_{x→y}(v) ∈ T_y S^7 (remains tangent) - Isometric: ||P_{x→y}(v)|| = ||v|| - Inner-product preserving: ⟨P(v), P(w)⟩ = ⟨v, w⟩ --- ## 2. Direction Representation ### 2.1 Defining a Direction A direction d at point S_0 ∈ S^7 is represented as: ``` struct Direction { d: R^8, // unit vector constraints: { ⟨d, S_0⟩ = 0, // tangent to S^7 ||d|| = 1, // unit length } } ``` ### 2.2 Creating a Direction from an Arbitrary Vector Given v ∈ R^8, construct a valid direction at S_0: ``` function make_direction(S_0: S^7, v: R^8) → Direction: v_tangent = v - ⟨v, S_0⟩ · S_0 // Project to T_{S_0} S^7 if ||v_tangent|| < ε: return null // v is parallel to S_0, no valid direction d = v_tangent / ||v_tangent|| // Normalize return Direction{d} ``` ### 2.3 Direction ↔ Spiral Correspondence Each direction d corresponds to a family of probability distributions along the geodesic. The spiral index function maps these back to N: ``` n(d, t) = spiral_index(γ_d(t)) = argmin_n ||f(n) - γ_d(t)||² ``` where f(n) = (√n · cos(nψ), √n · sin(nψ)) is the Φ-corkscrew. The **compression ratio** along direction d is: ``` C(d, t) = compression_ratio(n(d, t)) = original_size / |RLE(phinary(n(d, t)))| ``` --- ## 3. Core Algorithm: `geodesic_search` ### 3.1 Pseudocode ``` function geodesic_search( S_0: S^7, // Starting point on Fisher sphere N_directions: int, // Number of directions to explore T_steps: int, // Number of steps along each geodesic t_max: float = π/2, // Maximum geodesic distance step_pattern: str = "geometric" // "uniform", "geometric", "adaptive" ) → (d*, t*, S*, C*): // ---- Step 1: Sample directions ---- D = sample_directions(S_0, N_directions) // ---- Step 2: Define step sizes ---- if step_pattern == "uniform": T = linspace(0, t_max, T_steps) elif step_pattern == "geometric": T = t_max * (1 - exp(-3 * linspace(0, 1, T_steps))) / (1 - e^{-3}) elif step_pattern == "adaptive": T = golden_section_steps(0, t_max, T_steps) // ---- Step 3: Evaluate all (direction, step) pairs ---- best_C = -∞ best_triple = null for i = 1 to N_directions: d_i = D[i] // Compute t_max for this direction (boundary of simplex) t_boundary = compute_boundary_limit(S_0, d_i) T_eff = T ∩ [0, t_boundary] // Clip to valid range for j = 1 to |T_eff|: t_j = T_eff[j] S_ij = γ_{d_i}(t_j) = cos(t_j)·S_0 + sin(t_j)·d_i // (G1) // Map to probability simplex p_ij = S_ij ⊙ S_ij // element-wise square: p_k = (S_ij,k)² p_ij = p_ij / Σ_k p_ij,k // renormalize (numerical safety) // Encode and measure compression n_ij = spiral_index(S_ij) C_ij = compression_ratio(n_ij) // Update best if C_ij > best_C: best_C = C_ij best_triple = (d_i, t_j, S_ij, C_ij) // ---- Step 4: Return best ---- (d*, t*, S*, C*) = best_triple return (d*, t*, S*, C*) ``` ### 3.2 Step Size Patterns **Uniform**: `t_j = j · t_max / T_steps` Good for: Initial broad exploration. **Geometric** (default): `t_j = t_max · (1 - e^{-3j/T_steps}) / (1 - e^{-3})` Good for: More samples near S_0 where changes are smaller. **Adaptive (Golden Section)**: Uses golden-ratio spacing to concentrate samples where compression changes most rapidly. ### 3.3 Boundary Computation ``` function compute_boundary_limit(S_0: S^7, d: Direction) → float: """ Compute the maximum geodesic distance before leaving the positive orthant (p_i ≥ 0 for all i). """ t_max_valid = +∞ for i = 1 to 8: if d_i < 0 and S_0,i > 0: // Solve: cos(t)·S_0,i + sin(t)·d_i = 0 // => tan(t) = -S_0,i / d_i t_i = arctan(-S_0,i / d_i) if t_i > 0: t_max_valid = min(t_max_valid, t_i) return t_max_valid ``` --- ## 4. Gradient Computation: `compute_gradient` ### 4.1 Problem Statement We need the Riemannian gradient of the compression ratio C with respect to position on S^7. Since C is evaluated through the discrete `spiral_index` function (argmin), it is piecewise constant with jumps at boundaries. We use **finite differences on the sphere**. ### 4.2 Finite Difference Gradient ``` function compute_gradient( S: S^7, // Current point f: S^7 → R, // Objective function (here: f(x) = C(spiral_index(x))) ε: float = 1e-4, // Finite difference step method: str = "central" // "forward", "central", "complex" ) → ∇f ∈ T_S S^7: // Build orthonormal basis for T_S S^7 B = build_tangent_basis(S) // Returns {e_1, ..., e_7} ⊂ R^8 // Compute directional derivatives ∇f = 0 ∈ R^8 for k = 1 to 7: e_k = B[k] if method == "central": S_plus = exp_S(ε · e_k) = cos(ε)·S + sin(ε)·e_k S_minus = exp_S(-ε · e_k) = cos(ε)·S - sin(ε)·e_k ∂f/∂e_k = (f(S_plus) - f(S_minus)) / (2ε) elif method == "forward": S_plus = exp_S(ε · e_k) ∂f/∂e_k = (f(S_plus) - f(S)) / ε ∇f += (∂f/∂e_k) · e_k return ∇f // Lives in T_S S^7 by construction ``` ### 4.3 Tangent Basis Construction (Gram-Schmidt) ``` function build_tangent_basis(S: S^7) → List[R^8]: """ Build orthonormal basis for T_S S^7 using modified Gram-Schmidt. """ basis = [] for i = 1 to 8: e = standard_basis_vector(i) // (0,...,1,...,0) // Project to tangent space e = e - ⟨e, S⟩ · S // Orthogonalize against existing basis vectors for b in basis: e = e - ⟨e, b⟩ · b // Normalize (skip if numerically zero) if ||e|| > ε_machine: e = e / ||e|| basis.append(e) // Verify: should have exactly 7 basis vectors assert len(basis) == 7 return basis ``` ### 4.4 Gradient via Automatic Differentiation (Alternative) If the spiral_index function is differentiable (e.g., using a smooth relaxation), we can use automatic differentiation: ``` // Smooth relaxation of spiral_index: // Instead of argmin, use soft-argmin: n_smooth(x) = Σ_n n · softmax(-β·||f(n) - x||²)_n // Where β > 0 is an inverse temperature parameter. // As β → ∞, n_smooth → argmin. // Typical values: β ∈ [10, 1000] // Chain rule: ∇_x C = (∂C/∂n) · (∂n_smooth/∂x) // Project to tangent space: ∇_S C = Proj_{T_S}(∇_x C) ``` ### 4.5 Gradient with Respect to Direction When optimizing over directions d ∈ S^6 (the direction sphere at S_0): ``` function compute_direction_gradient( S_0: S^7, // Fixed starting point d: Direction, // Current direction f: S^7 → R, // f(x) = C(spiral_index(x)) ε: float = 1e-4 ) → ∇_d F ∈ T_d S^6: // Define F(d) = max_t f(γ_d(t)) [or f(γ_d(t*)) for fixed t*] // Build basis for T_d S^6 = {v ∈ T_{S_0} S^7 : ⟨v, d⟩ = 0} B = build_direction_basis(S_0, d) // Returns 6 orthonormal vectors // Finite differences on the direction sphere ∇F = 0 ∈ R^8 for k = 1 to 6: b_k = B[k] // Exponential map on S^6 (embedded in T_{S_0} S^7): d_plus = cos(ε)·d + sin(ε)·b_k d_minus = cos(ε)·d - sin(ε)·b_k F_plus = evaluate_along_geodesic(S_0, d_plus, f) F_minus = evaluate_along_geodesic(S_0, d_minus, f) ∂F/∂b_k = (F_plus - F_minus) / (2ε) ∇F += (∂F/∂b_k) · b_k return ∇F function build_direction_basis(S_0: S^7, d: Direction) → List[R^8]: """ Build orthonormal basis for T_d S^6, which is the subspace of T_{S_0} S^7 orthogonal to d. """ // Start with full tangent basis at S_0 B_full = build_tangent_basis(S_0) // 7 vectors // Remove component along d basis = [] for b in B_full: b_orth = b - ⟨b, d⟩ · d if ||b_orth|| > ε_machine: b_orth = b_orth / ||b_orth|| basis.append(b_orth) // Should now have 6 vectors assert len(basis) == 6 return basis ``` --- ## 5. Core Algorithm: `gradient_ascent_step` ### 5.1 Standard Riemannian Gradient Ascent ``` function gradient_ascent_step( S_current: S^7, // Current point f: S^7 → R, // Objective: compression ratio learning_rate: float, // Step size η grad_method: str = "fd" // "fd" or "smooth" ) → S_next: // Compute Riemannian gradient if grad_method == "fd": ∇f = compute_gradient(S_current, f, ε=1e-4, method="central") elif grad_method == "smooth": ∇f = compute_smooth_gradient(S_current, f, β=100) // Check for stationarity if ||∇f|| < ε_stationary: return S_current // Local maximum reached // Normalize for unit-step control (optional) // ∇f = ∇f / ||∇f|| // Riemannian gradient ascent via exponential map: v = learning_rate · ∇f S_next = exp_{S_current}(v) = cos(||v||) · S_current + sin(||v||) · (v / ||v||) // (EM) // Ensure we stay in the positive orthant if any(S_next,i < 0): // Project to nearest valid point on S^7 ∩ {x_i ≥ 0} S_next = project_to_simplex_sphere(S_next) return S_next ``` ### 5.2 Adaptive Step Size (Line Search) ``` function gradient_ascent_step_adaptive( S_current: S^7, f: S^7 → R, η_init: float = 0.1, α: float = 0.5, // Backtracking factor c: float = 1e-4, // Armijo constant max_backtrack: int = 10 ) → S_next: ∇f = compute_gradient(S_current, f) grad_norm = ||∇f|| if grad_norm < ε_stationary: return S_current // Unit-norm gradient direction g = ∇f / grad_norm η = η_init f_current = f(S_current) for step = 1 to max_backtrack: S_trial = cos(η·grad_norm)·S_current + sin(η·grad_norm)·g S_trial = project_to_simplex_sphere(S_trial) f_trial = f(S_trial) // Armijo condition: sufficient increase if f_trial ≥ f_current + c · η · grad_norm²: return S_trial η = α · η // Backtrack // No improvement found — return current point return S_current ``` ### 5.3 Projection to Valid Simplex Sphere ``` function project_to_simplex_sphere(y: R^8) → S^7: """ Project y to the nearest point on S^7 ∩ {x_i ≥ 0}. Uses alternating projection: sphere → orthant → sphere → ... """ x = y / ||y|| // Project to sphere for iter = 1 to max_proj_iter: // Clip negative components x_clipped = max(x, 0) // element-wise // If already in orthant, done if ||x_clipped - x|| < ε: return x_clipped / ||x_clipped|| // Project back to sphere x = x_clipped / ||x_clipped|| return x ``` ### 5.4 Momentum-Based Gradient Ascent on Manifold For faster convergence, use Riemannian momentum: ``` function gradient_ascent_momentum( S_current: S^7, v_prev: T_{S_prev} S^7, // Previous velocity (in old tangent space) S_prev: S^7, // Previous point f: S^7 → R, η: float, // Learning rate μ: float = 0.9 // Momentum coefficient ) → (S_next, v_next): // Compute gradient at current point ∇f = compute_gradient(S_current, f) // Transport previous velocity to current tangent space if v_prev is not null: v_transport = P_{S_prev → S_current}(v_prev) // (PT) else: v_transport = 0 // Update velocity: v = μ · v_transport + ∇f v = μ · v_transport + ∇f // in T_{S_current} S^7 // Riemannian gradient step with momentum S_next = exp_{S_current}(η · v) = cos(η·||v||)·S_current + sin(η·||v||)·(v/||v||) S_next = project_to_simplex_sphere(S_next) return (S_next, v) ``` --- ## 6. Core Algorithm: `sample_directions` ### 6.1 Uniform Sampling on T_{S_0} S^7 ``` function sample_directions( S_0: S^7, N: int, distribution: str = "uniform" // "uniform", "golden", "gradient_biased" ) → List[Direction]: if distribution == "uniform": return sample_uniform_directions(S_0, N) elif distribution == "golden": return sample_golden_directions(S_0, N) elif distribution == "gradient_biased": ∇f = compute_gradient(S_0, f) return sample_gradient_biased(S_0, N, ∇f) ``` ### 6.2 Uniform Random Directions ``` function sample_uniform_directions(S_0: S^7, N: int) → List[Direction]: """ Sample N i.i.d. uniformly distributed directions on the unit sphere S^6 ⊂ T_{S_0} S^7. Algorithm: Sample standard Gaussian in R^8, project to tangent space, normalize. This gives uniform distribution on S^6. """ directions = [] for _ = 1 to N: v = randn(8) // N(0, I_8) v_tangent = v - ⟨v, S_0⟩ · S_0 d = v_tangent / ||v_tangent|| directions.append(Direction{d}) return directions ``` ### 6.3 Golden-Angle Directions (Φ-Quasi-Monte Carlo) ``` function sample_golden_directions(S_0: S^7, N: int) → List[Direction]: """ Generate N directions using the generalized golden angle on S^6. Provides low-discrepancy coverage of the direction sphere. Uses the Fibonacci lattice generalized to 7 dimensions with the golden ratio φ = (1+√5)/2. """ φ = (1 + √5) / 2 // Build tangent basis B = build_tangent_basis(S_0) // {e_1, ..., e_7} directions = [] for n = 1 to N: // Generalized Fibonacci angles in 7D // Use the first 7 powers of φ modulo 1 angles = [] for k = 1 to 7: α_k = 2π · fractional_part(n · φ^{-k}) angles.append(α_k) // Convert angles to direction in the 7D tangent space // Use hyperspherical coordinates d_local = angles_to_unit_vector_7d(angles) // Map from local 7D coordinates to R^8 via basis d = Σ_{k=1}^7 d_local,k · e_k directions.append(Direction{d}) return directions ``` **Hyperspherical coordinate mapping (7D → unit vector)**: ``` function angles_to_unit_vector_7d(angles: [θ_1,...,θ_7]) → R^7: """ Convert 7 angles to unit vector in R^7. x_1 = cos(θ_1) x_2 = sin(θ_1)·cos(θ_2) x_3 = sin(θ_1)·sin(θ_2)·cos(θ_3) ... x_7 = sin(θ_1)·...·sin(θ_6)·cos(θ_7) x_8 = sin(θ_1)·...·sin(θ_6)·sin(θ_7) (We use 7 angles for S^6, embedded in R^7 ~ T_{S_0} S^7) """ x = zeros(7) sin_prod = 1.0 for k = 1 to 6: x[k] = sin_prod · cos(angles[k]) sin_prod = sin_prod · sin(angles[k]) x[7] = sin_prod · cos(angles[7]) // Note: for S^6 in R^7, we need 6 angles, not 7 // The above gives 7D; for S^6 we use 6 angles ``` *Correction*: S^6 requires 6 angles in the hyperspherical coordinate system. The generalized Fibonacci sequence on S^d uses d angles. ``` function sample_golden_directions_v2(S_0: S^7, N: int) → List[Direction]: φ = (1 + √5) / 2 B = build_tangent_basis(S_0) directions = [] for n = 1 to N: // 6D generalized Fibonacci point t = n / N // Use 6 angular coordinates coords = zeros(7) for k = 1 to 7: // Coordinate k = fractional part of n * φ^k, mapped to [-1, 1] coords[k] = 2 · fractional_part(n · φ^{(k-1)/7}) - 1 // Project to unit sphere in the tangent space d_local = coords / ||coords|| d = Σ_k d_local,k · B[k] directions.append(Direction{d}) return directions ``` ### 6.4 Gradient-Biased Direction Sampling ``` function sample_gradient_biased( S_0: S^7, N: int, ∇f: T_{S_0} S^7, concentration: float = 5.0 ) → List[Direction]: """ Sample directions biased toward the gradient direction. Uses a von Mises-Fisher distribution on S^6 with mean direction ∇f/||∇f|| and concentration parameter κ. """ μ = ∇f / ||∇f|| // Mean direction κ = concentration // Concentration directions = [] for _ = 1 to N: // Rejection sampling for von Mises-Fisher // or: sample from N(κ·μ, I), project to tangent space, normalize v = randn(8) + κ · μ v_tangent = v - ⟨v, S_0⟩ · S_0 d = v_tangent / ||v_tangent|| directions.append(Direction{d}) return directions ``` --- ## 7. Complete Self-Finding Loop Integration ### 7.1 Main Loop (from experiment spec) ``` function self_finding_loop( S_0: S^7, // Initial state max_iter: int = 100, // Maximum iterations N_directions: int = 32, // Directions per search T_steps: int = 16, // Steps per geodesic η: float = 0.05, // Learning rate ε_converge: float = 1e-6 // Convergence threshold ) → (S_opt, history): S = S_0 history = [] for iter = 1 to max_iter: // ---- Phase 1: Direction Search ---- (d*, t*, S_candidate, C*) = geodesic_search( S, N_directions, T_steps ) C_current = compression_ratio(spiral_index(S)) // ---- Phase 2: Gradient Ascent Step ---- if C* > C_current: // Move toward candidate f(x) = compression_ratio(spiral_index(x)) S = gradient_ascent_step(S, f, η, grad_method="fd") else: // Local maximum reached break // ---- Phase 3: Self-Reference (Strange Loop) ---- // Encode the search trajectory itself trajectory = history + [(d*, t*, C*)] n_exp = encode_trajectory(trajectory) // ---- Phase 4: Meta-Update ---- // Optionally: use the experiment encoding as next starting point // S_meta = spiral_to_sphere_point(n_exp) // This closes the strange loop history.append({ 'iteration': iter, 'S': S, 'C': C_current, 'd_best': d*, 't_best': t*, 'n_exp': n_exp }) // Convergence check if len(history) > 1: C_prev = history[-2]['C'] if |C_current - C_prev| / C_prev < ε_converge: break return (S, history) ``` ### 7.2 Encoding the Trajectory ``` function encode_trajectory(trajectory: List[(d, t, C)]) → int: """ Encode the search trajectory as a spiral index. This creates the 'strange loop' — the search encodes itself. The trajectory is a sequence of (direction, step, compression) tuples. We serialize this and encode as a DNA sequence via the Φ-corkscrew. """ // Serialize trajectory to a bit string bits = serialize(trajectory) // e.g., using IEEE 754 floats // Convert bit string to integer n_trajectory = bit_string_to_int(bits) // Map to spiral index n_exp = spiral_index_from_int(n_trajectory) return n_exp ``` --- ## 8. Convergence Analysis ### 8.1 Convergence of Geodesic Search **Theorem** (Local convergence). Let C: S^7 → R be the compression ratio function. If C is L-smooth on S^7 (Riemannian gradient is L-Lipschitz) and bounded above, then gradient ascent with step size η ∈ (0, 2/L) converges to a stationary point. **Proof sketch**: 1. The exponential map exp_x is the standard geodesic flow on S^7 2. The sectional curvature of S^7 is constant K = 1 > 0 3. For positively curved manifolds, gradient ascent with appropriate step size decreases the objective gradient norm to zero 4. By compactness of S^7 and continuity (of the smooth relaxation), the sequence has a convergent subsequence 5. The limit point satisfies ∇C = 0 (stationary point) ### 8.2 Rate of Convergence For gradient ascent on S^7 with constant step size η: ``` C(S_{k+1}) - C(S*) ≤ (1 - η·m) · (C(S_k) - C(S*)) [if m-strongly convex] ``` where m is the strong convexity parameter (lower bound on Hessian eigenvalues in the tangent space). For non-convex C (the general case): ``` min_{0≤k≤K} ||∇C(S_k)||² ≤ (C(S_0) - C(S*)) / (η·K) ``` O(1/√K) rate for gradient norm, O(1/K) for strongly convex regions. ### 8.3 Effect of Self-Reference When the trajectory encoding is used as the next starting point: ``` S_{k+1} = f(S_k, trajectory(S_0, ..., S_k)) ``` This creates a **non-Markovian** dynamical system. The convergence properties depend on the encoding map. Empirically, self-reference acts as a form of momentum: the encoded trajectory contains information about the "shape" of the compression landscape, which guides future searches. --- ## 9. Computational Complexity | Operation | Complexity | Notes | |-----------|-----------|-------| | `sample_directions` | O(N · d) | d=8, N directions | | `γ_d(t)` (geodesic) | O(d) | One evaluation | | `geodesic_search` | O(N · T · (d + T_spiral)) | T_spiral = spiral_index cost | | `compute_gradient` | O(2d · T_spiral) | 2·7 = 14 function evals | | `gradient_ascent_step` | O(2d · T_spiral) | Same as gradient | | `parallel_transport` | O(d) | Negligible | | Full iteration | O(N·T·T_spiral) | Dominated by spiral_index | Where `T_spiral` = cost of `spiral_index(x)` = O(N_spiral) for a search over N_spiral spiral points (can be accelerated with KD-trees to O(log N_spiral)). --- ## 10. Summary of Key Formulas | Formula | Equation | Description | |---------|----------|-------------| | Geodesic | γ_d(t) = cos(t)·x + sin(t)·d | (G1) Great circle on S^7 | | Exponential Map | exp_x(v) = cos(||v||)·x + sin(||v||)·(v/||v||) | (EM) Geodesic flow | | Logarithm Map | log_x(y) = arccos(⟨x,y⟩)·(y-⟨x,y⟩x)/||...|| | (LM) Inverse exponential | | Parallel Transport | P_{x→y}(v) = v - (⟨v,y⟩/(1+⟨x,y⟩))·(x+y) | (PT) Along geodesic | | Tangent Projection | Proj_{T_x}(v) = v - ⟨v,x⟩·x | Project to tangent space | | Boundary Limit | t_max = min_i arctan(-x_i/d_i) for d_i < 0 | (G2) Simplex boundary | | Gradient (FD) | ∇f = Σ_k [f(exp(ε·e_k)) - f(exp(-ε·e_k))]/(2ε) · e_k | Central differences | | Gradient Step | S_{next} = exp_S(η·∇f) | Riemannian ascent | --- ## 11. Interface Specification ### 11.1 Function Signatures ```python def geodesic_search( S_0: np.ndarray, # shape (8,), unit vector, S_0 >= 0 N_directions: int, # number of directions to sample T_steps: int, # number of steps per geodesic t_max: float = np.pi/2, # maximum geodesic distance step_pattern: str = "geometric", f: Callable = compression_ratio_at_point # f: S^7 -> R ) -> Tuple[np.ndarray, float, np.ndarray, float]: """Returns (d_best, t_best, S_best, C_best)""" def gradient_ascent_step( S_current: np.ndarray, # shape (8,), current point f: Callable, # objective function f: S^7 -> R learning_rate: float, # step size η grad_method: str = "fd", # "fd" or "smooth" ε: float = 1e-4 # finite difference step ) -> np.ndarray: """Returns S_next (shape (8,))""" def sample_directions( S_0: np.ndarray, # shape (8,), anchor point N: int, # number of directions distribution: str = "uniform" # "uniform", "golden", "gradient_biased" ) -> np.ndarray: """Returns directions (shape (N, 8)), each row is a unit tangent vector""" def compute_gradient( S: np.ndarray, # shape (8,), evaluation point f: Callable, # f: S^7 -> R ε: float = 1e-4, # FD step size method: str = "central" # "central", "forward" ) -> np.ndarray: """Returns gradient (shape (8,)), in T_S S^7""" def parallel_transport( v: np.ndarray, # shape (8,), tangent vector at x x: np.ndarray, # shape (8,), source point y: np.ndarray # shape (8,), target point ) -> np.ndarray: """Returns transported vector (shape (8,)), in T_y S^7""" ``` ### 11.2 Type Definitions ```python S7Point = np.ndarray # shape (8,), ||x|| = 1, x_i >= 0 Direction = np.ndarray # shape (8,), ||d|| = 1, = 0 TangentV = np.ndarray # shape (8,), = 0 ScalarFn = Callable[[S7Point], float] ``` --- ## 12. Testing Strategy ### 12.1 Unit Tests 1. **Geodesic stays on sphere**: For random x, d, t: verify ||γ_d(t)|| = 1 2. **Tangent constraint**: Verify ⟨d, S_0⟩ = 0 for all sampled directions 3. **Unit norm**: Verify ||d|| = 1 for all directions 4. **Gradient is tangent**: Verify ⟨∇f, S⟩ = 0 5. **Gradient ascent increases f**: For a known function, verify f(S_next) > f(S) 6. **Parallel transport is isometric**: Verify ||P(v)|| = ||v|| 7. **Exponential and log are inverses**: Verify log_x(exp_x(v)) ≈ v ### 12.2 Integration Tests 1. **Known optimum**: Test on a function with known maximum on S^7 2. **Symmetry**: Verify gradient is zero at symmetric points for symmetric functions 3. **Convergence**: Verify ||∇f|| → 0 as iterations increase 4. **Boundary handling**: Verify geodesics don't leave the simplex ### 12.3 End-to-End Test ``` S_0 = uniform point on S^7 ∩ {x_i >= 0} (d*, t*, S*, C*) = geodesic_search(S_0, N=64, T=32) assert C* >= compression_ratio(spiral_index(S_0)) assert ||S*|| ≈ 1 assert all(S* >= 0) ``` --- ## 13. References 1. Amari, S.-I. (2016). *Information Geometry and Its Applications*. Springer. 2. Chentsov, N.N. (1982). *Statistical Decision Rules and Optimal Inference*. 3. Absil, P.-A., Mahony, R., & Sepulchre, R. (2008). *Optimization Algorithms on Matrix Manifolds*. Princeton University Press. 4. Boumal, N. (2023). *An Introduction to Optimization on Smooth Manifolds*. Cambridge University Press. 5. Audet, C. & Dennis, J.E. (2006). "Mesh Adaptive Direct Search Algorithms for Constrained Optimization*. SIAM Journal on Optimization. --- *End of Geodesic Search Algorithm Specification*