Function compute_lambda_max 

pub fn compute_lambda_max(
    x: &Matrix,
    y: &[f64],
    alpha: f64,
    penalty_factor: Option<&[f64]>,
    intercept_col: Option<usize>,
) -> f64

Computes lambda_max: the smallest lambda for which all penalized coefficients are zero.

Formula:

lambda_max = max(|X_j^T y|/vp(j)) / max(alpha, 1e-3)

where g(j) = |X_j^T y| is the absolute correlation.

Key points:

• y is assumed to be STANDARDIZED to unit norm (||y|| = 1)
• X columns are assumed to be STANDARDIZED to unit norm
• The 1e-3 minimum prevents division issues for pure ridge

Arguments

• x - Standardized design matrix (n × p), first column is intercept if present
• y - Standardized response vector (||y|| = 1)
• alpha - Elastic net mixing parameter
• penalty_factor - Per-feature penalty factors (optional, defaults to all 1.0)
• intercept_col - Index of intercept column (typically 0, or None if no intercept)

Returns

The lambda_max value used by glmnet for lambda path construction.