Imported from gabrielmoreira/agent-skills-mirror (mirrors/repos/brycewang-stanford@Auto-Empirical-Research-Skills/skills/26-Data-Wise-scholar/skills/implementation/algorithm-designer/SKILL.md). Install upstream with npx skills add gabrielmoreira/agent-skills-mirror --skill algorithm-designer. Copyright stays with the author.

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

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

Pseudocode Conventions

Standard Format

ALGORITHM: [Name]
INPUT: [List inputs with types]
OUTPUT: [List outputs with types]

1. [Initialize]
2. [Main loop or procedure]
   2.1 [Sub-step]
   2.2 [Sub-step]
3. [Return]

Example: AIPW Estimator

ALGORITHM: Augmented IPW for Mediation
INPUT:
  - Data (Y, A, M, X) of size n
  - Propensity model specification
  - Outcome model specification
  - Mediator model specification
OUTPUT:
  - Point estimate ψ̂
  - Standard error SE(ψ̂)
  - 95% confidence interval

1. ESTIMATE NUISANCE FUNCTIONS
   1.1 Fit propensity score: π̂(x) = P̂(A=1|X=x)
   1.2 Fit mediator density: f̂(m|a,x)
   1.3 Fit outcome regression: μ̂(a,m,x) = Ê[Y|A=a,M=m,X=x]

2. COMPUTE PSEUDO-OUTCOMES
   For i = 1 to n:
     2.1 Compute IPW weight: w_i = A_i/π̂(X_i) + (1-A_i)/(1-π̂(X_i))
     2.2 Compute outcome prediction: μ̂_i = μ̂(A_i, M_i, X_i)
     2.3 Compute augmentation term
     2.4 φ_i = w_i(Y_i - μ̂_i) + [integration term]

3. ESTIMATE AND INFERENCE
   3.1 ψ̂ = n⁻¹ Σᵢ φ_i
   3.2 SE = √(n⁻¹ Σᵢ (φ_i - ψ̂)²)
   3.3 CI = [ψ̂ - 1.96·SE, ψ̂ + 1.96·SE]

4. RETURN (ψ̂, SE, CI)

Complexity Analysis

Big-O Notation Guide

Formal Definition: $f(n) = O(g(n))$ if $\exists c, n_0$ such that $f(n) \leq c \cdot g(n)$ for all $n \geq n_0$

Key Formulas

Master Theorem for recurrences $T(n) = aT(n/b) + f(n)$:

• If $f(n) = O(n^{\log_b a - \epsilon})$ then $T(n) = \Theta(n^{\log_b a})$
• If $f(n) = \Theta(n^{\log_b a})$ then $T(n) = \Theta(n^{\log_b a} \log n)$
• If $f(n) = \Omega(n^{\log_b a + \epsilon})$ then $T(n) = \Theta(f(n))$

Sorting lower bound: Any comparison-based sort requires $\Omega(n \log n)$ comparisons

Matrix operations:

• Naive multiplication: $O(n^3)$
• Strassen: $O(n^{2.807})$
• Matrix inversion: $O(n^3)$ (same as multiplication)

# Complexity analysis helper
analyze_complexity <- function(f, n_values = c(100, 500, 1000, 5000)) {
  times <- sapply(n_values, function(n) {
    system.time(f(n))[["elapsed"]]
  })

  # Fit log-log regression to estimate complexity
  fit <- lm(log(times) ~ log(n_values))
  estimated_power <- coef(fit)[2]

  list(
    times = data.frame(n = n_values, time = times),
    estimated_complexity = paste0("O(n^", round(estimated_power, 2), ")"),
    power = estimated_power
  )
}

Statistical Algorithm Complexities

Template for Analysis

TIME COMPLEXITY:
- Initialization: O(...)
- Per iteration: O(...)
- Total (T iterations): O(...)
- Convergence typically in T = O(...) iterations

SPACE COMPLEXITY:
- Data storage: O(n × p)
- Working memory: O(...)
- Output: O(...)

Convergence Analysis

Types of Convergence

1. Finite termination: Exact solution in finite steps
2. Linear: $|x_{k+1} - x^| \leq c|x_k - x^|$, $c < 1$
3. Superlinear: $|x_{k+1} - x^| / |x_k - x^| \to 0$
4. Quadratic: $|x_{k+1} - x^| \leq c|x_k - x^|^2$

Convergence Documentation Template

