diff --git a/Cargo.toml b/Cargo.toml
index a13676f..a3add5f 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -16,17 +16,20 @@ crate-type = ["cdylib", "rlib"]
 
 [dependencies]
 pyo3 = { version = "0.21", features = ["extension-module", "abi3-py38"] }
+numpy = "0.21"
 rand = "0.8"
 rand_distr = "0.4"
 ndarray = "0.15"
+ndarray-linalg = { version = "0.16", features = ["openblas-system"] }
 num-traits = "0.2"
-rayon = { version = "1.8", optional = true }
+rayon = "1.8"
 thiserror = "1.0"
 ordered-float = "4.2"
+statrs = "0.17"
 
 [features]
 default = []
-parallel = ["rayon"]
+parallel = []
 
 [dev-dependencies]
 criterion = "0.5"
diff --git a/python/optimizr/__init__.py b/python/optimizr/__init__.py
index 88e110b..c48ba42 100644
--- a/python/optimizr/__init__.py
+++ b/python/optimizr/__init__.py
@@ -16,6 +16,13 @@ from optimizr.core import (
     grid_search,
     mutual_information,
     shannon_entropy,
+    sparse_pca_py,
+    box_tao_decomposition_py,
+    elastic_net_py,
+    hurst_exponent_py,
+    compute_risk_metrics_py,
+    estimate_half_life_py,
+    bootstrap_returns_py,
 )
 
 __version__ = "0.1.0"
@@ -26,4 +33,11 @@ __all__ = [
     "grid_search",
     "mutual_information",
     "shannon_entropy",
+    "sparse_pca_py",
+    "box_tao_decomposition_py",
+    "elastic_net_py",
+    "hurst_exponent_py",
+    "compute_risk_metrics_py",
+    "estimate_half_life_py",
+    "bootstrap_returns_py",
 ]
diff --git a/python/optimizr/core.py b/python/optimizr/core.py
index d6be7ed..df3353f 100644
--- a/python/optimizr/core.py
+++ b/python/optimizr/core.py
@@ -14,6 +14,13 @@ try:
         grid_search as _rust_grid_search,
         mutual_information as _rust_mutual_information,
         shannon_entropy as _rust_shannon_entropy,
+        sparse_pca_py,
+        box_tao_decomposition_py,
+        elastic_net_py,
+        hurst_exponent_py,
+        compute_risk_metrics_py,
+        estimate_half_life_py,
+        bootstrap_returns_py,
     )
     RUST_AVAILABLE = True
 except ImportError:
