Research-Stack/docs/experiment_geodesic_search.md

# 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, <d, x> = 0
TangentV  = np.ndarray  # shape (8,), <v, x> = 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*