Random-key optimization can improve the efficiency of mixed-integer programs used in financial portfolio optimization.

Random-key optimization (RKO), when applied to mixed-integer programming (MIP) formulations of financial portfolio optimization problems, reduces computational solve time by ≥20% and/or improves solution quality (Sharpe ratio or objective value) by ≥5% compared to standard branch-and-bound MIP solvers (e.g., Gurobi, CPLEX) on benchmark portfolio instances with N ≥ 50 assets, cardinality constraints, and transaction cost terms, within a fixed computational budget of 3600 seconds per instance.

Minimum Viable Test (MVT): Compare RKO vs. Gurobi on 3 portfolio sizes (N=50, 200, 500) × 3 constraint types (cardinality-only, cardinality+transaction costs, cardinality+minimum lots) = 9 instance classes, 10 random instances each = 90 total instances. Measure: (a) solve time to 1% optimality gap, (b) best objective value at 3600s, (c) feasibility rate. Full Validation: Extend to N ∈ {50,100,200,500,1000}, 5 constraint types, 30 instances each = 750 instances. Include out-of-sample backtesting on real financial data (2010–2023 S&P 500 constituents). Statistical validation via bootstrap confidence intervals and Wilcoxon signed-rank tests.

# RKO Portfolio Optimization - Core Architecture

import numpy as np
from scipy.optimize import minimize
from dataclasses import dataclass
from typing import Tuple, List

@dataclass
class PortfolioInstance:
    mu: np.ndarray          # Expected returns [N]
    Sigma: np.ndarray       # Covariance matrix [N x N]
    k: int                  # Cardinality constraint
    w_min: float = 0.01     # Minimum weight per selected asset
    w_max: float = 0.30     # Maximum weight per selected asset
    transaction_costs: np.ndarray = None  # [N]
    w0: np.ndarray = None   # Current portfolio weights

class RKODecoder:
    """Maps continuous random keys in [0,1]^N to feasible portfolio"""
    
    def decode(self, keys: np.ndarray, instance: PortfolioInstance) -> np.ndarray:
        N = len(keys)
        # Step 1: Rank-based asset selection
        ranked_assets = np.argsort(keys)[::-1]  # Descending order
        selected = ranked_assets[:instance.k]
        
        # Step 2: Solve continuous QP for weights on selected assets
        weights = self._solve_subproblem(selected, instance)
        
        # Step 3: Map back to full weight vector
        w_full = np.zeros(N)
        w_full[selected] = weights
        return w_full
    
    def _solve_subproblem(self, selected: np.ndarray, 
                          instance: PortfolioInstance) -> np.ndarray:
        k = len(selected)
        Sigma_sub = instance.Sigma[np.ix_(selected, selected)]
        mu_sub = instance.mu[selected]
        lambda_risk = 0.5  # Risk aversion parameter
        
        def objective(w):
            portfolio_var = w @ Sigma_sub @ w
            portfolio_ret = mu_sub @ w
            tc = 0.0
            if instance.transaction_costs is not None:
                w0_sub = instance.w0[selected] if instance.w0 is not None \
                         else np.zeros(k)
                tc = instance.transaction_costs[selected] @ np.abs(w - w0_sub)
            return lambda_risk * portfolio_var - portfolio_ret + tc
        
        constraints = [{'type': 'eq', 'fun': lambda w: np.sum(w) - 1.0}]
        bounds = [(instance.w_min, instance.w_max)] * k
        w0 = np.ones(k) / k
        
        result = minimize(objective, w0, method='SLSQP',
                         bounds=bounds, constraints=constraints,
                         options={'ftol': 1e-9, 'maxiter': 1000})
        return result.x if result.success else w0

