/algorithm-designer

name: algorithm-designer
description: Design and document statistical algorithms with pseudocode and complexity analysis

Algorithm Designer

You are an expert in designing and documenting statistical algorithms.

Algorithm Documentation Standards

Required Components

1. **Purpose**: What problem does this solve? 2. **Input/Output**: Precise specifications 3. **Pseudocode**: Language-agnostic description 4. **Complexity**: Time and space analysis 5. **Convergence**: Conditions and guarantees 6. **Implementation notes**: Practical considerations

Input/Output Specification

Formal Specification Template

Every algorithm must have precise input/output documentation:

INPUT SPECIFICATION:
- Data: D = {(Y_i, A_i, M_i, X_i)}_{i=1}^n where:
  - Y_i ∈ ℝ (continuous outcome)
  - A_i ∈ {0,1} (binary treatment)
  - M_i ∈ ℝ^d (d-dimensional mediator)
  - X_i ∈ ℝ^p (p covariates)
- Parameters: θ ∈ Θ ⊆ ℝ^k (parameter space)
- Tolerance: ε > 0 (convergence criterion)
- Max iterations: T_max ∈ ℕ

OUTPUT SPECIFICATION:
- Estimate: θ̂ ∈ ℝ^k (point estimate)
- Variance: V̂ ∈ ℝ^{k×k} (covariance matrix)
- Convergence: boolean (did algorithm converge?)
- Iterations: t ∈ ℕ (iterations used)

# R implementation of formal I/O specification
define_algorithm_io <- function() {
  list(
    input = list(
      data = "data.frame with columns Y, A, M, X",
      params = "list(tol = 1e-6, max_iter = 1000)",
      models = "list(outcome_formula, mediator_formula, propensity_formula)"
    ),
    output = list(
      estimate = "numeric vector of parameter estimates",
      se = "numeric vector of standard errors",
      vcov = "variance-covariance matrix",
      converged = "logical indicating convergence",
      iterations = "integer count of iterations"
    ),
    complexity = list(
      time = "O(n * p^2) per iteration",
      space = "O(n * p)",
      iterations = "O(log(1/epsilon)) for Newton-type"
    )
  )
}

---

Convergence Criteria

Standard Convergence Conditions

| Criterion | Formula | Use Case | |-----------|---------|----------| | Absolute | $\|\theta^{(t+1)} - \theta^{(t)}\| < \varepsilon$ | Parameter convergence | | Relative | $\|\theta^{(t+1)} - \theta^{(t)}\|/\|\theta^{(t)}\| < \varepsilon$ | Scale-invariant | | Gradient | $\|\nabla L(\theta^{(t)})\| < \varepsilon$ | Optimization | | Function | $\|L(\theta^{(t+1)}) - L(\theta^{(t)})\| < \varepsilon$ | Objective convergence | | Cauchy | $\max_{i} |\theta_i^{(t+1)} - \theta_i^{(t)}| < \varepsilon$ | Component-wise |

Mathematical Formulation

**Convergence tolerance**: $\varepsilon = 10^{-6}$ (typical default)

**Standard tolerances by application**:

• Numerical optimization: $\varepsilon = 10^{-8}$
• Statistical estimation: $\varepsilon = 10^{-6}$
• Approximate methods: $\varepsilon = 10^{-4}$

Complexity Formulas

**Linear complexity** $O(n)$: Operations grow proportionally to input size $$T(n) = c \cdot n + O(1)$$

**Quadratic complexity** $O(n^2)$: Nested iterations over input $$T(n) = c \cdot n^2 + O(n)$$

**Linearithmic complexity** $O(n \log n)$: Divide-and-conquer with linear work per level $$T(n) = c \cdot n \log_2 n + O(n)$$

**Space-Time Tradeoff**: $$\text{Time} \times \text{Space} \geq \Omega(\text{Information Content})$$

**Convergence rate analysis**:

• Linear convergence: $\|\theta^{(t)} - \theta^*\| \leq C \cdot \rho^t$ where $0 < \rho < 1$
• Quadratic convergence: $\|\theta^{(t+1)} - \theta^*\| \leq C \cdot \|\theta^{(t)} - \theta^*\|^2$
• Superlinear: $\lim_{t \to \infty} \frac{\|\theta^{(t+1)} - \theta^*\|}{\|\theta^{(t)} - \theta^*\|} = 0$

# Comprehensive convergence checking
check_convergence <- function(theta_new, theta_old, gradient = NULL,
                              objective_new = NULL, objective_old = NULL,
                              tol = 1e-6, method = "relative") {
  switch(method,
    "absolute" = {
      # |θ^(t+1) - θ^t| < ε
      converged <- max(abs(theta_new - theta_old)) < tol
      criterion <- max(abs(theta_new - theta_old))
    },
    "relative" = {
      # |θ^(t+1) - θ^t| / |θ^t| < ε
      denom <- pmax(abs(theta_old), 1)  # Avoid division by zero
      converged <- max(abs(theta_new - theta_old) / denom) < tol
      criterion <- max(abs(theta_new - theta_old) / denom)
    },
    "gradient" = {
      # |∇L(θ)| < ε
      stopifnot(!is.null(gradient))
      converged <- sqrt(sum(gradient^2)) < tol
      criterion <- sqrt(sum(gradient^2))
    },
    "objective" = {
      # |L(θ^(t+1)) - L(θ^t)| < ε
      stopifnot(!is.null(objective_new), !is.null(objective_old))
      converged <- abs(objective_new - objective_old) < tol
      criterion <- abs(objective_new - objective_old)
    }
  )

  list(converged = converged, criterion = criterion, method = method)
}

# Newton-Raphson with convergence monitoring
newton_raphson <- function(f, grad, hess, theta0, tol = 1e-6, max_iter = 100) {
  theta <- theta0
  history <- list()

  for (t in 1:max_iter) {
    g <- grad(theta)
    H <- hess(theta)

    # Newton step: θ^(t+1) = θ^t - H^(-1) * g
    # Time complexity: O(p^3) for matrix inversion
    delta <- solve(H, g)
    theta_new <- theta - delta

    # Check convergence
    conv <- check_convergence(theta_new, theta, gradient = g, tol = tol)
    history[[t]] <- list(theta = theta, gradient_norm = sqrt(sum(g^2)))

    if (conv$converged) {
      return(list(
        estimate = theta_new,
        iterations = t,
        converged = TRUE,
        history = history
      ))
    }

    theta <- theta_new
  }

  list(estimate = theta, iterations = max_iter, converged = FALSE, history = history)
}

Complexity and Convergence Relationship

| Algorithm | Convergence Rate | Iterations to $\varepsilon$ | |-----------|-----------------|----------------------------| | Gradient Descent | $O(1/t)$ | $O(1/\varepsilon)$ | | Accelerated GD | $O(1/t^2)$ | $O(1/\sqrt{\varepsilon})$ | | Newton-Raphson | Quadratic | $O(\log\log(1/\varepsilon))$ | | EM Algorithm | Li