CONVERGENCE:
- Type: [Linear/Superlinear/Quadratic]
- Rate: [Expression]
- Conditions: [What must hold]
- Stopping criterion: [When to stop]
- Typical iterations: [Order of magnitude]

Optimization Algorithms

Gradient-Based Methods

ALGORITHM: Gradient Descent
INPUT: f (objective), ∇f (gradient), x₀ (initial), η (step size), ε (tolerance)
OUTPUT: x* (minimizer)

1. k ← 0
2. WHILE ‖∇f(xₖ)‖ > ε:
   2.1 xₖ₊₁ ← xₖ - η∇f(xₖ)
   2.2 k ← k + 1
3. RETURN xₖ

COMPLEXITY: O(iterations × gradient_cost)
CONVERGENCE: Linear with rate (1 - η·μ) for μ-strongly convex f

Newton's Method

ALGORITHM: Newton-Raphson
INPUT: f, ∇f, ∇²f, x₀, ε
OUTPUT: x*

1. k ← 0
2. WHILE ‖∇f(xₖ)‖ > ε:
   2.1 Solve ∇²f(xₖ)·d = -∇f(xₖ) for direction d
   2.2 xₖ₊₁ ← xₖ + d
   2.3 k ← k + 1
3. RETURN xₖ

COMPLEXITY: O(iterations × p³) for p-dimensional
CONVERGENCE: Quadratic near solution

EM Algorithm Template

ALGORITHM: Expectation-Maximization
INPUT: Data Y, model parameters θ₀, tolerance ε
OUTPUT: MLE θ̂

1. θ ← θ₀
2. REPEAT:
   2.1 E-STEP: Compute Q(θ'|θ) = E[log L(θ'|Y,Z) | Y, θ]
   2.2 M-STEP: θ_new ← argmax_θ' Q(θ'|θ)
   2.3 Δ ← |θ_new - θ|
   2.4 θ ← θ_new
3. UNTIL Δ < ε
4. RETURN θ

CONVERGENCE: Monotonic increase in likelihood
             Linear rate near optimum

Bootstrap Algorithms

Nonparametric Bootstrap

ALGORITHM: Nonparametric Bootstrap
INPUT: Data X of size n, statistic T, B (number of replicates)
OUTPUT: SE estimate, CI

1. FOR b = 1 to B:
   1.1 Draw X*_b by sampling n observations with replacement from X
   1.2 Compute T*_b = T(X*_b)
2. SE_boot ← SD({T*_1, ..., T*_B})
3. CI_percentile ← [quantile(T*, 0.025), quantile(T*, 0.975)]
4. RETURN (SE_boot, CI_percentile)

COMPLEXITY: O(B × cost(T))
NOTES: B ≥ 1000 for SE, B ≥ 10000 for percentile CI

Parametric Bootstrap

ALGORITHM: Parametric Bootstrap
INPUT: Data X, parametric model M, B replicates
OUTPUT: SE estimate

1. Fit θ̂ = MLE(X, M)
2. FOR b = 1 to B:
   2.1 Generate X*_b ~ M(θ̂)
   2.2 Compute θ̂*_b = MLE(X*_b, M)
3. SE_boot ← SD({θ̂*_1, ..., θ̂*_B})
4. RETURN SE_boot

Numerical Stability Notes

Common Issues

1. Overflow/Underflow: Work on log scale
2. Cancellation: Reformulate subtractions
3. Ill-conditioning: Use regularization or pivoting
4. Convergence: Add damping or line search

Stability Techniques

# Log-sum-exp trick
log_sum_exp <- function(x) {
  max_x <- max(x)
  max_x + log(sum(exp(x - max_x)))
}

# Numerically stable variance
stable_var <- function(x) {
  n <- length(x)
  m <- mean(x)
  sum((x - m)^2) / (n - 1)  # One-pass with correction
}

Implementation Checklist

Before Coding

• Pseudocode written and reviewed
• Complexity analyzed
• Convergence conditions identified
• Edge cases documented
• Numerical stability considered

During Implementation

• Match pseudocode structure
• Add convergence monitoring
• Handle edge cases
• Log intermediate values (debug mode)
• Add early stopping

After Implementation

• Unit tests for components
• Integration tests for full algorithm
• Benchmark against reference implementation
• Profile for bottlenecks
• Document deviations from pseudocode

Key References

• CLRS
• Numerical Recipes