class BRKGA:
    """Biased Random-Key Genetic Algorithm for Portfolio Optimization"""
    
    def __init__(self, N: int, P: int = 100, rho_e: float = 0.15,
                 rho_m: float = 0.10, p_e: float = 0.70):
        self.N = N          # Number of assets (chromosome length)
        self.P = P          # Population size
        self.rho_e = rho_e  # Elite fraction
        self.rho_m = rho_m  # Mutant fraction
        self.p_e = p_e      # Elite inheritance probability
        self.n_elite = max(1, int(P * rho_e))
        self.n_mutant = max(1, int(P * rho_m))
        self.n_crossover = P - self.n_elite - self.n_mutant
        self.decoder = RKODecoder()
    
    def initialize_population(self) -> np.ndarray:
        return np.random.uniform(0, 1, (self.P, self.N))
    
    def evaluate_population(self, population: np.ndarray,
                            instance: PortfolioInstance) -> np.ndarray:
        fitness = np.zeros(self.P)
        for i in range(self.P):
            w = self.decoder.decode(population[i], instance)
            fitness[i] = self._compute_fitness(w, instance)
        return fitness
    
    def _compute_fitness(self, w: np.ndarray, 
                         instance: PortfolioInstance) -> float:
        # Negative because BRKGA minimizes
        portfolio_var = w @ instance.Sigma @ w
        portfolio_ret = instance.mu @ w
        return -(portfolio_ret - 0.5 * portfolio_var)  # Minimize negative utility
    
    def evolve(self, population: np.ndarray, 
               fitness: np.ndarray) -> np.ndarray:
        # Sort by fitness (ascending = better for minimization)
        sorted_idx = np.argsort(fitness)
        elite = population[sorted_idx[:self.n_elite]]
        
        new_pop = np.zeros_like(population)
        # Keep elite
        new_pop[:self.n_elite] = elite
        # Add mutants
        new_pop[self.n_elite:self.n_elite+self.n_mutant] = \
            np.random.uniform(0, 1, (self.n_mutant, self.N))
        # Biased crossover
        for i in range(self.n_crossover):
            elite_parent = elite[np.random.randint(self.n_elite)]
            non_elite_idx = sorted_idx[self.n_elite + 
                                       np.random.randint(self.P - self.n_elite)]
            non_elite_parent = population[non_elite_idx]
            mask = np.random.uniform(0, 1, self.N) < self.p_e
            child = np.where(mask, elite_parent, non_elite_parent)
            new_pop[self.n_elite + self.n_mutant + i] = child
        return new_pop
    
    def optimize(self, instance: PortfolioInstance, 
                 max_time: float = 3600.0,
                 max_generations: int = 10000,
                 stagnation_limit: int = 50) -> Tuple[np.ndarray, float, List]:
        import time
        population = self.initialize_population()
        fitness = self.evaluate_population(population, instance)
        
        best_fitness = np.min(fitness)
        best_weights = self.decoder.decode(
            population[np.argmin(fitness)], instance)
        convergence_curve = [(0.0, best_fitness)]
        stagnation_count = 0
        start_time = time.time()
        
        for gen in range(max_generations):
            elapsed = time.time() - start_time
            if elapsed >= max_time:
                break
            
            population = self.evolve(population, fitness)
            fitness = self.evaluate_population(population, instance)
            
            gen_best = np.min(fitness)
            if gen_best < best_fitness - 1e-8:
                best_fitness = gen_best
                best_weights = self.decoder.decode(
                    population[np.argmin(fitness)], instance)
                stagnation_count = 0
            else:
                stagnation_count += 1
            
            convergence_curve.append((elapsed, best_fitness))
            
            # Restart on stagnation
            if stagnation_count >= stagnation_limit:
                population[self.n_elite:] = np.random.uniform(
                    0, 1, (self.P - self.n_elite, self.N))
                stagnation_count = 0
        
        return best_weights, best_fitness, convergence_curve

# Benchmark comparison harness
def run_benchmark(instance: PortfolioInstance, 
                  time_limit: float = 3600.0) -> dict:
    import gurobipy as gp
    from gurobipy import GRB
    
    results = {}
    
    # --- RKO Run ---
    brkga = BRKGA(N=len(instance.mu), P=100)
    rko_weights, rko_obj, rko_curve = brkga.optimize(
        instance, max_time=time_limit)
    results['rko'] = {
        'weights': rko_weights,
        'objective': rko_obj,
        'convergence': rko_curve,
        'feasible': validate_feasibility(rko_weights, instance)
    }
    
    # --- Gurobi MIP Run ---
    N = len(instance.mu)
    m = gp.Model()
    m.setParam('TimeLimit', time_limit)
    m.setParam('MIPGap', 0.01)
    m.setParam('Threads', 4)
    
    w = m.addVars(N, lb=0, ub=instance.w_max, name='w')
    z = m.addVars(N, vtype=GRB.BINARY, name='z')
    
    # Objective: minimize variance - lambda * return
    obj = gp.quicksum(instance.Sigma[i,j] * w[i] * w[j] 
                      for i in range(N) for j in range(N))
    obj -= 2.0 * gp.quicksum(instance.mu[i] * w[i] for i in range(N))
    m.setObjective(obj, GRB.MINIMIZE)
    
    m.addConstr(gp.quicksum(w[i] for i in range(N)) == 1)
    m.addConstr(gp.quicksum(z[i] for i in range(N)) == instance.k)
    for i in range(N):
        m.addConstr(w[i] >= instance.w_min * z[i])
        m.addConstr(w[i] <= instance.w_max * z[i])
    
    m.optimize()
    results['gurobi'] = {
        'weights': np.array([w[i].X for i in range(N)]),
        'objective': m.ObjVal,
        'mip_gap': m.MIPGap,
        'solve_time': m.Runtime
    }
    return results

def validate_feasibility(weights: np.ndarray, 
                         instance: PortfolioInstance) -> bool:
    tol = 1e-4
    budget_ok = abs(np.sum(weights) - 1.0) < tol
    cardinality_ok = np.sum(weights > tol) == instance.k
    bounds_ok = all(w == 0 or (instance.w_min - tol <= w <= instance.w_max + tol)
                    for w in weights)
    return budget_ok and cardinality_ok and bounds_ok