diff --git a/crates/ml-alpha/cuda/mamba2_alpha_kernel.cu b/crates/ml-alpha/cuda/mamba2_alpha_kernel.cu index 45af591a2..37e80a2a5 100644 --- a/crates/ml-alpha/cuda/mamba2_alpha_kernel.cu +++ b/crates/ml-alpha/cuda/mamba2_alpha_kernel.cu @@ -204,6 +204,57 @@ extern "C" __global__ void mamba2_alpha_reduce_d_proj( } +/* --------------------------------------------------------------------- + * AdamW step — bias-corrected decoupled-weight-decay update for one + * parameter tensor. Identical to ml-core's adamw_update kernel, included + * here so ml-alpha's cubin is self-contained (no cross-crate cubin loads). + * + * grad_scale: 1.0 disables gradient clipping; <1.0 applied as a scale. + * t: 1-indexed training step (for bias correction). + * + * Grid: ceil(n / 256), Block: 256. + * --------------------------------------------------------------------- */ +extern "C" __global__ void mamba2_alpha_adamw_step( + float* __restrict__ param, + const float* __restrict__ grad, + float* __restrict__ m, + float* __restrict__ v, + float lr, + float beta1, + float beta2, + float epsilon, + float weight_decay, + float grad_scale, + int t, + int n +) { + int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i >= n) return; + + float g = grad[i] * grad_scale; + float p = param[i]; + + /* Decoupled weight decay. */ + p = p * (1.0f - lr * weight_decay); + + /* Moment updates. */ + float mi = beta1 * m[i] + (1.0f - beta1) * g; + float vi = beta2 * v[i] + (1.0f - beta2) * g * g; + m[i] = mi; + v[i] = vi; + + /* Bias correction. */ + float bc1 = 1.0f - powf(beta1, (float)t); + float bc2 = 1.0f - powf(beta2, (float)t); + float m_hat = mi / bc1; + float v_hat = vi / bc2; + + /* Parameter update. */ + p = p - lr * m_hat / (sqrtf(v_hat) + epsilon); + param[i] = p; +} + + /* --------------------------------------------------------------------- * Reduction kernel: sum d_w_c_per_sample[N, sh2, state_d] across N → * d_w_c[sh2, state_d]. diff --git a/crates/ml-alpha/src/mamba2_block.rs b/crates/ml-alpha/src/mamba2_block.rs index c966bbbc1..26ebef9a1 100644 --- a/crates/ml-alpha/src/mamba2_block.rs +++ b/crates/ml-alpha/src/mamba2_block.rs @@ -168,6 +168,9 @@ pub struct Mamba2Block { /// Reduces `d_w_c_per_sample[N, sh2, state_d]` across N. pub kernel_reduce_d_w_c: CudaFunction, + /// AdamW step kernel — applied per parameter tensor. + pub kernel_adamw: CudaFunction, + /// cuBLAS handle for the four linear projections (input / A / B / output). pub cublas: CudaBlas, } @@ -198,6 +201,9 @@ impl Mamba2Block { let kernel_reduce_d_w_c = module .load_function("mamba2_alpha_reduce_d_w_c") .map_err(|e| anyhow!("Mamba2Block: d_w_c reduction kernel resolve: {e}"))?; + let kernel_adamw = module + .load_function("mamba2_alpha_adamw_step") + .map_err(|e| anyhow!("Mamba2Block: AdamW kernel resolve: {e}"))?; // cuBLAS handle for the four GEMM projections. let cublas = CudaBlas::new(Arc::clone(&stream)) @@ -241,6 +247,7 @@ impl Mamba2Block { kernel_bwd, kernel_reduce_d_proj, kernel_reduce_d_w_c, + kernel_adamw, cublas, }) } @@ -556,6 +563,210 @@ impl Mamba2Block { } } +// ===================================================================== +// AdamW optimizer for Mamba2Block +// ===================================================================== + +/// AdamW hyperparameters. +#[derive(Debug, Clone)] +pub struct Mamba2AdamWConfig { + pub lr: f32, + pub beta1: f32, + pub beta2: f32, + pub epsilon: f32, + pub weight_decay: f32, + pub grad_clip_max_norm: Option, +} + +impl Default for Mamba2AdamWConfig { + fn default() -> Self { + Self { + lr: 1e-3, + beta1: 0.9, + beta2: 0.999, + epsilon: 1e-8, + weight_decay: 1e-2, + grad_clip_max_norm: Some(1.0), + } + } +} + +/// Per-parameter Adam moment buffers (m, v) on GPU. Allocated once at +/// optimizer construction; reused across all training steps. +struct AdamState { + m: CudaSlice, + v: CudaSlice, +} + +/// GPU-native AdamW optimizer specialized for Mamba2Block's nine parameter +/// tensors. Holds the (m, v) moment state for each parameter; `step` applies +/// the AdamW update kernel to every weight + bias in place. +pub struct Mamba2AdamW { + pub config: Mamba2AdamWConfig, + pub step_count: i32, + stream: Arc, + kernel: CudaFunction, + + // (m, v) state for each of the 9 Mamba2 parameter tensors. + s_w_in: AdamState, + s_b_in: AdamState, + s_w_a: AdamState, + s_b_a: AdamState, + s_w_b: AdamState, + s_b_b: AdamState, + s_w_c: AdamState, + s_w_out: AdamState, + s_b_out: AdamState, +} + +impl Mamba2AdamW { + /// Allocate optimizer state matching the shapes of `block`'s parameters. + /// State buffers are zero-initialised; first AdamW step is unbiased. + pub fn new(block: &Mamba2Block, config: Mamba2AdamWConfig) -> Result { + let stream = Arc::clone(&block.stream); + let kernel = block.kernel_adamw.clone(); + + // Build all AdamStates first (immutable borrow of `stream` during + // allocation), then move `stream` into Self. + let alloc = |n: usize, name: &str| -> Result { + let m = stream.alloc_zeros::(n) + .map_err(|e| anyhow!("alloc Adam m for {name}: {e}"))?; + let v = stream.alloc_zeros::(n) + .map_err(|e| anyhow!("alloc Adam v for {name}: {e}"))?; + Ok(AdamState { m, v }) + }; + let s_w_in = alloc(block.w_in.weight.len(), "w_in.weight")?; + let s_b_in = alloc(block.w_in.bias.len(), "w_in.bias")?; + let s_w_a = alloc(block.w_a.weight.len(), "w_a.weight")?; + let s_b_a = alloc(block.w_a.bias.len(), "w_a.bias")?; + let s_w_b = alloc(block.w_b.weight.len(), "w_b.weight")?; + let s_b_b = alloc(block.w_b.bias.len(), "w_b.bias")?; + let s_w_c = alloc(block.w_c.len(), "w_c")?; + let s_w_out = alloc(block.w_out.weight.len(), "w_out.weight")?; + let s_b_out = alloc(block.w_out.bias.len(), "w_out.bias")?; + + Ok(Self { + stream, + kernel, + step_count: 0, + s_w_in, s_b_in, s_w_a, s_b_a, s_w_b, s_b_b, s_w_c, s_w_out, s_b_out, + config, + }) + } + + /// One AdamW step. Applies grad_clip_max_norm if configured (computed + /// from the L2 norm of the concatenated parameter gradient vector). + pub fn step(&mut self, block: &mut Mamba2Block, grads: &Mamba2BackwardGrads) -> Result<()> { + self.step_count += 1; + let t = self.step_count; + + // Compute grad-clip scale on host. The L2 norm requires a sum-of- + // squares reduction across all 9 tensors; cheapest is a host-side + // reduction over the 9 individual norms (since each is small). + let grad_scale = if let Some(max_norm) = self.config.grad_clip_max_norm { + // Each tensor's L2 squared norm = sum of (g_i * g_i). For Phase + // 1d.1 sizes the total parameter count is ~10k, so the host + // dtoh + sum is ~tens of microseconds. Re-evaluate at scale if + // it shows up in the training profile. + let mut total_sq = 0.0_f32; + for slice in [ + grads.dw_in.cuda_data(), grads.db_in.cuda_data(), + grads.dw_a.cuda_data(), grads.db_a.cuda_data(), + grads.dw_b.cuda_data(), grads.db_b.cuda_data(), + &grads.dw_c, + grads.dw_out.cuda_data(),grads.db_out.cuda_data(), + ] { + let mut host = vec![0.0_f32; slice.len()]; + self.stream.memcpy_dtoh(slice, &mut host) + .map_err(|e| anyhow!("grad-norm dtoh: {e}"))?; + for g in &host { + total_sq += g * g; + } + } + let norm = total_sq.sqrt(); + if norm > max_norm { max_norm / norm } else { 1.0 } + } else { + 1.0 + }; + + // Borrow the immutable bits as locals so we can mutably borrow + // self.s_* fields independently. + let stream = &self.stream; + let kernel = &self.kernel; + let cfg = &self.config; + + // Launch the AdamW kernel for each parameter tensor. + adamw_apply(stream, kernel, cfg, block.w_in.weight.len(), + &mut block.w_in.weight, grads.dw_in.cuda_data(), &mut self.s_w_in, t, grad_scale)?; + adamw_apply(stream, kernel, cfg, block.w_in.bias.len(), + &mut block.w_in.bias, grads.db_in.cuda_data(), &mut self.s_b_in, t, grad_scale)?; + adamw_apply(stream, kernel, cfg, block.w_a.weight.len(), + &mut block.w_a.weight, grads.dw_a.cuda_data(), &mut self.s_w_a, t, grad_scale)?; + adamw_apply(stream, kernel, cfg, block.w_a.bias.len(), + &mut block.w_a.bias, grads.db_a.cuda_data(), &mut self.s_b_a, t, grad_scale)?; + adamw_apply(stream, kernel, cfg, block.w_b.weight.len(), + &mut block.w_b.weight, grads.dw_b.cuda_data(), &mut self.s_w_b, t, grad_scale)?; + adamw_apply(stream, kernel, cfg, block.w_b.bias.len(), + &mut block.w_b.bias, grads.db_b.cuda_data(), &mut self.s_b_b, t, grad_scale)?; + adamw_apply(stream, kernel, cfg, block.w_c.len(), + &mut block.w_c, &grads.dw_c, &mut self.s_w_c, t, grad_scale)?; + adamw_apply(stream, kernel, cfg, block.w_out.weight.len(), + &mut block.w_out.weight, grads.dw_out.cuda_data(), &mut self.s_w_out, t, grad_scale)?; + adamw_apply(stream, kernel, cfg, block.w_out.bias.len(), + &mut block.w_out.bias, grads.db_out.cuda_data(), &mut self.s_b_out, t, grad_scale)?; + + Ok(()) + } + + pub fn set_learning_rate(&mut self, lr: f32) { + self.config.lr = lr; + } +} + +/// One AdamW kernel launch for a single parameter tensor. Free function +/// (rather than a method) so the caller can mutably borrow distinct fields +/// of `Mamba2AdamW` (the per-param Adam state) while sharing immutable +/// references to the stream + kernel + hyperparameters. +fn adamw_apply( + stream: &Arc, + kernel: &CudaFunction, + cfg: &Mamba2AdamWConfig, + n: usize, + param: &mut CudaSlice, + grad: &CudaSlice, + state: &mut AdamState, + t: i32, + grad_scale: f32, +) -> Result<()> { + let block_threads: u32 = 256; + let blocks: u32 = ((n + block_threads as usize - 1) / block_threads as usize) as u32; + let launch_cfg = LaunchConfig { + grid_dim: (blocks, 1, 1), + block_dim: (block_threads, 1, 1), + shared_mem_bytes: 0, + }; + let n_i32 = n as i32; + unsafe { + stream + .launch_builder(kernel) + .arg(param) + .arg(grad) + .arg(&mut state.m) + .arg(&mut state.v) + .arg(&cfg.lr) + .arg(&cfg.beta1) + .arg(&cfg.beta2) + .arg(&cfg.epsilon) + .arg(&cfg.weight_decay) + .arg(&grad_scale) + .arg(&t) + .arg(&n_i32) + .launch(launch_cfg) + .map_err(|e| anyhow!("mamba2_alpha_adamw_step launch (n={n}): {e}"))?; + } + Ok(()) +} + #[cfg(test)] mod tests { use super::*; @@ -584,6 +795,98 @@ mod tests { assert!(cfg.validate().is_err(), "seq_len > 32 must be rejected (backward x_hist limit)"); } + /// THE end-to-end correctness check: build a tiny Mamba2Block, run + /// 20 AdamW steps on a fixed input with a fixed binary label, and + /// assert mean BCE loss falls monotonically. If analytical backward + /// has a sign flip or wrong scale somewhere, this test fails by + /// observing the loss either flatten or rise. + #[test] + fn test_mamba2_block_training_loop_decreases_loss() { + let stream = match cuda_stream_or_skip() { + Some(s) => s, + None => return, + }; + let cfg = Mamba2BlockConfig { + in_dim: 4, hidden_dim: 8, state_dim: 4, seq_len: 8, + }; + let n_batch = 4; + let mut block = Mamba2Block::new(cfg.clone(), Arc::clone(&stream)).expect("init"); + let mut opt = Mamba2AdamW::new(&block, Mamba2AdamWConfig { + lr: 1e-2, // bigger lr — small model, 20-step convergence + ..Default::default() + }).expect("opt init"); + + // Fixed input + fixed binary labels (half +1, half 0). + let n = n_batch * cfg.seq_len * cfg.in_dim; + let h: Vec = (0..n).map(|i| ((i as f32) * 0.137).sin()).collect(); + let dev_in = stream.clone_htod(&h).expect("htod"); + let input = GpuTensor::new(dev_in, vec![n_batch, cfg.seq_len, cfg.in_dim]).expect("input"); + let labels: Vec = (0..n_batch).map(|i| if i < n_batch / 2 { 1.0 } else { 0.0 }).collect(); + + let bce_loss = |logits_host: &[f32], labels_host: &[f32]| -> f32 { + // Numerically-stable BCE-with-logits, mean over batch. + let mut s = 0.0_f32; + let eps = 1e-7_f32; + for (&z, &y) in logits_host.iter().zip(labels_host.iter()) { + let z = z.clamp(-50.0, 50.0); + let p = (1.0 / (1.0 + (-z).exp())).clamp(eps, 1.0 - eps); + s += -(y * p.ln() + (1.0 - y) * (1.0 - p).ln()); + } + s / logits_host.len() as f32 + }; + + let mut last_loss = f32::INFINITY; + let mut decrease_streak = 0; + let mut total_decreases = 0; + for step in 0..20 { + let (logit, cache) = block.forward_train(&input).expect("forward_train"); + let logit_host = logit.to_host(&stream).expect("logit dtoh"); + let loss = bce_loss(&logit_host, &labels); + + // d_logit = (sigmoid(z) - y) / N for mean-BCE-with-logits. + let d_logit_host: Vec = logit_host.iter().zip(labels.iter()) + .map(|(&z, &y)| { + let p = 1.0 / (1.0 + (-z.clamp(-50.0, 50.0)).exp()); + (p - y) / (n_batch as f32) + }) + .collect(); + let d_logit_dev = stream.clone_htod(&d_logit_host).expect("htod d_logit"); + let d_logit = GpuTensor::new(d_logit_dev, vec![n_batch, 1]).expect("d_logit"); + + let grads = block.backward(&cache, &d_logit).expect("backward"); + opt.step(&mut block, &grads).expect("opt step"); + + // Track convergence behaviour. + if loss < last_loss { + decrease_streak += 1; + total_decreases += 1; + } else { + decrease_streak = 0; + } + eprintln!("step={step:2} loss={loss:.5} streak={decrease_streak}"); + last_loss = loss; + } + + // Strict invariants: + // - Loss must have decreased a majority of the steps (allow some + // non-monotonicity from the AdamW momentum + small batch noise). + assert!( + total_decreases >= 15, + "expected ≥15 decreasing steps out of 20, got {total_decreases} — backward likely broken" + ); + // - Final loss must be meaningfully below the initial random-init loss + // (which starts near ln(2) ≈ 0.693 for balanced labels). Initial + // loss after random Xavier init varies; we just require we got + // below the chance baseline. + let (final_logit, _) = block.forward_train(&input).expect("final forward"); + let final_logit_host = final_logit.to_host(&stream).expect("final dtoh"); + let final_loss = bce_loss(&final_logit_host, &labels); + assert!( + final_loss < 0.65, + "final loss {final_loss} did not drop below 0.65 — training did not converge" + ); + } + #[test] fn test_mamba2_block_backward_returns_finite_grads_with_correct_shapes() { let stream = match cuda_stream_or_skip() {