diff --git a/src/core.rs b/src/core.rs
index 923893a..8e17606 100644
--- a/src/core.rs
+++ b/src/core.rs
@@ -11,6 +11,9 @@ pub enum OptimizrError {
     #[error("Invalid parameter: {0}")]
     InvalidParameter(String),
     
+    #[error("Invalid input: {0}")]
+    InvalidInput(String),
+    
     #[error("Dimension mismatch: expected {expected}, got {actual}")]
     DimensionMismatch { expected: usize, actual: usize },
     
@@ -30,6 +33,9 @@ pub enum OptimizrError {
 /// Result type for OptimizR operations
 pub type Result<T> = std::result::Result<T, OptimizrError>;
 
+/// Alias for OptimizR result type (legacy compatibility)
+pub type OptimizrResult<T> = Result<T>;
+
 /// Trait for optimization algorithms
 pub trait Optimizer {
     type Config;
diff --git a/src/hmm/emission.rs b/src/hmm/emission.rs
index f271f96..cad4878 100644
--- a/src/hmm/emission.rs
+++ b/src/hmm/emission.rs
@@ -3,7 +3,7 @@
 //! This module defines the EmissionModel trait and provides
 //! implementations for different probability distributions.
 
-use crate::core::{OptimizrError, Result};
+use crate::core::Result;
 use std::f64;
 
 /// Trait for emission probability models
diff --git a/src/hmm/viterbi.rs b/src/hmm/viterbi.rs
index 46e726c..ba1efb2 100644
--- a/src/hmm/viterbi.rs
+++ b/src/hmm/viterbi.rs
@@ -2,7 +2,7 @@
 //!
 //! Implements the Viterbi algorithm for HMM inference.
 
-use super::config::HMMConfig;
+
 use super::emission::EmissionModel;
 use super::model::HMM;
 use crate::core::Result;
diff --git a/src/lib.rs b/src/lib.rs
index 8aacb46..08caf22 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -21,6 +21,8 @@
 //! - `differential_evolution`: Global optimization algorithm
 //! - `grid_search`: Exhaustive parameter space search
 //! - `information_theory`: Mutual information and entropy calculations
+//! - `sparse_optimization`: Sparse PCA, Box-Tao, Elastic Net
+//! - `risk_metrics`: Portfolio risk analysis and Hurst exponent
 
 use pyo3::prelude::*;
 use pyo3::types::PyModule;
@@ -33,6 +35,8 @@ pub mod functional;
 pub mod hmm;
 pub mod mcmc;
 pub mod de;
+pub mod sparse_optimization;
+pub mod risk_metrics;
 
 // Legacy modules for backward compatibility
 mod hmm_legacy;
@@ -72,5 +76,18 @@ fn _core(_py: Python, m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(information_theory::mutual_information, m)?)?;
     m.add_function(wrap_pyfunction!(information_theory::shannon_entropy, m)?)?;
     
+    // ===== New Optimization Algorithms =====
+    
+    // Sparse optimization functions
+    m.add_function(wrap_pyfunction!(sparse_optimization::sparse_pca_py, m)?)?;
+    m.add_function(wrap_pyfunction!(sparse_optimization::box_tao_decomposition_py, m)?)?;
+    m.add_function(wrap_pyfunction!(sparse_optimization::elastic_net_py, m)?)?;
+    
+    // Risk metrics functions
+    m.add_function(wrap_pyfunction!(risk_metrics::hurst_exponent_py, m)?)?;
+    m.add_function(wrap_pyfunction!(risk_metrics::compute_risk_metrics_py, m)?)?;
+    m.add_function(wrap_pyfunction!(risk_metrics::estimate_half_life_py, m)?)?;
+    m.add_function(wrap_pyfunction!(risk_metrics::bootstrap_returns_py, m)?)?;
+    
     Ok(())
 }
diff --git a/src/risk_metrics.rs b/src/risk_metrics.rs
new file mode 100644
index 0000000..80bd0f2
--- /dev/null
+++ b/src/risk_metrics.rs
@@ -0,0 +1,577 @@
+//! Risk Metrics and Portfolio Analysis
+//!
+//! Comprehensive risk metrics for portfolio evaluation:
+//! - Hurst exponent (mean-reversion detection)
+//! - Sharpe, Sortino, Calmar ratios
+//! - VaR and CVaR (Expected Shortfall)
+//! - Maximum Drawdown and recovery metrics
+//! - Information Ratio and tracking error
+//! - Monte Carlo simulation utilities
+
+use ndarray::{Array1, Array2};
+use crate::core::{OptimizrError, OptimizrResult};
+use rayon::prelude::*;
+
+/// Result from Hurst exponent calculation
+#[derive(Debug, Clone)]
+pub struct HurstResult {
+    pub hurst_exponent: f64,           // H < 0.5 = mean-reverting
+    pub confidence_interval: (f64, f64), // 95% CI
+    pub standard_error: f64,           // SE of estimate
+    pub is_mean_reverting: bool,       // True if CI upper < 0.5
+    pub window_sizes: Vec<usize>,      // Window sizes used
+    pub rs_values: Vec<f64>,           // R/S values
+}
+
+/// Hurst Exponent via Rescaled Range (R/S) Analysis
+///
+/// Estimates the Hurst exponent H:
+/// - H = 0.5: Random walk (Brownian motion)
+/// - H < 0.5: Mean-reverting (anti-persistent)
+/// - H > 0.5: Trending (persistent)
+///
+/// # Arguments
+/// * `series` - Time series data
+/// * `window_sizes` - Window sizes for R/S analysis (e.g., [8, 16, 32, 64])
+///
+/// # Returns
+/// `HurstResult` with Hurst exponent and statistics
+pub fn hurst_exponent(
+    series: &Array1<f64>,
+    window_sizes: &[usize],
+) -> OptimizrResult<HurstResult> {
+    let n = series.len();
+    
+    if n < 20 {
+        return Err(OptimizrError::InvalidInput(
+            "Series too short for Hurst analysis (need >= 20)".to_string()
+        ));
+    }
+    
+    if window_sizes.is_empty() {
+        return Err(OptimizrError::InvalidInput(
+            "window_sizes cannot be empty".to_string()
+        ));
+    }
+    
+    let mut rs_values = Vec::with_capacity(window_sizes.len());
+    let mut log_lags = Vec::with_capacity(window_sizes.len());
+    
+    for &lag in window_sizes {
+        if lag >= n {
+            continue;
+        }
+        
+        let n_windows = n / lag;
+        let mut rs_window_values = Vec::with_capacity(n_windows);
+        
+        for w in 0..n_windows {
+            let start = w * lag;
+            let end = start + lag;
+            let window = series.slice(ndarray::s![start..end]);
+            
+            // Compute mean
+            let mean = window.mean().unwrap_or(0.0);
+            
+            // Cumulative deviations from mean
+            let mut y = Vec::with_capacity(lag);
+            let mut cumsum = 0.0;
+            for &x in window.iter() {
+                cumsum += x - mean;
+                y.push(cumsum);
+            }
+            
+            // Range R
+            let max_y = y.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
+            let min_y = y.iter().cloned().fold(f64::INFINITY, f64::min);
+            let range = max_y - min_y;
+            
+            // Standard deviation S
+            let variance: f64 = window.iter()
+                .map(|&x| (x - mean).powi(2))
+                .sum::<f64>() / lag as f64;
+            let std = variance.sqrt();
+            
+            if std > 1e-10 {
+                rs_window_values.push(range / std);
+            }
+        }
+        
+        if !rs_window_values.is_empty() {
+            let mean_rs: f64 = rs_window_values.iter().sum::<f64>() / rs_window_values.len() as f64;
+            rs_values.push(mean_rs);
+            log_lags.push((lag as f64).ln());
+        }
+    }
+    
+    if rs_values.len() < 2 {
+        return Err(OptimizrError::ComputationError(
+            "Not enough valid window sizes for regression".to_string()
+        ));
+    }
+    
+    // Linear regression: log(R/S) = H * log(lag) + c
+    let log_rs: Vec<f64> = rs_values.iter().map(|x| x.ln()).collect();
+    let (slope, _intercept, se) = linear_regression(&log_lags, &log_rs)?;
+    
+    let hurst = slope;
+    
+    // 95% confidence interval
+    let t_critical = 1.96; // Approximate for large n
+    let ci_lower = hurst - t_critical * se;
+    let ci_upper = hurst + t_critical * se;
+    
+    let is_mean_reverting = ci_upper < 0.5;
+    
+    Ok(HurstResult {
+        hurst_exponent: hurst,
+        confidence_interval: (ci_lower, ci_upper),
+        standard_error: se,
+        is_mean_reverting,
+        window_sizes: window_sizes.to_vec(),
+        rs_values,
+    })
+}
+
+/// Simple linear regression
+fn linear_regression(x: &[f64], y: &[f64]) -> OptimizrResult<(f64, f64, f64)> {
+    if x.len() != y.len() || x.is_empty() {
+        return Err(OptimizrError::InvalidInput(
+            "x and y must have same non-zero length".to_string()
+        ));
+    }
+    
+    let n = x.len() as f64;
+    let x_mean: f64 = x.iter().sum::<f64>() / n;
+    let y_mean: f64 = y.iter().sum::<f64>() / n;
+    
+    let mut numerator = 0.0;
+    let mut denominator = 0.0;
+    
+    for i in 0..x.len() {
+        let x_dev = x[i] - x_mean;
+        let y_dev = y[i] - y_mean;
+        numerator += x_dev * y_dev;
+        denominator += x_dev * x_dev;
+    }
+    
+    if denominator.abs() < 1e-10 {
+        return Err(OptimizrError::ComputationError(
+            "Degenerate regression (zero variance in x)".to_string()
+        ));
+    }
+    
+    let slope = numerator / denominator;
+    let intercept = y_mean - slope * x_mean;
+    
+    // Standard error of slope
+    let mut residual_sum = 0.0;
+    for i in 0..x.len() {
+        let predicted = slope * x[i] + intercept;
+        residual_sum += (y[i] - predicted).powi(2);
+    }
+    
+    let mse = residual_sum / (n - 2.0).max(1.0);
+    let se = (mse / denominator).sqrt();
+    
+    Ok((slope, intercept, se))
+}
+
+/// Portfolio risk metrics
+#[derive(Debug, Clone)]
+pub struct RiskMetrics {
+    pub sharpe_ratio: f64,          // Risk-adjusted return
+    pub sortino_ratio: f64,         // Downside risk-adjusted return
+    pub calmar_ratio: f64,          // Return / Max Drawdown
+    pub max_drawdown: f64,          // Maximum peak-to-trough decline
+    pub max_drawdown_duration: usize, // Days in max drawdown
+    pub var_95: f64,                // Value at Risk (95%)
+    pub cvar_95: f64,               // Conditional VaR (Expected Shortfall)
+    pub volatility: f64,            // Annualized volatility
+    pub downside_deviation: f64,    // Downside risk
+    pub skewness: f64,              // Return distribution skewness
+    pub kurtosis: f64,              // Return distribution kurtosis
+}
+
+/// Compute comprehensive risk metrics for a return series
+///
+/// # Arguments
+/// * `returns` - Daily returns (not log returns)
+/// * `risk_free_rate` - Annual risk-free rate (e.g., 0.02 for 2%)
+/// * `periods_per_year` - 252 for daily, 52 for weekly, 12 for monthly
+///
+/// # Returns
+/// `RiskMetrics` with all computed metrics
+pub fn compute_risk_metrics(
+    returns: &Array1<f64>,
+    risk_free_rate: f64,
+    periods_per_year: f64,
+) -> OptimizrResult<RiskMetrics> {
+    if returns.is_empty() {
+        return Err(OptimizrError::InvalidInput("Returns array is empty".to_string()));
+    }
+    
+    let n = returns.len() as f64;
+    let daily_rf = risk_free_rate / periods_per_year;
+    
+    // Mean return
+    let mean_return = returns.mean().unwrap_or(0.0);
+    let excess_return = mean_return - daily_rf;
+    
+    // Volatility
+    let variance = returns.iter()
+        .map(|&r| (r - mean_return).powi(2))
+        .sum::<f64>() / n;
+    let volatility = variance.sqrt() * periods_per_year.sqrt();
+    
+    // Sharpe ratio
+    let sharpe_ratio = if volatility > 1e-10 {
+        excess_return * periods_per_year.sqrt() / volatility
+    } else {
+        0.0
+    };
+    
+    // Downside deviation (semi-deviation)
+    let downside_returns: Vec<f64> = returns.iter()
+        .filter(|&&r| r < 0.0)
+        .copied()
+        .collect();
+    
+    let downside_deviation = if !downside_returns.is_empty() {
+        let dd_variance = downside_returns.iter()
+            .map(|&r| r.powi(2))
+            .sum::<f64>() / downside_returns.len() as f64;
+        dd_variance.sqrt() * periods_per_year.sqrt()
+    } else {
+        0.0
+    };
+    
+    // Sortino ratio
+    let sortino_ratio = if downside_deviation > 1e-10 {
+        excess_return * periods_per_year.sqrt() / downside_deviation
+    } else {
+        0.0
+    };
+    
+    // Maximum drawdown
+    let cum_returns = returns.iter()
+        .scan(1.0, |state, &r| {
+            *state *= 1.0 + r;
+            Some(*state)
+        })
+        .collect::<Vec<f64>>();
+    
+    let mut running_max = cum_returns[0];
+    let mut max_dd: f64 = 0.0;
+    let mut current_dd_duration = 0;
+    let mut max_dd_duration = 0;
+    let mut in_drawdown = false;
+    
+    for &cum_ret in &cum_returns {
+        if cum_ret > running_max {
+            running_max = cum_ret;
+            if in_drawdown {
+                max_dd_duration = max_dd_duration.max(current_dd_duration);
+                current_dd_duration = 0;
+                in_drawdown = false;
+            }
+        } else {
+            let dd = (running_max - cum_ret) / running_max;
+            max_dd = max_dd.max(dd);
+            in_drawdown = true;
+            current_dd_duration += 1;
+        }
+    }
+    
+    if in_drawdown {
+        max_dd_duration = max_dd_duration.max(current_dd_duration);
+    }
+    
+    // Calmar ratio
+    let annual_return = mean_return * periods_per_year;
+    let calmar_ratio = if max_dd > 1e-10 {
+        annual_return / max_dd
+    } else {
+        0.0
+    };
+    
+    // VaR and CVaR (95% confidence)
+    let mut sorted_returns = returns.to_vec();
+    sorted_returns.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
+    
+    let var_95_idx = (n * 0.05) as usize;
+    let var_95 = -sorted_returns[var_95_idx.min(sorted_returns.len() - 1)];
+    
+    let cvar_95 = if var_95_idx > 0 {
+        -sorted_returns[..var_95_idx].iter().sum::<f64>() / var_95_idx as f64
+    } else {
+        var_95
+    };
+    
+    // Skewness and kurtosis
+    let skewness = if n > 2.0 {
+        let m3: f64 = returns.iter()
+            .map(|&r| (r - mean_return).powi(3))
+            .sum::<f64>() / n;
+        m3 / variance.powf(1.5)
+    } else {
+        0.0
+    };
+    
+    let kurtosis = if n > 3.0 {
+        let m4: f64 = returns.iter()
+            .map(|&r| (r - mean_return).powi(4))
+            .sum::<f64>() / n;
+        m4 / variance.powi(2) - 3.0  // Excess kurtosis
+    } else {
+        0.0
+    };
+    
+    Ok(RiskMetrics {
+        sharpe_ratio,
+        sortino_ratio,
+        calmar_ratio,
+        max_drawdown: max_dd,
+        max_drawdown_duration: max_dd_duration,
+        var_95,
+        cvar_95,
+        volatility,
+        downside_deviation,
+        skewness,
+        kurtosis,
+    })
+}
+
+/// Half-life estimation for mean-reverting series
+///
+/// Estimates the half-life of mean reversion using AR(1) model:
+/// ΔX_t = θ X_{t-1} + ε_t
+///
+/// Half-life = -ln(2) / ln(|θ|)
+///
+/// # Arguments
+/// * `series` - Time series (typically portfolio value or spread)
+///
+/// # Returns
+/// Half-life in time periods (e.g., days if daily data)
+pub fn estimate_half_life(series: &Array1<f64>) -> OptimizrResult<f64> {
+    if series.len() < 10 {
+        return Err(OptimizrError::InvalidInput(
+            "Series too short for half-life estimation".to_string()
+        ));
+    }
+    
+    let mean = series.mean().unwrap_or(0.0);
+    let deviations: Vec<f64> = series.iter().map(|&x| x - mean).collect();
+    
+    // AR(1) regression: dev[t] = theta * dev[t-1] + error
+    let x: Vec<f64> = deviations[..deviations.len()-1].to_vec();
+    let y: Vec<f64> = deviations[1..].to_vec();
+    
+    let (theta, _, _) = linear_regression(&x, &y)?;
+    
+    if theta.abs() >= 1.0 || theta.abs() < 1e-10 {
+        // Not mean-reverting or degenerate
+        return Ok(f64::INFINITY);
+    }
+    
+    let half_life = -2.0_f64.ln() / theta.abs().ln();
+    
+    Ok(half_life)
+}
+
+/// Monte Carlo bootstrap simulation
+///
+/// Performs bootstrap resampling to estimate uncertainty in metrics
+///
+/// # Arguments
+/// * `returns` - Historical returns
+/// * `n_simulations` - Number of bootstrap samples
+/// * `block_size` - Optional block size for block bootstrap (preserves autocorrelation)
+///
+/// # Returns
+/// Vector of simulated return series
+pub fn bootstrap_returns(
+    returns: &Array1<f64>,
+    n_simulations: usize,
+    block_size: Option<usize>,
+) -> Vec<Array1<f64>> {
+    use rand::Rng;
+    use rand::SeedableRng;
+    
+    let n = returns.len();
+    let block_size = block_size.unwrap_or(1);
+    
+    (0..n_simulations)
+        .into_par_iter()
+        .map(|seed| {
+            let mut rng = rand::rngs::StdRng::seed_from_u64(seed as u64);
+            let mut sim_returns = Vec::with_capacity(n);
+            
+            if block_size == 1 {
+                // Simple bootstrap
+                while sim_returns.len() < n {
+                    let idx = rng.gen_range(0..n);
+                    sim_returns.push(returns[idx]);
+                }
+            } else {
+                // Block bootstrap
+                let n_blocks = (n as f64 / block_size as f64).ceil() as usize;
+                
+                for _ in 0..n_blocks {
+                    let start_idx = rng.gen_range(0..=(n.saturating_sub(block_size)));
+                    let end_idx = (start_idx + block_size).min(n);
+                    
+                    for idx in start_idx..end_idx {
+                        if sim_returns.len() < n {
+                            sim_returns.push(returns[idx]);
+                        }
+                    }
+                }
+            }
+            
+            Array1::from_vec(sim_returns[..n].to_vec())
+        })
+        .collect()
+}
+
+// ============================================================================
+// Python Bindings
+// ============================================================================
+
+use pyo3::prelude::*;
+use pyo3::exceptions::PyValueError;
+use pyo3::types::PyDict;
+use numpy::{PyArray2, PyReadonlyArray1};
+
+#[pyfunction]
+#[pyo3(signature = (series, window_sizes=None))]
+pub fn hurst_exponent_py(
+    py: Python,
+    series: PyReadonlyArray1<f64>,
+    window_sizes: Option<Vec<usize>>,
+) -> PyResult<PyObject> {
+    let series_arr = series.as_array().to_owned();
+    let windows = window_sizes.unwrap_or_else(|| vec![8, 16, 32, 64, 128]);
+    
+    let result = hurst_exponent(&series_arr, &windows)
+        .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+    
+    let dict = PyDict::new_bound(py);
+    dict.set_item("hurst_exponent", result.hurst_exponent)?;
+    dict.set_item("confidence_interval", result.confidence_interval)?;
+    dict.set_item("standard_error", result.standard_error)?;
+    dict.set_item("is_mean_reverting", result.is_mean_reverting)?;
+    dict.set_item("window_sizes", result.window_sizes)?;
+    dict.set_item("rs_values", result.rs_values)?;
+    
+    Ok(dict.into())
+}
+
+#[pyfunction]
+#[pyo3(signature = (returns, risk_free_rate=0.02, periods_per_year=252.0))]
+pub fn compute_risk_metrics_py(
+    py: Python,
+    returns: PyReadonlyArray1<f64>,
+    risk_free_rate: f64,
+    periods_per_year: f64,
+) -> PyResult<PyObject> {
+    let returns_arr = returns.as_array().to_owned();
+    
+    let result = compute_risk_metrics(&returns_arr, risk_free_rate, periods_per_year)
+        .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+    
+    let dict = PyDict::new_bound(py);
+    dict.set_item("sharpe_ratio", result.sharpe_ratio)?;
+    dict.set_item("sortino_ratio", result.sortino_ratio)?;
+    dict.set_item("calmar_ratio", result.calmar_ratio)?;
+    dict.set_item("max_drawdown", result.max_drawdown)?;
+    dict.set_item("max_drawdown_duration", result.max_drawdown_duration)?;
+    dict.set_item("var_95", result.var_95)?;
+    dict.set_item("cvar_95", result.cvar_95)?;
+    dict.set_item("volatility", result.volatility)?;
+    dict.set_item("downside_deviation", result.downside_deviation)?;
+    dict.set_item("skewness", result.skewness)?;
+    dict.set_item("kurtosis", result.kurtosis)?;
+    
+    Ok(dict.into())
+}
+
+#[pyfunction]
+pub fn estimate_half_life_py(
+    series: PyReadonlyArray1<f64>,
+) -> PyResult<f64> {
+    let series_arr = series.as_array().to_owned();
+    
+    estimate_half_life(&series_arr)
+        .map_err(|e| PyValueError::new_err(format!("{}", e)))
+}
+
+#[pyfunction]
+#[pyo3(signature = (returns, n_simulations=1000, block_size=None))]
+pub fn bootstrap_returns_py(
+    py: Python,
+    returns: PyReadonlyArray1<f64>,
+    n_simulations: usize,
+    block_size: Option<usize>,
+) -> PyResult<PyObject> {
+    let returns_arr = returns.as_array().to_owned();
+    
+    let simulations = bootstrap_returns(&returns_arr, n_simulations, block_size);
+    
+    // Convert Vec<Array1> to 2D array
+    let n = returns_arr.len();
+    let mut result = Array2::zeros((n_simulations, n));
+    
+    for (i, sim) in simulations.iter().enumerate() {
+        for (j, &val) in sim.iter().enumerate() {
+            result[[i, j]] = val;
+        }
+    }
+    
+    Ok(PyArray2::from_owned_array_bound(py, result).into())
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use approx::assert_relative_eq;
+    
+    #[test]
+    fn test_hurst_random_walk() {
+        // Random walk should have H ≈ 0.5
+        let n = 1000;
+        let mut series = vec![0.0];
+        for i in 1..n {
+            series.push(series[i-1] + if i % 2 == 0 { 1.0 } else { -1.0 });
+        }
+        
+        let series = Array1::from_vec(series);
+        let result = hurst_exponent(&series, &[8, 16, 32, 64]).unwrap();
+        
+        // Should be close to 0.5
+        assert!((result.hurst_exponent - 0.5).abs() < 0.2);
+    }
+    
+    #[test]
+    fn test_risk_metrics() {
+        // Simple return series
+        let returns = Array1::from_vec(vec![0.01, -0.005, 0.015, -0.01, 0.02]);
+        let metrics = compute_risk_metrics(&returns, 0.02, 252.0).unwrap();
+        
+        assert!(metrics.volatility > 0.0);
+        assert!(metrics.max_drawdown >= 0.0);
+        assert!(metrics.var_95 >= 0.0);
+    }
+    
+    #[test]
+    fn test_half_life() {
+        // Mean-reverting series
+        let series = Array1::from_vec(vec![
+            1.0, 0.8, 0.6, 0.7, 0.9, 1.0, 1.1, 0.9, 1.0, 0.95
+        ]);
+        
+        let half_life = estimate_half_life(&series).unwrap();
+        assert!(half_life > 0.0 && half_life < 100.0);
+    }
+}
diff --git a/src/sparse_optimization.rs b/src/sparse_optimization.rs
new file mode 100644
index 0000000..9322bb8
--- /dev/null
+++ b/src/sparse_optimization.rs
@@ -0,0 +1,573 @@
+//! Sparse Optimization Algorithms
+//! 
+//! Generic implementations of sparse optimization and decomposition methods:
+//! - Sparse PCA with L1 regularization
+//! - Box & Tao decomposition (Robust PCA)
+//! - Elastic Net for sparse linear models
+//! - ADMM (Alternating Direction Method of Multipliers)
+//! 
+//! Based on:
+//! - d'Aspremont (2011): "Identifying Small Mean Reverting Portfolios"
+//! - Candès et al. (2011): "Robust Principal Component Analysis?"
+//! - Zou & Hastie (2005): "Regularization and Variable Selection via Elastic Net"
+
+use ndarray::{Array1, Array2, Axis};
+use ndarray_linalg::{SVD, Norm};
+use crate::core::{OptimizrError, OptimizrResult};
+
+/// Soft thresholding operator
+/// 
+/// S_λ(x) = sign(x) * max(|x| - λ, 0)
+#[inline]
+pub fn soft_threshold(x: f64, lambda: f64) -> f64 {
+    let abs_x = x.abs();
+    if abs_x <= lambda {
+        0.0
+    } else {
+        x.signum() * (abs_x - lambda)
+    }
+}
+
+/// Vectorized soft thresholding
+pub fn soft_threshold_array(arr: &Array1<f64>, lambda: f64) -> Array1<f64> {
+    arr.mapv(|x| soft_threshold(x, lambda))
+}
+
+/// Soft thresholding for matrices
+pub fn soft_threshold_matrix(mat: &Array2<f64>, lambda: f64) -> Array2<f64> {
+    mat.mapv(|x| soft_threshold(x, lambda))
+}
+
+/// SVD soft thresholding for nuclear norm regularization
+/// 
+/// Returns U * S_λ(Σ) * V^T where S_λ is soft thresholding on singular values
+pub fn svd_soft_threshold(mat: &Array2<f64>, lambda: f64) -> OptimizrResult<Array2<f64>> {
+    let (u, s, vt) = mat.svd(true, true)
+        .map_err(|e| OptimizrError::ComputationError(format!("SVD failed: {}", e)))?;
+    
+    let u = u.ok_or_else(|| OptimizrError::ComputationError("SVD U is None".to_string()))?;
+    let vt = vt.ok_or_else(|| OptimizrError::ComputationError("SVD Vt is None".to_string()))?;
+    
+    // Soft threshold singular values
+    let s_thresh = s.mapv(|val| soft_threshold(val, lambda));
+    
+    // Reconstruct: U * diag(s_thresh) * V^T
+    let s_diag = Array2::from_diag(&s_thresh);
+    let us = u.dot(&s_diag);
+    Ok(us.dot(&vt))
+}
+
+/// Result from Sparse PCA
+#[derive(Debug, Clone)]
+pub struct SparsePCAResult {
+    pub weights: Array2<f64>,              // (n_components, n_features) sparse weights
+    pub variance_explained: Array1<f64>,   // Variance explained per component
+    pub sparsity: Array1<f64>,             // Sparsity level per component
+    pub iterations: Array1<usize>,         // Iterations to converge per component
+    pub converged: Vec<bool>,              // Convergence flag per component
+}
+
+/// Sparse Principal Component Analysis with L1 regularization
+/// 
+/// Finds sparse principal components that maximize variance with sparsity constraint:
+/// 
+/// max_w  w^T Σ w - λ ||w||_1  s.t. ||w||_2 = 1
+/// 
+/// # Arguments
+/// * `covariance` - Covariance matrix (n_features, n_features)
+/// * `n_components` - Number of sparse components to extract
+/// * `lambda` - Sparsity parameter (larger = sparser)
+/// * `max_iter` - Maximum iterations per component
+/// * `tol` - Convergence tolerance
+/// 
+/// # Returns
+/// `SparsePCAResult` with sparse principal components
+pub fn sparse_pca(
+    covariance: &Array2<f64>,
+    n_components: usize,
+    lambda: f64,
+    max_iter: usize,
+    tol: f64,
+) -> OptimizrResult<SparsePCAResult> {
+    let n_features = covariance.nrows();
+    
+    if covariance.ncols() != n_features {
+        return Err(OptimizrError::InvalidInput(
+            "Covariance matrix must be square".to_string()
+        ));
+    }
+    
+    if n_components > n_features {
+        return Err(OptimizrError::InvalidInput(
+            format!("n_components ({}) > n_features ({})", n_components, n_features)
+        ));
+    }
+    
+    let mut weights = Array2::zeros((n_components, n_features));
+    let mut variance_explained = Array1::zeros(n_components);
+    let mut sparsity = Array1::zeros(n_components);
+    let mut iterations_vec = Array1::zeros(n_components);
+    let mut converged_vec = vec![false; n_components];
+    
+    let mut residual_cov = covariance.clone();
+    
+    for comp in 0..n_components {
+        // Initialize with leading eigenvector
+        let (_, s, vt) = residual_cov.svd(false, true)
+            .map_err(|e| OptimizrError::ComputationError(format!("SVD failed: {}", e)))?;
+        
+        let vt = vt.ok_or_else(|| OptimizrError::ComputationError("SVD Vt is None".to_string()))?;
+        let mut w = vt.row(0).to_owned();
+        
+        let mut converged = false;
+        let mut iter = 0;
+        
+        // Iterative soft-thresholding
+        for _ in 0..max_iter {
+            let w_old = w.clone();
+            
+            // Update: w_new = Σ * w
+            let mut w_new = residual_cov.dot(&w);
+            
+            // Apply soft thresholding
+            w_new = soft_threshold_array(&w_new, lambda);
+            
+            // Normalize
+            let norm = w_new.norm_l2();
+            if norm > 1e-10 {
+                w_new /= norm;
+            } else {
+                // Degenerate case - use previous
+                w_new = w_old.clone();
+                break;
+            }
+            
+            // Check convergence
+            let diff = (&w_new - &w_old).norm_l2();
+            if diff < tol {
+                converged = true;
+                w = w_new;
+                break;
+            }
+            
+            w = w_new;
+            iter += 1;
+        }
+        
+        // Store results for this component
+        weights.row_mut(comp).assign(&w);
+        
+        // Variance explained
+        let var_exp = w.dot(&residual_cov.dot(&w));
+        variance_explained[comp] = var_exp.max(0.0);
+        
+        // Sparsity (proportion of non-zero elements)
+        let non_zero = w.iter().filter(|&&x| x.abs() > 1e-10).count();
+        sparsity[comp] = 1.0 - (non_zero as f64 / n_features as f64);
+        
+        iterations_vec[comp] = iter as f64;
+        converged_vec[comp] = converged;
+        
+        // Deflate covariance matrix
+        let w_outer = outer_product(&w, &w);
+        residual_cov = &residual_cov - &(w_outer * var_exp);
+    }
+    
+    Ok(SparsePCAResult {
+        weights,
+        variance_explained,
+        sparsity,
+        iterations: iterations_vec.mapv(|x| x as usize),
+        converged: converged_vec,
+    })
+}
+
+/// Helper: outer product of two vectors
+fn outer_product(a: &Array1<f64>, b: &Array1<f64>) -> Array2<f64> {
+    let n = a.len();
+    let m = b.len();
+    let mut result = Array2::zeros((n, m));
+    
+    for i in 0..n {
+        for j in 0..m {
+            result[[i, j]] = a[i] * b[j];
+        }
+    }
+    
+    result
+}
+
+/// Result from Box & Tao Decomposition
+#[derive(Debug, Clone)]
+pub struct BoxTaoResult {
+    pub low_rank: Array2<f64>,       // Low-rank component (common factors)
+    pub sparse: Array2<f64>,         // Sparse component (idiosyncratic)
+    pub noise: Array2<f64>,          // Noise/residual
+    pub rank: usize,                 // Rank of low-rank component
+    pub sparsity: f64,               // Sparsity of sparse component
+    pub iterations: usize,           // Number of ADMM iterations
+    pub converged: bool,             // Convergence flag
+    pub objective_values: Vec<f64>,  // Objective function history
+}
+
+/// Box & Tao Decomposition (Robust PCA)
+/// 
+/// Decomposes matrix into low-rank + sparse + noise:
+/// 
+/// X = L + S + N
+/// 
+/// min_{L,S}  ||L||_* + λ ||S||_1  s.t. ||X - L - S||_F ≤ ε
+/// 
+/// Solved via ADMM (Alternating Direction Method of Multipliers)
+/// 
+/// # Arguments
+/// * `matrix` - Input matrix to decompose
+/// * `lambda` - Sparsity parameter for sparse component
+/// * `mu` - Penalty parameter for ADMM
+/// * `max_iter` - Maximum ADMM iterations
+/// * `tol` - Convergence tolerance
+/// 
+/// # Returns
+/// `BoxTaoResult` with decomposed components
+pub fn box_tao_decomposition(
+    matrix: &Array2<f64>,
+    lambda: f64,
+    mu: f64,
+    max_iter: usize,
+    tol: f64,
+) -> OptimizrResult<BoxTaoResult> {
+    let (m, n) = matrix.dim();
+    
+    // Initialize L, S, Y (dual variable)
+    let mut low_rank = Array2::<f64>::zeros((m, n));
+    let mut sparse = Array2::<f64>::zeros((m, n));
+    let mut dual = Array2::<f64>::zeros((m, n));
+    
+    let mut objective_values = Vec::with_capacity(max_iter);
+    let mut converged = false;
+    let mut iter = 0;
+    
+    let rho = mu;  // ADMM penalty parameter
+    
+    for _ in 0..max_iter {
+        // Update L: SVD soft-thresholding
+        let l_update = matrix - &sparse + &(&dual / rho);
+        low_rank = svd_soft_threshold(&l_update, 1.0 / rho)?;
+        
+        // Update S: Element-wise soft-thresholding  
+        let s_update = matrix - &low_rank + &(&dual / rho);
+        sparse = soft_threshold_matrix(&s_update, lambda / rho);
+        
+        // Update dual variable Y
+        let residual = matrix - &low_rank - &sparse;
+        let primal_residual = residual.norm_l2();
+        dual = &dual + &(residual * rho);
+        
+        // Compute objective
+        let nuclear_norm = low_rank.svd(false, false)
+            .map(|(_, s, _)| s.sum())
+            .unwrap_or(0.0);
+        let l1_norm = sparse.iter().map(|x| x.abs()).sum::<f64>();
+        let objective = nuclear_norm + lambda * l1_norm;
+        objective_values.push(objective);
+        
+        // Check convergence
+        let dual_residual = if iter > 0 {
+            rho * (&sparse - matrix).norm_l2()
+        } else {
+            f64::INFINITY
+        };
+        
+        if primal_residual < tol && dual_residual < tol {
+            converged = true;
+            break;
+        }
+        
+        iter += 1;
+    }
+    
+    let noise = matrix - &low_rank - &sparse;
+    
+    // Compute rank of low-rank component
+    let rank = low_rank.svd(false, false)
+        .map(|(_, s, _)| s.iter().filter(|&&x| x > 1e-10).count())
+        .unwrap_or(0);
+    
+    // Compute sparsity of sparse component
+    let total_elements = (m * n) as f64;
+    let non_zero = sparse.iter().filter(|&&x| x.abs() > 1e-10).count() as f64;
+    let sparsity = 1.0 - (non_zero / total_elements);
+    
+    Ok(BoxTaoResult {
+        low_rank,
+        sparse,
+        noise,
+        rank,
+        sparsity,
+        iterations: iter,
+        converged,
+        objective_values,
+    })
+}
+
+/// Result from Elastic Net regression
+#[derive(Debug, Clone)]
+pub struct ElasticNetResult {
+    pub weights: Array1<f64>,        // Sparse regression coefficients
+    pub intercept: f64,              // Intercept term
+    pub sparsity: f64,               // Sparsity level
+    pub iterations: usize,           // Iterations to converge
+    pub converged: bool,             // Convergence flag
+    pub objective_value: f64,        // Final objective value
+}
+
+/// Elastic Net Regression
+/// 
+/// Sparse linear regression with L1 + L2 regularization:
+/// 
+/// min_w  (1/2n) ||y - Xw||_2^2 + λ₁ ||w||_1 + (λ₂/2) ||w||_2^2
+/// 
+/// Combines LASSO (L1) sparsity with Ridge (L2) smoothness
+/// 
+/// # Arguments
+/// * `x` - Feature matrix (n_samples, n_features)
+/// * `y` - Target vector (n_samples,)
+/// * `lambda_l1` - L1 regularization (sparsity)
+/// * `lambda_l2` - L2 regularization (smoothness)
+/// * `max_iter` - Maximum iterations
+/// * `tol` - Convergence tolerance
+/// 
+/// # Returns
+/// `ElasticNetResult` with sparse coefficients
+pub fn elastic_net(
+    x: &Array2<f64>,
+    y: &Array1<f64>,
+    lambda_l1: f64,
+    lambda_l2: f64,
+    max_iter: usize,
+    tol: f64,
+) -> OptimizrResult<ElasticNetResult> {
+    let (n_samples, n_features) = x.dim();
+    
+    if y.len() != n_samples {
+        return Err(OptimizrError::InvalidInput(
+            format!("y length ({}) != n_samples ({})", y.len(), n_samples)
+        ));
+    }
+    
+    // Center data
+    let x_mean = x.mean_axis(Axis(0)).unwrap();
+    let y_mean = y.mean().unwrap();
+    
+    let x_centered = x - &x_mean;
+    let y_centered = y - y_mean;
+    
+    // Initialize weights
+    let mut w = Array1::zeros(n_features);
+    let mut converged = false;
+    let mut iter = 0;
+    
+    // Coordinate descent
+    for _ in 0..max_iter {
+        let w_old = w.clone();
+        
+        for j in 0..n_features {
+            // Compute residual excluding feature j
+            let mut residual = y_centered.clone();
+            for k in 0..n_features {
+                if k != j {
+                    let x_k = x_centered.column(k);
+                    residual = &residual - &(&x_k * w[k]);
+                }
+            }
+            
+            // Update weight j
+            let x_j = x_centered.column(j);
+            let rho = x_j.dot(&residual) / n_samples as f64;
+            let z_j = x_j.dot(&x_j) / n_samples as f64 + lambda_l2;
+            
+            w[j] = soft_threshold(rho, lambda_l1) / z_j;
+        }
+        
+        // Check convergence
+        let diff = (&w - &w_old).norm_l2();
+        if diff < tol {
+            converged = true;
+            break;
+        }
+        
+        iter += 1;
+    }
+    
+    // Compute intercept
+    let intercept = y_mean - x_mean.dot(&w);
+    
+    // Compute sparsity
+    let non_zero = w.iter().filter(|&&x| x.abs() > 1e-10).count() as f64;
+    let sparsity = 1.0 - (non_zero / n_features as f64);
+    
+    // Compute objective value
+    let predictions = &x_centered.dot(&w) + y_mean;
+    let residuals = y - &predictions;
+    let mse = residuals.dot(&residuals) / (2.0 * n_samples as f64);
+    let l1_penalty = lambda_l1 * w.iter().map(|x| x.abs()).sum::<f64>();
+    let l2_penalty = 0.5 * lambda_l2 * w.dot(&w);
+    let objective_value = mse + l1_penalty + l2_penalty;
+    
+    Ok(ElasticNetResult {
+        weights: w,
+        intercept,
+        sparsity,
+        iterations: iter,
+        converged,
+        objective_value,
+    })
+}
+
+// ============================================================================
+// Python Bindings
+// ============================================================================
+
+use pyo3::prelude::*;
+use pyo3::exceptions::PyValueError;
+use pyo3::types::PyDict;
+use numpy::{PyArray1, PyArray2, PyReadonlyArray1, PyReadonlyArray2};
+
+#[pyfunction]
+#[pyo3(signature = (covariance, n_components=1, lambda=0.1, max_iter=1000, tol=1e-6))]
+pub fn sparse_pca_py(
+    py: Python,
+    covariance: PyReadonlyArray2<f64>,
+    n_components: usize,
+    lambda: f64,
+    max_iter: usize,
+    tol: f64,
+) -> PyResult<PyObject> {
+    let cov = covariance.as_array().to_owned();
+    
+    let result = sparse_pca(&cov, n_components, lambda, max_iter, tol)
+        .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+    
+    let dict = PyDict::new_bound(py);
+    dict.set_item("weights", PyArray2::from_owned_array_bound(py, result.weights))?;
+    dict.set_item("variance_explained", PyArray1::from_owned_array_bound(py, result.variance_explained))?;
+    dict.set_item("sparsity", PyArray1::from_owned_array_bound(py, result.sparsity))?;
+    dict.set_item("iterations", result.iterations.to_vec())?;
+    dict.set_item("converged", result.converged)?;
+    
+    Ok(dict.into())
+}
+
+#[pyfunction]
+#[pyo3(signature = (matrix, lambda=0.1, mu=1.0, max_iter=500, tol=1e-5))]
+pub fn box_tao_decomposition_py(
+    py: Python,
+    matrix: PyReadonlyArray2<f64>,
+    lambda: f64,
+    mu: f64,
+    max_iter: usize,
+    tol: f64,
+) -> PyResult<PyObject> {
+    let mat = matrix.as_array().to_owned();
+    
+    let result = box_tao_decomposition(&mat, lambda, mu, max_iter, tol)
+        .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+    
+    let dict = PyDict::new_bound(py);
+    dict.set_item("low_rank", PyArray2::from_owned_array_bound(py, result.low_rank))?;
+    dict.set_item("sparse", PyArray2::from_owned_array_bound(py, result.sparse))?;
+    dict.set_item("noise", PyArray2::from_owned_array_bound(py, result.noise))?;
+    dict.set_item("rank", result.rank)?;
+    dict.set_item("sparsity", result.sparsity)?;
+    dict.set_item("iterations", result.iterations)?;
+    dict.set_item("converged", result.converged)?;
+    dict.set_item("objective_values", result.objective_values)?;
+    
+    Ok(dict.into())
+}
+
+#[pyfunction]
+#[pyo3(signature = (x, y, lambda_l1=0.1, lambda_l2=0.1, max_iter=1000, tol=1e-6))]
+pub fn elastic_net_py(
+    py: Python,
+    x: PyReadonlyArray2<f64>,
+    y: PyReadonlyArray1<f64>,
+    lambda_l1: f64,
+    lambda_l2: f64,
+    max_iter: usize,
+    tol: f64,
+) -> PyResult<PyObject> {
+    let x_arr = x.as_array().to_owned();
+    let y_arr = y.as_array().to_owned();
+    
+    let result = elastic_net(&x_arr, &y_arr, lambda_l1, lambda_l2, max_iter, tol)
+        .map_err(|e| PyValueError::new_err(format!("{}", e)))?;
+    
+    let dict = PyDict::new_bound(py);
+    dict.set_item("weights", PyArray1::from_owned_array_bound(py, result.weights))?;
+    dict.set_item("intercept", result.intercept)?;
+    dict.set_item("sparsity", result.sparsity)?;
+    dict.set_item("iterations", result.iterations)?;
+    dict.set_item("converged", result.converged)?;
+    dict.set_item("objective_value", result.objective_value)?;
+    
+    Ok(dict.into())
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use approx::assert_relative_eq;
+    
+    #[test]
+    fn test_soft_threshold() {
+        assert_relative_eq!(soft_threshold(3.0, 1.0), 2.0);
+        assert_relative_eq!(soft_threshold(-3.0, 1.0), -2.0);
+        assert_relative_eq!(soft_threshold(0.5, 1.0), 0.0);
+        assert_relative_eq!(soft_threshold(-0.5, 1.0), 0.0);
+    }
+    
+    #[test]
+    fn test_sparse_pca_simple() {
+        // Simple 3x3 covariance matrix
+        let cov = Array2::from_shape_vec(
+            (3, 3),
+            vec![
+                4.0, 2.0, 1.0,
+                2.0, 3.0, 1.0,
+                1.0, 1.0, 2.0,
+            ]
+        ).unwrap();
+        
+        let result = sparse_pca(&cov, 1, 0.1, 100, 1e-6).unwrap();
+        
+        assert_eq!(result.weights.nrows(), 1);
+        assert_eq!(result.weights.ncols(), 3);
+        assert!(result.variance_explained[0] > 0.0);
+        assert!(result.converged[0]);
+    }
+    
+    #[test]
+    fn test_elastic_net_simple() {
+        // Simple linear problem: y = 2*x1 + 3*x2
+        let x = Array2::from_shape_vec(
+            (5, 2),
+            vec![
+                1.0, 1.0,
+                2.0, 1.0,
+                3.0, 2.0,
+                4.0, 3.0,
+                5.0, 4.0,
+            ]
+        ).unwrap();
+        
+        let y = Array1::from_vec(vec![5.0, 7.0, 12.0, 17.0, 22.0]);
+        
+        let result = elastic_net(&x, &y, 0.01, 0.01, 1000, 1e-6).unwrap();
+        
+        assert!(result.converged);
+        // Coefficients should be approximately [2, 3]
+        assert!((result.weights[0] - 2.0).abs() < 0.5);
+        assert!((result.weights[1] - 3.0).abs() < 0.5);
+    }
+}