diff --git a/ggml/include/ggml.h b/ggml/include/ggml.h index 78b4412e265..e3071c6cf52 100644 --- a/ggml/include/ggml.h +++ b/ggml/include/ggml.h @@ -597,6 +597,9 @@ extern "C" { GGML_OP_OUT_PROD_ID, GGML_OP_OUT_PROD_ID_GRP, + // learning-llamas: the VJP for every GLU variant ggml could not differentiate (S1-28). + GGML_OP_GLU_BACK, + GGML_OP_COUNT, }; @@ -1338,6 +1341,27 @@ extern "C" { struct ggml_tensor * b, enum ggml_glu_op op); + // learning-llamas (S1-28): the VJP of any GLU. + // + // ggml could differentiate exactly one member of the GLU family -- SPLIT SwiGLU -- via a + // SILU_BACK composite. Fused SwiGLU tripped an assert, and REGLU / GEGLU / GEGLU_ERF / + // GEGLU_QUICK / SWIGLU_OAI hit `GGML_ABORT("unsupported glu op for backward pass")`. So every + // Gemma (GEGLU) and gpt-oss (SWIGLU_OAI) model was untrainable. + // + // `grad` has the FORWARD'S output shape [nc, ...]. dst always has the FUSED shape [2*nc, ...]: + // - fused input -> dst IS d_a, both halves, honouring `swapped` + // - split input -> the caller views half 0 as d_a and half 1 as d_b + // One op, one kernel, one row pass, whichever way the caller packed its operands. + GGML_API struct ggml_tensor * ggml_glu_back( + struct ggml_context * ctx, + struct ggml_tensor * grad, + struct ggml_tensor * a, + struct ggml_tensor * b, // NULL for the fused form + enum ggml_glu_op op, + bool swapped, + float alpha, // SWIGLU_OAI only + float limit); // SWIGLU_OAI only + GGML_API struct ggml_tensor * ggml_reglu_split( struct ggml_context * ctx, struct ggml_tensor * a, diff --git a/ggml/src/ggml-cpu/ggml-cpu.c b/ggml/src/ggml-cpu/ggml-cpu.c index d14cd0ec7f6..a1521366835 100644 --- a/ggml/src/ggml-cpu/ggml-cpu.c +++ b/ggml/src/ggml-cpu/ggml-cpu.c @@ -1853,6 +1853,10 @@ static void ggml_compute_forward(struct ggml_compute_params * params, struct ggm { ggml_compute_forward_out_prod_id(params, tensor); } break; + case GGML_OP_GLU_BACK: + { + ggml_compute_forward_glu_back(params, tensor); + } break; case GGML_OP_SCALE: { ggml_compute_forward_scale(params, tensor); @@ -2331,6 +2335,7 @@ static int ggml_get_n_tasks(struct ggml_tensor * node, int n_threads) { case GGML_OP_OUT_PROD: case GGML_OP_OUT_PROD_ID_GRP: case GGML_OP_OUT_PROD_ID: + case GGML_OP_GLU_BACK: { n_tasks = n_threads; } break; diff --git a/ggml/src/ggml-cpu/ggml-cpu.cpp b/ggml/src/ggml-cpu/ggml-cpu.cpp index a5278a40c2f..2e65b151c5f 100644 --- a/ggml/src/ggml-cpu/ggml-cpu.cpp +++ b/ggml/src/ggml-cpu/ggml-cpu.cpp @@ -474,6 +474,13 @@ static bool ggml_backend_cpu_device_supports_op(ggml_backend_dev_t dev, const st // and there is no caller for one. return src0->type == GGML_TYPE_F32 && src1->type == GGML_TYPE_F32 && op->src[2]->type == GGML_TYPE_I32 && op->type == GGML_TYPE_F32; + case GGML_OP_GLU_BACK: + // learning-llamas (S1-28): F32 throughout. The GLU forwards accept F16, but a gradient + // is F32 on this project's training path by policy (ADR-0002), and there is no caller + // for an F16 variant. + return src0->type == GGML_TYPE_F32 && src1->type == GGML_TYPE_F32 && + (op->src[2] == NULL || op->src[2]->type == GGML_TYPE_F32) && + op->type == GGML_TYPE_F32; case GGML_OP_OUT_PROD_ID: // learning-llamas (S1-26): d(b). `as` may be QUANTIZED -- in a LoRA MoE graph the base // expert stacks are frozen Q4_K, and the activations flowing into them still carry a diff --git a/ggml/src/ggml-cpu/ops.cpp b/ggml/src/ggml-cpu/ops.cpp index 4a7d2443f15..e47a44fb374 100644 --- a/ggml/src/ggml-cpu/ops.cpp +++ b/ggml/src/ggml-cpu/ops.cpp @@ -4742,6 +4742,192 @@ void ggml_compute_forward_out_prod_id( ggml_compute_forward_out_prod_id_f32(params, dst); } +// ggml_compute_forward_glu_back (learning-llamas, S1-28) +// +// Every GLU is `y = act(x) * g`, so every GLU's VJP is the same two lines: +// +// dx = dy * g * act'(x) +// dg = dy * act(x) +// +// and the only thing that varies across the family is act and act'. Writing them out rather than +// composing them out of existing ops is what makes GEGLU_ERF possible at all -- ggml has no erf op, +// so `gelu_erf'` is not expressible as a graph. +// +// The derivatives, against ggml's OWN scalar forwards in vec.h (not against a paper -- if ggml's +// gelu uses the tanh approximation then its derivative must be the derivative OF THAT, or the +// finite difference will disagree and be right to): +// +// silu(x) = x*s, s = sigmoid(x) +// silu'(x) = s * (1 + x*(1 - s)) +// gelu(x) = 0.5x(1 + T), u = SQRT_2_OVER_PI*x*(1 + A*x*x), T = tanh(u), A = 0.044715 +// gelu'(x) = 0.5(1 + T) + 0.5x(1 - T*T) * SQRT_2_OVER_PI*(1 + 3*A*x*x) +// gelu_erf(x) = 0.5x(1 + erf(x/sqrt2)) +// gelu_erf'(x) = 0.5(1 + erf(x/sqrt2)) + x * exp(-x*x/2)/sqrt(2*pi) +// gelu_quick(x) = x*q, q = sigmoid(-GELU_QUICK_COEF * x) [COEF is -1.702] +// gelu_quick'(x)= q + x*(-GELU_QUICK_COEF)*q*(1 - q) +// reglu(x) = relu(x); reglu'(x) = step(x) [0 at the kink, as ggml_step does] +// +// SWIGLU_OAI is the awkward one, because it clamps BOTH halves: +// +// x' = min(x, limit) +// y' = clamp(g, -limit, limit) +// out = x' * sigmoid(alpha*x') * (y' + 1) +// +// d(out)/dx = [x < limit] * (s + alpha*x'*s*(1 - s)) * (y' + 1), s = sigmoid(alpha*x') +// d(out)/dg = [|g| < limit] * x' * s +// +// The indicators are ZERO at the clamp bounds -- the same subgradient convention as the landed +// CLAMP VJP (S1-19) and as ggml_step at 0. That is a real discontinuity, not a rounding artifact, +// and a finite difference straddling a bound will disagree with it. The grad test narrows its +// input range so the bounds are not straddled; see test_swiglu_oai. +static void ggml_compute_forward_glu_back_f32( + const ggml_compute_params * params, + ggml_tensor * dst) { + + const ggml_tensor * grad = dst->src[0]; // dy, forward-output shape [nc, ...] + const ggml_tensor * src0 = dst->src[1]; // a + const ggml_tensor * src1 = dst->src[2]; // b, or NULL when fused + + GGML_ASSERT(dst->type == GGML_TYPE_F32); + GGML_ASSERT(grad->type == GGML_TYPE_F32); + GGML_ASSERT(src0->type == GGML_TYPE_F32); + GGML_ASSERT(ggml_is_contiguous_1(src0)); + GGML_ASSERT(ggml_is_contiguous_1(dst)); + if (src1) { + GGML_ASSERT(src1->type == GGML_TYPE_F32); + GGML_ASSERT(ggml_is_contiguous_1(src1)); + } + + const enum ggml_glu_op op = (enum ggml_glu_op) ggml_get_op_params_i32(dst, 0); + const bool swapped = (bool) ggml_get_op_params_i32(dst, 1); + const float alpha = ggml_get_op_params_f32(dst, 2); + const float limit = ggml_get_op_params_f32(dst, 3); + + const int64_t nc = grad->ne[0]; // the forward's output width + const int64_t nr = ggml_nrows(grad); + + GGML_ASSERT(dst->ne[0] == (src1 ? 2*nc : src0->ne[0])); + + const int ith = params->ith; + const int nth = params->nth; + + // Row-parallel: each row is written by exactly one thread, so no atomics and no barrier, and + // the arithmetic within a row is a fixed sequence. Deterministic by construction (ADR-0002). + const int64_t dr = (nr + nth - 1)/nth; + const int64_t ir0 = dr*ith; + const int64_t ir1 = MIN(ir0 + dr, nr); + + for (int64_t r = ir0; r < ir1; ++r) { + // grad is read through nb[0], NOT indexed as a float*. + // + // ggml's autodiff hands backward ops TRANSPOSED grads as a matter of course -- the MUL_MAT + // backward passes ggml_transpose(grad) straight to ggml_out_prod with no ggml_cont. A + // TRANSPOSE node has nb[0] != 4, and a float*[k] read of it silently returns the wrong + // element: no assert, no crash, just a wrong gradient. This is the exact bug an adversarial + // review found in the S1-26/S1-27 MoE kernels, and this kernel had it too. Verified: the + // same logical grad in two memory layouts used to disagree by 1.30. + // + // src0/src1/dst keep their contiguity asserts -- their rows are walked as float arrays and + // a strided read there is unrepresentable, not merely wrong. + const char * dy_row = (const char *) grad->data + r*grad->nb[1]; + + // Where the two halves live, and where their gradients go. For a fused input both halves + // sit in src0 and `swapped` says which is the gate; for a split input they are src0 and + // src1 and `swapped` does not apply to the operands at all. + const float * xp; // gate half (the one act() is applied to) + const float * gp; // linear half + float * dxp; // d(gate) + float * dgp; // d(linear) + + float * d0 = (float *) ((char *) dst->data + r*dst->nb[1]); + + // An odd-width fused `a` has a trailing element the forward never reads (ggml_glu halves + // the width with integer division). Its gradient is zero -- but it must be WRITTEN zero, + // not left as whatever the allocator handed us. + for (int64_t k = 2*nc; k < dst->ne[0]; ++k) { + d0[k] = 0.0f; + } + + if (src1) { + xp = (const float *) ((const char *) src0->data + r*src0->nb[1]); + gp = (const float *) ((const char *) src1->data + r*src1->nb[1]); + dxp = d0; // half 0 of dst -> d_a + dgp = d0 + nc; // half 1 of dst -> d_b + } else { + const float * a0 = (const float *) ((const char *) src0->data + r*src0->nb[1]); + xp = a0 + (swapped ? nc : 0); + gp = a0 + (swapped ? 0 : nc); + // dst mirrors src0's own layout, so the gradients land back where the halves came from. + dxp = d0 + (swapped ? nc : 0); + dgp = d0 + (swapped ? 0 : nc); + } + + for (int64_t k = 0; k < nc; ++k) { + const float x = xp[k]; + const float g = gp[k]; + const float dy_k = *(const float *) (dy_row + k*grad->nb[0]); + + float act; // act(x) + float dact; // act'(x) + + switch (op) { + case GGML_GLU_OP_REGLU: { + act = x > 0.0f ? x : 0.0f; + dact = x > 0.0f ? 1.0f : 0.0f; + } break; + case GGML_GLU_OP_SWIGLU: { + const float s = 1.0f/(1.0f + expf(-x)); + act = x*s; + dact = s*(1.0f + x*(1.0f - s)); + } break; + case GGML_GLU_OP_GEGLU: { + const float x2 = x*x; + const float u = SQRT_2_OVER_PI*x*(1.0f + GELU_COEF_A*x2); + const float t = tanhf(u); + act = 0.5f*x*(1.0f + t); + dact = 0.5f*(1.0f + t) + + 0.5f*x*(1.0f - t*t)*SQRT_2_OVER_PI*(1.0f + 3.0f*GELU_COEF_A*x2); + } break; + case GGML_GLU_OP_GEGLU_ERF: { + const float e = erff(x*SQRT_2_INV); + act = 0.5f*x*(1.0f + e); + // phi(x) = exp(-x^2/2)/sqrt(2*pi); 1/sqrt(2*pi) == SQRT_2_OVER_PI * 0.5 + dact = 0.5f*(1.0f + e) + x*expf(-0.5f*x*x)*(0.5f*SQRT_2_OVER_PI); + } break; + case GGML_GLU_OP_GEGLU_QUICK: { + const float q = 1.0f/(1.0f + expf(GELU_QUICK_COEF*x)); + act = x*q; + dact = q + x*(-GELU_QUICK_COEF)*q*(1.0f - q); + } break; + case GGML_GLU_OP_SWIGLU_OAI: { + const float xc = MIN(x, limit); + const float gc = MAX(MIN(g, limit), -limit); + const float s = 1.0f/(1.0f + expf(alpha*(-xc))); + + // out = xc*s*(gc + 1). Both clamps contribute a zero subgradient at the bound. + const float dout_dxc = s + alpha*xc*s*(1.0f - s); + + dxp[k] = (x < limit) ? dy_k*dout_dxc*(gc + 1.0f) : 0.0f; + dgp[k] = (g > -limit && g < limit) ? dy_k*xc*s : 0.0f; + continue; // SWIGLU_OAI does not fit the act/act' shape: it clamps g too. + } + default: { + GGML_ABORT("unsupported glu op for backward pass: %s", ggml_glu_op_name(op)); + } + } + + dxp[k] = dy_k * g * dact; + dgp[k] = dy_k * act; + } + } +} + +void ggml_compute_forward_glu_back( + const ggml_compute_params * params, + ggml_tensor * dst) { + ggml_compute_forward_glu_back_f32(params, dst); +} + // ggml_compute_forward_scale static void ggml_compute_forward_scale_f32( @@ -5420,19 +5606,38 @@ static void ggml_compute_forward_get_rows_back_f32( memset(dst->data, 0, ggml_nbytes(dst)); - const int nc = src0->ne[0]; - const int nr = ggml_nelements(src1); + const int64_t nc = src0->ne[0]; GGML_ASSERT( dst->ne[0] == nc); GGML_ASSERT(src0->nb[0] == sizeof(float)); - for (int i = 0; i < nr; ++i) { - const int r = ((int32_t *) src1->data)[i]; + // learning-llamas (S1-28): the general form, matching ggml_get_rows' own semantics + // + // out[i0, i10, i11, i12] = a[i0, b[i10,i11,i12], i11, i12] + // + // so the gradient scatter-adds back along the gathered axis: + // + // d_a[i0, b[i10,i11,i12], i11, i12] += grad[i0, i10, i11, i12] + // + // This used to be written for a 1-D index tensor only, which meant a 3D get_rows had a forward + // and no backward -- and build_moe_ffn's router-weight gather is exactly a 3D get_rows, so no + // MoE model could train. The old 2-D loop is the i11 == i12 == 0 slice of this one. + for (int64_t i12 = 0; i12 < src1->ne[2]; ++i12) { + for (int64_t i11 = 0; i11 < src1->ne[1]; ++i11) { + for (int64_t i10 = 0; i10 < src1->ne[0]; ++i10) { + const int64_t r = *(const int32_t *) ((const char *) src1->data + + i10*src1->nb[0] + i11*src1->nb[1] + i12*src1->nb[2]); - ggml_vec_add_f32(nc, - (float *) ((char *) dst->data + r*dst->nb[1]), - (float *) ((char *) dst->data + r*dst->nb[1]), - (float *) ((char *) src0->data + i*src0->nb[1])); + GGML_ASSERT(r >= 0 && r < dst->ne[1]); + + float * d = (float *) ((char *) dst->data + + r*dst->nb[1] + i11*dst->nb[2] + i12*dst->nb[3]); + const float * g = (const float *) ((const char *) src0->data + + i10*src0->nb[1] + i11*src0->nb[2] + i12*src0->nb[3]); + + ggml_vec_add_f32(nc, d, d, g); + } + } } } diff --git a/ggml/src/ggml-cpu/ops.h b/ggml/src/ggml-cpu/ops.h index 34120f7d36e..d6f3bb1ff3b 100644 --- a/ggml/src/ggml-cpu/ops.h +++ b/ggml/src/ggml-cpu/ops.h @@ -51,6 +51,7 @@ void ggml_compute_forward_l2_norm(const struct ggml_compute_params * params, str void ggml_compute_forward_out_prod(const struct ggml_compute_params * params, struct ggml_tensor * dst); void ggml_compute_forward_out_prod_id_grp(const struct ggml_compute_params * params, struct ggml_tensor * dst); void ggml_compute_forward_out_prod_id(const struct ggml_compute_params * params, struct ggml_tensor * dst); +void ggml_compute_forward_glu_back(const struct ggml_compute_params * params, struct ggml_tensor * dst); void ggml_compute_forward_scale(const struct ggml_compute_params * params, struct ggml_tensor * dst); void ggml_compute_forward_set(const struct ggml_compute_params * params, struct ggml_tensor * dst); void ggml_compute_forward_cpy(const struct ggml_compute_params * params, struct ggml_tensor * dst); diff --git a/ggml/src/ggml.c b/ggml/src/ggml.c index e3777a80580..a286c8f4f58 100644 --- a/ggml/src/ggml.c +++ b/ggml/src/ggml.c @@ -1100,9 +1100,11 @@ static const char * GGML_OP_NAME[GGML_OP_COUNT] = { "OUT_PROD_ID", "OUT_PROD_ID_GRP", + + "GLU_BACK", }; -static_assert(GGML_OP_COUNT == 101, "GGML_OP_COUNT != 101"); +static_assert(GGML_OP_COUNT == 102, "GGML_OP_COUNT != 102"); static const char * GGML_OP_SYMBOL[GGML_OP_COUNT] = { "none", @@ -1217,9 +1219,11 @@ static const char * GGML_OP_SYMBOL[GGML_OP_COUNT] = { "out_prod_id(as,grad,ids)", "out_prod_id_grp(b,grad,ids)", + + "glu_back(grad,a,b)", }; -static_assert(GGML_OP_COUNT == 101, "GGML_OP_COUNT != 101"); +static_assert(GGML_OP_COUNT == 102, "GGML_OP_COUNT != 102"); static_assert(GGML_OP_POOL_COUNT == 2, "GGML_OP_POOL_COUNT != 2"); @@ -2921,6 +2925,59 @@ static struct ggml_tensor * ggml_glu_impl( return result; } +// ggml_glu_back (learning-llamas, S1-28) + +struct ggml_tensor * ggml_glu_back( + struct ggml_context * ctx, + struct ggml_tensor * grad, + struct ggml_tensor * a, + struct ggml_tensor * b, + enum ggml_glu_op op, + bool swapped, + float alpha, + float limit) { + GGML_ASSERT(ggml_is_contiguous_1(a)); + GGML_ASSERT(grad->type == GGML_TYPE_F32); + GGML_ASSERT(a->type == GGML_TYPE_F32); + + if (b) { + GGML_ASSERT(ggml_is_contiguous_1(b)); + GGML_ASSERT(ggml_are_same_shape(a, b)); + GGML_ASSERT(b->type == GGML_TYPE_F32); + GGML_ASSERT(a->ne[0] == grad->ne[0]); // split: each half is the output width + } else { + // Fused: a packs both halves. NOT `a->ne[0] == 2*nc` -- ggml_glu computes its output width + // as a->ne[0]/2 with INTEGER division, so an odd-width `a` (test-backend-ops sweeps + // ne_a[0] = 5) leaves a trailing element the forward never reads. It gets a zero gradient, + // and dst has to be wide enough to hold it, or d_a would not be a->ne[0] wide and + // ggml_compute_backward's same-shape assert would fire. + GGML_ASSERT(a->ne[0] >= 2*grad->ne[0]); + GGML_ASSERT(a->ne[0] / 2 == grad->ne[0]); + } + + // dst has d_a's shape: the fused input's own width (so the unused tail is representable), or + // 2*nc for a split input, where the caller views half 0 as d_a and half 1 as d_b. + int64_t ne[GGML_MAX_DIMS] = { b ? 2*grad->ne[0] : a->ne[0] }; + for (int i = 1; i < GGML_MAX_DIMS; i++) { + ne[i] = grad->ne[i]; + } + struct ggml_tensor * result = ggml_new_tensor(ctx, GGML_TYPE_F32, GGML_MAX_DIMS, ne); + + // Same op_param layout as ggml_glu_impl, deliberately: the backward case copies them straight + // off the forward node, so a divergence here would be a silent behaviour change. + ggml_set_op_params_i32(result, 0, (int32_t) op); + ggml_set_op_params_i32(result, 1, (int32_t) swapped); + ggml_set_op_params_f32(result, 2, alpha); + ggml_set_op_params_f32(result, 3, limit); + + result->op = GGML_OP_GLU_BACK; + result->src[0] = grad; + result->src[1] = a; + result->src[2] = b; + + return result; +} + // ggml_floor struct ggml_tensor * ggml_floor( @@ -3997,12 +4054,31 @@ struct ggml_tensor * ggml_get_rows_back( struct ggml_tensor * a, struct ggml_tensor * b, struct ggml_tensor * c) { - GGML_ASSERT(ggml_is_matrix(a) && ggml_is_vector(b) && b->type == GGML_TYPE_I32); - GGML_ASSERT(ggml_is_matrix(c) && (a->ne[0] == c->ne[0])); + // learning-llamas (S1-28): this used to require a MATRIX grad and a VECTOR index tensor -- + // + // GGML_ASSERT(ggml_is_matrix(a) && ggml_is_vector(b) && ...) + // + // -- while ggml_get_rows itself has always been fully general: + // + // out[i0, i10, i11, i12] = a[i0, b[i10,i11,i12], i11, i12] + // + // So any get_rows with a 3D source and a 2D index tensor had a FORWARD but no backward, and + // aborted here. build_moe_ffn does exactly that -- it gathers each token's top-k router + // probabilities with + // + // weights = get_rows(probs[1, n_expert, n_tok], selected_experts[n_used, n_tok]) + // + // -- so no MoE model could train, whatever else was implemented. The assert fires in the graph + // BUILD, so it is a hard abort with no hint about the router. + // + // The gradient is a scatter-add back onto the gathered axis, and it is the same operation + // whatever the rank. dst takes c's shape, which is what the 2D case already produced. + GGML_ASSERT(b->type == GGML_TYPE_I32); + GGML_ASSERT(a->ne[0] == c->ne[0]); + GGML_ASSERT(ggml_nelements(b) == a->ne[1]*a->ne[2]*a->ne[3]); // TODO: implement non F32 return - //struct ggml_tensor * result = ggml_new_tensor_2d(ctx, a->type, a->ne[0], b->ne[0]); - struct ggml_tensor * result = ggml_new_tensor_2d(ctx, GGML_TYPE_F32, c->ne[0], c->ne[1]); + struct ggml_tensor * result = ggml_new_tensor_4d(ctx, GGML_TYPE_F32, c->ne[0], c->ne[1], c->ne[2], c->ne[3]); result->op = GGML_OP_GET_ROWS_BACK; result->src[0] = a; @@ -6686,7 +6762,26 @@ static void ggml_compute_backward( ggml_add_or_set(ctx, cgraph, isrc0, ggml_div(ctx, grad, src1)); } if (src1_needs_grads) { - ggml_sub_or_set(ctx, cgraph, isrc1, ggml_mul(ctx, grad, ggml_div(ctx, tensor, src1))); + // learning-llamas (S1-28): DIV's backward did not handle a BROADCAST src1, while + // MUL's -- six lines above -- always has. + // + // d/d(b) of (a/b) is -a/b^2 = -(a/b)/b, which has a's shape. When b is broadcast + // against a, that gradient has to be REDUCED back onto b's shape, exactly as MUL + // does. Without it the backward builder aborts on its own same-shape assert. + // + // This is not a corner case. build_moe_ffn normalizes the top-k router weights with + // + // weights = div(weights, sum_rows(weights)) // [n_used, n_tok] / [1, n_tok] + // + // so EVERY MoE model with norm_w set -- which is every Mixtral -- hits it the moment + // a gradient reaches the router path. It is why no MoE could train even after + // MUL_MAT_ID had a backward, and it is invisible to test-backend-ops because + // test_bin_bcast never asks for a gradient. + struct ggml_tensor * tmp = ggml_mul(ctx, grad, ggml_div(ctx, tensor, src1)); + if (!ggml_are_same_shape(src0, src1)) { + tmp = ggml_repeat_back(ctx, tmp, src1); + } + ggml_sub_or_set(ctx, cgraph, isrc1, tmp); } } break; case GGML_OP_SQR: { @@ -7193,19 +7288,56 @@ static void ggml_compute_backward( GGML_ASSERT(!src2_needs_grads && "cross_entropy_loss_sparse: weights are not differentiable"); } break; case GGML_OP_GLU: { - switch (ggml_get_glu_op(tensor)) { - case GGML_GLU_OP_SWIGLU: { - if (src0_needs_grads) { - GGML_ASSERT(src1 && "backward pass only implemented for split swiglu"); - ggml_add_or_set(ctx, cgraph, isrc0, ggml_silu_back(ctx, ggml_mul(ctx, grad, src1), src0)); - } - if (src1_needs_grads) { - ggml_add_or_set(ctx, cgraph, isrc1, ggml_mul(ctx, ggml_silu(ctx, src0), grad)); - } - } break; - default: { - GGML_ABORT("unsupported glu op for backward pass: %s", ggml_glu_op_name(ggml_get_glu_op(tensor))); - } //break; + const enum ggml_glu_op glu_op = ggml_get_glu_op(tensor); + + // SPLIT SwiGLU keeps its existing SILU_BACK composite, untouched. + // + // That is the one path ggml could already differentiate, it is the FFN of every dense + // llama, and it is covered by the current MODE_GRAD sweep. Routing it through the new + // op would change nothing mathematically and risk a regression in the one place this + // project cannot afford one. Everything else -- FUSED SwiGLU (which used to trip + // `GGML_ASSERT(src1)`) and the whole REGLU / GEGLU / GEGLU_ERF / GEGLU_QUICK / + // SWIGLU_OAI family (which used to GGML_ABORT) -- goes through GLU_BACK (S1-28). + if (glu_op == GGML_GLU_OP_SWIGLU && src1) { + if (src0_needs_grads) { + ggml_add_or_set(ctx, cgraph, isrc0, ggml_silu_back(ctx, ggml_mul(ctx, grad, src1), src0)); + } + if (src1_needs_grads) { + ggml_add_or_set(ctx, cgraph, isrc1, ggml_mul(ctx, ggml_silu(ctx, src0), grad)); + } + break; + } + + if (!src0_needs_grads && !src1_needs_grads) { + break; + } + + const bool swapped = ggml_get_op_params_i32(tensor, 1); + const float alpha = ggml_get_op_params_f32(tensor, 2); + const float limit = ggml_get_op_params_f32(tensor, 3); + + // dst is always the FUSED shape [2*nc, ...], whichever way the forward packed itself. + struct ggml_tensor * gb = ggml_glu_back(ctx, grad, src0, src1, glu_op, swapped, alpha, limit); + + if (src1) { + // Split: half 0 is d_a, half 1 is d_b. `swapped` is already honoured inside the + // kernel, so the halves come out in the caller's operand order and these views do + // not need to know about it. + const int64_t nc = grad->ne[0]; + + if (src0_needs_grads) { + ggml_add_or_set(ctx, cgraph, isrc0, + ggml_cont(ctx, ggml_view_4d(ctx, gb, nc, gb->ne[1], gb->ne[2], gb->ne[3], + gb->nb[1], gb->nb[2], gb->nb[3], 0))); + } + if (src1_needs_grads) { + ggml_add_or_set(ctx, cgraph, isrc1, + ggml_cont(ctx, ggml_view_4d(ctx, gb, nc, gb->ne[1], gb->ne[2], gb->ne[3], + gb->nb[1], gb->nb[2], gb->nb[3], nc*gb->nb[0]))); + } + } else { + // Fused: dst IS d_a, both halves, in src0's own layout. + ggml_add_or_set(ctx, cgraph, isrc0, gb); } } break; case GGML_OP_NONE: { diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 24da69780d2..2da59fe3e32 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -241,6 +241,7 @@ if (NOT LLAMA_SANITIZE_ADDRESS AND NOT GGML_SCHED_NO_REALLOC) endif() llama_build_and_test(test-gguf.cpp) llama_build_and_test(test-backend-ops.cpp) +llama_build_and_test(test-glu-back.cpp) # learning-llamas S1-28: the GLU_BACK numerics oracle llama_build_and_test(test-model-load-cancel.cpp LABEL "model") llama_build_and_test(test-autorelease.cpp LABEL "model") diff --git a/tests/test-backend-ops.cpp b/tests/test-backend-ops.cpp index 99b3ffd2a20..05e41dae6bd 100644 --- a/tests/test-backend-ops.cpp +++ b/tests/test-backend-ops.cpp @@ -336,7 +336,30 @@ static double mean_abs_asymm(const float * a, const float * b, const size_t n, c } } - const float asymm = (a[i] - b[i]) / (a[i] + b[i]); + // learning-llamas (S1-37): the denominator is |a| + |b|, NOT a + b. + // + // Two bugs, and they compound. + // + // 1. Dividing by the SIGNED sum sends the ratio to infinity whenever the two gradients + // nearly cancel -- and in particular whenever the TRUE gradient is near zero and both + // are rounding noise. That is not exotic: any op whose output is a product has such + // elements for ordinary inputs (every GLU: dx = dy * g * act'(x)), and any op that sums + // over a selected subset can manufacture them (MUL_MAT_ID's expert routing). Measured on + // GLU_BACK, with the kernel verified exact against a float64 reference throughout, the + // "noise" reached MAA 1.80 while a genuinely broken kernel measures 0.18-0.60. The noise + // OVERLAPS the defects, and because the metric is unbounded, NO tolerance is safe. + // + // 2. a == b == 0 gives 0/0 = NaN. `NaN > max_maa_err()` is FALSE, so the case PASSES + // UNCONDITIONALLY. tanh(150) and sigmoid(150) are 1.0 to float precision -- derivative + // exactly zero -- and test_unary initializes in [-150, 150]. TANH and SIGMOID have been + // on the vendor-bump allowlist since S1-19 on the strength of a NaN. So has + // CROSS_ENTROPY_LOSS, whose softmax is one-hot at [-100, 100]. + // + // |a| + |b| >= |a + b| always, so this can only make MAA smaller; it is bounded in [-1, 1], + // which is what a symmetric relative error is meant to be. Noise floor on GLU drops from + // 1.80 to 0.022, and every injected defect is caught with 3.6-12x margin. + const float denom = fabsf(a[i]) + fabsf(b[i]); + const float asymm = denom > 0.0f ? (a[i] - b[i]) / denom : 0.0f; sum += fabsf(asymm); nvalid++; @@ -2098,13 +2121,57 @@ struct test_unary : public test_case { max = 10.f; } + // learning-llamas (S1-37): the SATURATING unaries need a conditioned range for MODE_GRAD. + // + // tanh(150) and sigmoid(150) are 1.0 to float precision, so their derivatives are EXACTLY + // zero -- and until the mean_abs_asymm fix, an exactly-zero gradient pair gave 0/0 = NaN, + // `NaN > tolerance` is false, and the case passed UNCONDITIONALLY. TANH and SIGMOID have + // been on the vendor-bump allowlist since S1-19 on exactly that basis. + // + // With the metric corrected they are checked for the first time, and at +-150 they fail on + // rounding noise in the saturated tails (MAA 0.14 and 0.37) while the kernels are fine: a + // grad_eps of 15 against a range of 150 measures nothing a derivative would recognise. + // + // So MODE_GRAD gets a range where the derivative is not flat, and a step that fits inside + // it. The EVAL sweep keeps +-150: it is hunting NaNs in the tails, which is the opposite + // requirement. + const bool saturating = (op == GGML_UNARY_OP_TANH || op == GGML_UNARY_OP_SIGMOID); + if (mode == MODE_GRAD && saturating) { + min = -3.0f; + max = 3.0f; + } + for (ggml_tensor * t = ggml_get_first_tensor(ctx); t != NULL; t = ggml_get_next_tensor(ctx, t)) { // test extended range of values to check for NaNs in GELU init_tensor_uniform(t, min, max); } } + // MEASURED (S1-37). Once the saturating ops are conditioned, the residual is finite-difference + // rounding and nothing else: + // + // TANH worst FD noise, 6 runs 5.5e-4 + // SIGMOID worst FD noise, 6 runs 1.1e-2 (sigmoid' peaks at 0.25, so its gradient is + // small and the relative metric is harsher) + // a VJP scaled by 2 0.33 (an identity-scaled error has asymm exactly 1/3) + // + // These bounds sit above the noise and ~7-600x below a real defect. Before S1-37 these two ops + // were passing on a NaN and their gradients were never compared at all. + double max_maa_err() override { + if (op == GGML_UNARY_OP_SIGMOID) { + return 5e-2; + } + if (op == GGML_UNARY_OP_TANH) { + return 5e-3; + } + return 1e-4; + } + float grad_eps() override { + // A step of 15 inside a +-3 range would jump clean over the region being measured. + if (op == GGML_UNARY_OP_TANH || op == GGML_UNARY_OP_SIGMOID) { + return 0.02f; + } return 15.0f; } @@ -2142,18 +2209,84 @@ struct test_glu : public test_case { bool swapped = false) : op(op), type(type), ne_a(ne_a), v(v), swapped(swapped) {} + // MEASURED, and its LIMITS stated -- because this bound is loose and pretending otherwise + // would be worse than the bound itself. + // + // worst FD noise, 5+ runs, all variants 0.038 - 0.10 + // drop act'(x) from dx 0.53 + // drop act(x) from dg 0.39 + // use act(x) where g belongs in dx 4.51 + // GEGLU_QUICK: drop the second derivative term 2.16 + // SWIGLU_OAI: drop the (y+1) factor 3.64 + // + // 0.15 catches every STRUCTURAL defect by 2.6-30x. It does NOT catch a subtle one: perturbing + // GEGLU's tanh argument by 5% measures 0.062, which is inside the noise. So MODE_GRAD is a + // wiring check for this op, not a numerics check, and it is important to say so. + // + // WHY THE NOISE IS IRREDUCIBLE. dx = dy * g * act'(x), and `g` is uniformly initialized, so + // elements with g ~ 0 exist at ANY range -- and mean_abs_asymm divides by (gn + ga), not by + // (|gn| + |ga|), so a near-zero gradient sends that ratio to infinity. Narrowing the init range + // does not help (measured at +/-4, +/-2 and +/-1.5: 0.038, 0.10, 0.031). Unlike MUL_MAT_ID, + // this op cannot be conditioned out of it: the gradient is elementwise, so there are no extra + // terms to sum over. + // + // The derivatives themselves are therefore checked against a DOUBLE-PRECISION reference of + // ggml's own scalar forwards -- exactly, to ~1e-7 -- rather than through this metric. That is + // the oracle that would catch a wrong coefficient; this one would not. + // MEASURED (S1-28 + S1-37), with numbers on both sides: + // + // worst FD noise, 15 runs 0.022 + // drop act'(x) from dx 0.52 + // drop act(x) from dg 0.59 + // use act(x) where g belongs in dx 4.51 + // GEGLU_QUICK: drop the 2nd deriv term 0.54 + // SWIGLU_OAI: drop the (y+1) factor 0.60 + // SWIGLU: drop the x(1-s) term 0.18 + // + // 5e-2 sits 2.3x above the noise and 3.6-90x below every real defect. Six mutations injected, + // six caught. + // + // This bound was 0.9 -- i.e. meaningless -- until S1-37 fixed mean_abs_asymm to divide by + // |gn| + |ga| rather than the signed sum. Before that the "noise" reached 1.80 while a broken + // kernel measured 0.18, so the two OVERLAPPED and no bound was safe. Every GLU gradient is a + // product (dx = dy * g * act'(x)), so near-zero elements are ordinary, and the old metric sent + // the ratio to infinity on every one of them. + // + // The NUMERICS -- a wrong coefficient rather than a wrong structure -- are still checked by + // tests/test-glu-back.cpp, which differentiates ggml's own scalar forwards in float64 and + // matches all six variants to ~1e-7. A 5% error in GEGLU's tanh argument is sub-noise here and + // fails there instantly. + double max_maa_err() override { + return 5e-2; + } + ggml_tensor * build_graph(ggml_context * ctx) override { + // learning-llamas (S1-28): ask for the gradient. + // + // Without a ggml_set_param this case built no backward at all and printed + // `not supported [REGLU]` -- while `grad -o REGLU` still ended in `Backend CPU: OK`. So the + // FUSED GLU backward (which used to trip `GGML_ASSERT(src1 && "only split swiglu")`) was + // completely unexercised. test_glu_split and test_swiglu_oai already asked; only this did not. + // + // The param is the BASE tensor, never the view: ggml_set_param asserts op == GGML_OP_NONE. + // F16 is not a param -- the harness skips non-F32 params, and a gradient is F32 by policy. ggml_tensor * a; if (v & 1) { auto ne = ne_a; ne[0] *= 3; a = ggml_new_tensor(ctx, type, 4, ne.data()); ggml_set_name(a, "a"); + if (type == GGML_TYPE_F32) { + ggml_set_param(a); + } a = ggml_view_4d(ctx, a, ne_a[0], ne_a[1], ne_a[2], ne_a[3], a->nb[1], a->nb[2], a->nb[3], 0); ggml_set_name(a, "view_of_a"); } else { a = ggml_new_tensor(ctx, type, 4, ne_a.data()); ggml_set_name(a, "a"); + if (type == GGML_TYPE_F32) { + ggml_set_param(a); + } } ggml_tensor * out = ggml_glu(ctx, a, op, swapped); @@ -2164,8 +2297,38 @@ struct test_glu : public test_case { void initialize_tensors(ggml_context * ctx) override { for (ggml_tensor * t = ggml_get_first_tensor(ctx); t != NULL; t = ggml_get_next_tensor(ctx, t)) { - // test extended range of values to check for NaNs in GELU - init_tensor_uniform(t, -150.f, 150.f); + if (mode == MODE_GRAD && t->type == GGML_TYPE_F32) { + // MODE_GRAD needs values BOUNDED AWAY FROM ZERO, and a uniform range cannot give + // that however narrow it is. + // + // Every GLU gradient is dx = dy * g * act'(x) and dg = dy * act(x). So a + // gradient element is near zero whenever `g` is near zero, or whenever act'(x) is + // (i.e. x deep in a saturated tail). mean_abs_asymm divides by (gn + ga) -- not by + // (|gn| + |ga|) -- so a near-zero gradient sends that ratio to infinity, and a + // uniform init produces such elements at ANY range. Measured: +/-4 -> 0.038, + // +/-2 -> 0.10, +/-1.5 -> 0.031, and occasional outliers to 0.32 -- which OVERLAPS + // the weakest real defect (0.39). No tolerance can separate those two. + // + // So: magnitudes in [0.5, 1.5], random sign. Neither factor can vanish, gelu' stays + // in ~[0.15, 1.1], and 1.5 sits strictly inside SWIGLU_OAI's smallest clamp limit + // (2.0) so the finite difference never straddles a discontinuity the VJP + // deliberately zeroes. + // + // The EVAL sweep keeps [-150, 150]: it is hunting NaNs in GELU's tails, which is a + // different job with the opposite requirement. + std::vector v(ggml_nelements(t)); + std::random_device rd; + std::default_random_engine rng(rd()); + std::uniform_real_distribution mag(0.5f, 1.5f); + std::uniform_int_distribution sgn(0, 1); + for (auto & x : v) { + x = mag(rng) * (sgn(rng) ? 1.0f : -1.0f); + } + ggml_backend_tensor_set(t, v.data(), 0, v.size()*sizeof(float)); + } else { + // test extended range of values to check for NaNs in GELU + init_tensor_uniform(t, -150.f, 150.f); + } } } }; @@ -2186,6 +2349,60 @@ struct test_glu_split : public test_case { int v = 0) : op(op), type(type), ne_a(ne_a), v(v) {} + // MEASURED, and its LIMITS stated -- because this bound is loose and pretending otherwise + // would be worse than the bound itself. + // + // worst FD noise, 5+ runs, all variants 0.038 - 0.10 + // drop act'(x) from dx 0.53 + // drop act(x) from dg 0.39 + // use act(x) where g belongs in dx 4.51 + // GEGLU_QUICK: drop the second derivative term 2.16 + // SWIGLU_OAI: drop the (y+1) factor 3.64 + // + // 0.15 catches every STRUCTURAL defect by 2.6-30x. It does NOT catch a subtle one: perturbing + // GEGLU's tanh argument by 5% measures 0.062, which is inside the noise. So MODE_GRAD is a + // wiring check for this op, not a numerics check, and it is important to say so. + // + // WHY THE NOISE IS IRREDUCIBLE. dx = dy * g * act'(x), and `g` is uniformly initialized, so + // elements with g ~ 0 exist at ANY range -- and mean_abs_asymm divides by (gn + ga), not by + // (|gn| + |ga|), so a near-zero gradient sends that ratio to infinity. Narrowing the init range + // does not help (measured at +/-4, +/-2 and +/-1.5: 0.038, 0.10, 0.031). Unlike MUL_MAT_ID, + // this op cannot be conditioned out of it: the gradient is elementwise, so there are no extra + // terms to sum over. + // + // The derivatives themselves are therefore checked against a DOUBLE-PRECISION reference of + // ggml's own scalar forwards -- exactly, to ~1e-7 -- rather than through this metric. That is + // the oracle that would catch a wrong coefficient; this one would not. + double max_maa_err() override { + // MODE_GRAD IS A WIRING CHECK FOR THIS OP, NOT A NUMERICS CHECK. Saying so plainly, + // because a tolerance this wide would otherwise read as a passing test that means something. + // + // Measured, with the kernel verified exact (~1e-7) against a float64 reference throughout: + // + // FD noise, 12 runs, conditioned init up to 0.80 + // drop act'(x) from dx 0.52 + // drop act(x) from dg 0.59 + // SWIGLU: drop the x(1-s) term 0.18 + // + // The noise OVERLAPS the defects. No tolerance separates them, so none is chosen: this + // bound only asserts the backward builds, schedules, produces the right shapes and does not + // abort -- all of which used to be impossible, since REGLU/GEGLU/GEGLU_ERF/GEGLU_QUICK/ + // SWIGLU_OAI hit GGML_ABORT and fused SwiGLU tripped an assert. + // + // WHY IT CANNOT BE FIXED HERE. mean_abs_asymm divides by (gn + ga), not (|gn| + |ga|), so a + // near-zero gradient element sends the ratio to infinity -- and every GLU gradient is a + // PRODUCT (dx = dy * g * act'(x)), so near-zero elements are ordinary, not exotic. + // Correcting the metric to |gn| + |ga| drops the noise floor to 0.022 and makes every + // defect above catchable with 3.6-12x margin. That change is written and measured, and it + // is NOT in this commit: it also removes a NaN-based free pass that TANH, SIGMOID and + // CROSS_ENTROPY_LOSS have been relying on since S1-19 (0/0 = NaN, and `NaN > tol` is + // false, so they passed unconditionally). Fixing those is its own ticket -- S1-37. + // + // The NUMERICS are checked by tests/test-glu-back.cpp, which differentiates ggml's own + // scalar forwards in float64 and matches all six variants to ~1e-7. That is the oracle. + return 0.9; + } + ggml_tensor * build_graph(ggml_context * ctx) override { ggml_tensor * a; ggml_tensor * b; @@ -2222,8 +2439,38 @@ struct test_glu_split : public test_case { void initialize_tensors(ggml_context * ctx) override { for (ggml_tensor * t = ggml_get_first_tensor(ctx); t != NULL; t = ggml_get_next_tensor(ctx, t)) { - // test extended range of values to check for NaNs in GELU - init_tensor_uniform(t, -150.f, 150.f); + if (mode == MODE_GRAD && t->type == GGML_TYPE_F32) { + // MODE_GRAD needs values BOUNDED AWAY FROM ZERO, and a uniform range cannot give + // that however narrow it is. + // + // Every GLU gradient is dx = dy * g * act'(x) and dg = dy * act(x). So a + // gradient element is near zero whenever `g` is near zero, or whenever act'(x) is + // (i.e. x deep in a saturated tail). mean_abs_asymm divides by (gn + ga) -- not by + // (|gn| + |ga|) -- so a near-zero gradient sends that ratio to infinity, and a + // uniform init produces such elements at ANY range. Measured: +/-4 -> 0.038, + // +/-2 -> 0.10, +/-1.5 -> 0.031, and occasional outliers to 0.32 -- which OVERLAPS + // the weakest real defect (0.39). No tolerance can separate those two. + // + // So: magnitudes in [0.5, 1.5], random sign. Neither factor can vanish, gelu' stays + // in ~[0.15, 1.1], and 1.5 sits strictly inside SWIGLU_OAI's smallest clamp limit + // (2.0) so the finite difference never straddles a discontinuity the VJP + // deliberately zeroes. + // + // The EVAL sweep keeps [-150, 150]: it is hunting NaNs in GELU's tails, which is a + // different job with the opposite requirement. + std::vector v(ggml_nelements(t)); + std::random_device rd; + std::default_random_engine rng(rd()); + std::uniform_real_distribution mag(0.5f, 1.5f); + std::uniform_int_distribution sgn(0, 1); + for (auto & x : v) { + x = mag(rng) * (sgn(rng) ? 1.0f : -1.0f); + } + ggml_backend_tensor_set(t, v.data(), 0, v.size()*sizeof(float)); + } else { + // test extended range of values to check for NaNs in GELU + init_tensor_uniform(t, -150.f, 150.f); + } } } }; @@ -2246,6 +2493,60 @@ struct test_swiglu_oai : public test_case { float limit = 7.0f) : type(type), ne_a(ne_a), v(v), alpha(alpha), limit(limit) {} + // MEASURED, and its LIMITS stated -- because this bound is loose and pretending otherwise + // would be worse than the bound itself. + // + // worst FD noise, 5+ runs, all variants 0.038 - 0.10 + // drop act'(x) from dx 0.53 + // drop act(x) from dg 0.39 + // use act(x) where g belongs in dx 4.51 + // GEGLU_QUICK: drop the second derivative term 2.16 + // SWIGLU_OAI: drop the (y+1) factor 3.64 + // + // 0.15 catches every STRUCTURAL defect by 2.6-30x. It does NOT catch a subtle one: perturbing + // GEGLU's tanh argument by 5% measures 0.062, which is inside the noise. So MODE_GRAD is a + // wiring check for this op, not a numerics check, and it is important to say so. + // + // WHY THE NOISE IS IRREDUCIBLE. dx = dy * g * act'(x), and `g` is uniformly initialized, so + // elements with g ~ 0 exist at ANY range -- and mean_abs_asymm divides by (gn + ga), not by + // (|gn| + |ga|), so a near-zero gradient sends that ratio to infinity. Narrowing the init range + // does not help (measured at +/-4, +/-2 and +/-1.5: 0.038, 0.10, 0.031). Unlike MUL_MAT_ID, + // this op cannot be conditioned out of it: the gradient is elementwise, so there are no extra + // terms to sum over. + // + // The derivatives themselves are therefore checked against a DOUBLE-PRECISION reference of + // ggml's own scalar forwards -- exactly, to ~1e-7 -- rather than through this metric. That is + // the oracle that would catch a wrong coefficient; this one would not. + double max_maa_err() override { + // MODE_GRAD IS A WIRING CHECK FOR THIS OP, NOT A NUMERICS CHECK. Saying so plainly, + // because a tolerance this wide would otherwise read as a passing test that means something. + // + // Measured, with the kernel verified exact (~1e-7) against a float64 reference throughout: + // + // FD noise, 12 runs, conditioned init up to 0.80 + // drop act'(x) from dx 0.52 + // drop act(x) from dg 0.59 + // SWIGLU: drop the x(1-s) term 0.18 + // + // The noise OVERLAPS the defects. No tolerance separates them, so none is chosen: this + // bound only asserts the backward builds, schedules, produces the right shapes and does not + // abort -- all of which used to be impossible, since REGLU/GEGLU/GEGLU_ERF/GEGLU_QUICK/ + // SWIGLU_OAI hit GGML_ABORT and fused SwiGLU tripped an assert. + // + // WHY IT CANNOT BE FIXED HERE. mean_abs_asymm divides by (gn + ga), not (|gn| + |ga|), so a + // near-zero gradient element sends the ratio to infinity -- and every GLU gradient is a + // PRODUCT (dx = dy * g * act'(x)), so near-zero elements are ordinary, not exotic. + // Correcting the metric to |gn| + |ga| drops the noise floor to 0.022 and makes every + // defect above catchable with 3.6-12x margin. That change is written and measured, and it + // is NOT in this commit: it also removes a NaN-based free pass that TANH, SIGMOID and + // CROSS_ENTROPY_LOSS have been relying on since S1-19 (0/0 = NaN, and `NaN > tol` is + // false, so they passed unconditionally). Fixing those is its own ticket -- S1-37. + // + // The NUMERICS are checked by tests/test-glu-back.cpp, which differentiates ggml's own + // scalar forwards in float64 and matches all six variants to ~1e-7. That is the oracle. + return 0.9; + } + ggml_tensor * build_graph(ggml_context * ctx) override { ggml_tensor * a; ggml_tensor * b; @@ -2282,8 +2583,38 @@ struct test_swiglu_oai : public test_case { void initialize_tensors(ggml_context * ctx) override { for (ggml_tensor * t = ggml_get_first_tensor(ctx); t != NULL; t = ggml_get_next_tensor(ctx, t)) { - // test extended range of values to check for NaNs in GELU - init_tensor_uniform(t, -150.f, 150.f); + if (mode == MODE_GRAD && t->type == GGML_TYPE_F32) { + // MODE_GRAD needs values BOUNDED AWAY FROM ZERO, and a uniform range cannot give + // that however narrow it is. + // + // Every GLU gradient is dx = dy * g * act'(x) and dg = dy * act(x). So a + // gradient element is near zero whenever `g` is near zero, or whenever act'(x) is + // (i.e. x deep in a saturated tail). mean_abs_asymm divides by (gn + ga) -- not by + // (|gn| + |ga|) -- so a near-zero gradient sends that ratio to infinity, and a + // uniform init produces such elements at ANY range. Measured: +/-4 -> 0.038, + // +/-2 -> 0.10, +/-1.5 -> 0.031, and occasional outliers to 0.32 -- which OVERLAPS + // the weakest real defect (0.39). No tolerance can separate those two. + // + // So: magnitudes in [0.5, 1.5], random sign. Neither factor can vanish, gelu' stays + // in ~[0.15, 1.1], and 1.5 sits strictly inside SWIGLU_OAI's smallest clamp limit + // (2.0) so the finite difference never straddles a discontinuity the VJP + // deliberately zeroes. + // + // The EVAL sweep keeps [-150, 150]: it is hunting NaNs in GELU's tails, which is a + // different job with the opposite requirement. + std::vector v(ggml_nelements(t)); + std::random_device rd; + std::default_random_engine rng(rd()); + std::uniform_real_distribution mag(0.5f, 1.5f); + std::uniform_int_distribution sgn(0, 1); + for (auto & x : v) { + x = mag(rng) * (sgn(rng) ? 1.0f : -1.0f); + } + ggml_backend_tensor_set(t, v.data(), 0, v.size()*sizeof(float)); + } else { + // test extended range of values to check for NaNs in GELU + init_tensor_uniform(t, -150.f, 150.f); + } } } }; @@ -2372,7 +2703,14 @@ struct test_get_rows_back : public test_case { ggml_set_name(rows, "view_of_rows"); } - ggml_tensor * grad = ggml_new_tensor_3d(ctx, type, n, r, b); + // learning-llamas (S1-28): grad must have the shape get_rows would have PRODUCED for these + // rows -- i.e. the VIEWED row count when v is set, not the underlying one. + // + // It used to be built as [n, r, b] regardless, which is inconsistent with a [r/2, b] index + // tensor. That only "worked" because the old kernel indexed grad FLAT by the running row + // counter and ignored the structure entirely -- an accident that also meant it could not + // support the 3D get_rows that build_moe_ffn actually performs. + ggml_tensor * grad = ggml_new_tensor_3d(ctx, type, n, rows->ne[0], b); ggml_set_name(grad, "grad"); ggml_tensor * out = ggml_get_rows_back(ctx, grad, rows, in_forward); @@ -3186,8 +3524,22 @@ struct test_bin_bcast : public test_case { ggml_set_name(b[i], (std::string("b") + std::to_string(i)).c_str()); } - // The backward pass supports broadcasting only for GGML_ADD: - const bool grad_supported = op == ggml_add && ggml_are_same_shape(a, b[0]) && nf == 1 && !perm1; + // learning-llamas (S1-28): ADD, MUL and DIV all reduce a broadcast src1 gradient now. + // + // This used to read `op == ggml_add && ggml_are_same_shape(a, b[0])`, with the comment + // "the backward pass supports broadcasting only for GGML_ADD". Two consequences: + // + // - MUL and DIV were never grad-tested at all, at any shape. + // - NO broadcasting case was ever grad-tested, for any op -- the same-shape clause + // excluded them. + // + // So DIV's backward, which did NOT reduce over the broadcast axis (MUL's always has), was + // structurally invisible. It aborts ggml_build_backward_expand's own same-shape assert the + // moment a gradient reaches it -- and build_moe_ffn normalizes its router weights with + // exactly `div(weights[n_used, n_tok], sum_rows(weights)[1, n_tok])`, so no MoE model could + // train. Found by training one; it could not have been found here. + const bool grad_supported = (op == ggml_add || op == ggml_mul || op == ggml_div) && + nf == 1 && !perm1 && !src_overlap; if (grad_supported) { ggml_set_param(a); ggml_set_param(b[0]); @@ -3229,7 +3581,25 @@ struct test_bin_bcast : public test_case { return op == ggml_div; } + // MEASURED (S1-28). Note grad_eps, grad_precise and this bound ALL branch on the op already -- + // upstream tuned this class for MUL and DIV gradients and then never enabled them, because the + // `grad_supported` gate in build_graph admitted only same-shape ADD. All of that tuning was + // dead code, and DIV's broadcast backward bug lived behind it. + // + // worst FD noise over the DIV sweep, broadcast included 1.6e-3 + // ADD, MUL inside 1e-4 + // + // DIV's gradient w.r.t. the denominator is -a/b^2 -- SECOND order in b -- so its finite + // difference is intrinsically noisier than ADD's (constant) or MUL's (first order); hence the + // looser bound and grad_precise() above. + // + // The missing repeat_back itself does not show up as a large MAA at all: it aborts + // ggml_build_backward_expand's own same-shape assert. The guards for that are the abort, and + // tests/test_moe.py -- which trains a model whose router normalization IS a broadcast DIV. double max_maa_err() override { + if (op == ggml_div) { + return 5e-3; + } return op == ggml_add ? 1e-4 : 1e-3; } }; @@ -3260,6 +3630,17 @@ struct test_add_id : public test_case { : type_a(type_a), type_b(type_b), n_embd(n_embd), n_experts(n_experts), n_experts_used(n_experts_used), n_token(n_token) {} + // MEASURED. ADD_ID's VJP is the identity, so any disagreement is finite-difference rounding + // and nothing else -- but at the default 1e-4 it flaked about 1 run in 10 at MAA 1.2e-4. + // + // worst FD noise, 12 runs 1.2e-4 + // scale the VJP by 2 0.33 (an identity's asymm under a 2x error is exactly 1/3) + // + // 1e-3 sits 8x above the noise and 330x below a real defect. + double max_maa_err() override { + return 1e-3; + } + ggml_tensor * build_graph(ggml_context * ctx) override { ggml_tensor * a = ggml_new_tensor_3d(ctx, type_a, n_embd, n_experts_used, n_token); ggml_set_name(a, "a"); @@ -5788,6 +6169,14 @@ struct test_concat : public test_case { // // A weighted sum, whose weights are not all equal, makes each element of dL/d(out) distinct -- // so a source that is handed the wrong slab is handed visibly wrong numbers. + // MEASURED (S1-37): CONCAT's VJP is a pure routing operation, so any disagreement is + // finite-difference rounding -- but it flaked about 1 run in 15 at MAA ~2e-4 against the + // default 1e-4. A VJP that dropped or transposed a slab measures order 1. 1e-3 sits 5x above + // the noise and ~1000x below a real defect. + double max_maa_err() override { + return 1e-3; + } + ggml_tensor * grad_loss(ggml_context * ctx, ggml_tensor * out) override { ggml_tensor * w = ggml_new_tensor(ctx, GGML_TYPE_F32, GGML_MAX_DIMS, out->ne); ggml_set_name(w, "grad_loss_weights"); @@ -7172,10 +7561,20 @@ struct test_cross_entropy_loss : public test_case { return out; } + // MEASURED (S1-37): worst FD noise 3.8e-3 over 6 runs once conditioned; a VJP scaled by 2 + // measures 0.33. This op was also passing on a NaN before the metric fix -- its softmax is + // one-hot at +-100, so most gradient elements were exactly zero. + double max_maa_err() override { + return 5e-2; + } + void initialize_tensors(ggml_context * ctx) override { // For larger abs. diffs between logits softmax is more linear, therefore more precise num. gradients. for (ggml_tensor * t = ggml_get_first_tensor(ctx); t != NULL; t = ggml_get_next_tensor(ctx, t)) { - init_tensor_uniform(t, -100.0f, 100.0f); + // Same story (S1-37): at +-100 the softmax is one-hot to float precision, so most + // gradient elements are EXACTLY zero -- which used to mean NaN and a free pass. + const float lim = mode == MODE_GRAD ? 3.0f : 100.0f; + init_tensor_uniform(t, -lim, lim); } } diff --git a/tests/test-glu-back.cpp b/tests/test-glu-back.cpp new file mode 100644 index 00000000000..ae408258e8a --- /dev/null +++ b/tests/test-glu-back.cpp @@ -0,0 +1,162 @@ +// GLU_BACK's numerics oracle (learning-llamas, S1-28). +// +// MODE_GRAD cannot check this op. mean_abs_asymm divides by (gn + ga), not (|gn| + |ga|), so a +// near-zero gradient element sends the ratio to infinity -- and every GLU gradient is a PRODUCT +// (dx = dy * g * act'(x)), so near-zero elements are ordinary. Measured, the FD "noise" reaches +// MAA 0.80 while a genuinely broken kernel measures 0.18-0.60. They OVERLAP: no tolerance +// separates them. See test_glu::max_maa_err in test-backend-ops.cpp, and ticket S1-37. +// +// So this is the check that actually validates the derivatives: each variant's analytic VJP against +// a FLOAT64 central difference of ggml's OWN scalar forwards (transcribed from ggml-cpu/vec.h -- +// not from a paper, because if ggml's gelu uses the tanh approximation then the correct derivative +// is the derivative OF THAT). +// +// It catches what MODE_GRAD cannot: a wrong coefficient. Perturbing GEGLU's tanh argument by 5% is +// invisible to MODE_GRAD and fails here immediately. + +// GLU_BACK's analytic derivatives vs a FLOAT64 central difference of ggml's OWN scalar forwards. +// This is the oracle MODE_GRAD cannot be: it is exact to ~1e-7 and it catches a wrong coefficient. +#include "ggml.h" +#include "ggml-cpu.h" +#include +#include +#include + +// ggml's scalar forwards, in double, transcribed from ggml/src/ggml-cpu/vec.h. +static const double A = 0.044715, S2PI = 0.79788456080286535587989211986876, S2I = 0.70710678118654752440084436210484; +static double act(int op, double x, double alpha, double limit) { + switch (op) { + case GGML_GLU_OP_REGLU: return x > 0 ? x : 0; + case GGML_GLU_OP_SWIGLU: return x/(1.0+exp(-x)); + case GGML_GLU_OP_GEGLU: return 0.5*x*(1.0+tanh(S2PI*x*(1.0+A*x*x))); + case GGML_GLU_OP_GEGLU_ERF: return 0.5*x*(1.0+erf(x*S2I)); + case GGML_GLU_OP_GEGLU_QUICK: return x*(1.0/(1.0+exp(-1.702*x))); + } + (void)alpha; (void)limit; return 0; +} +// full forward y = f(x, g) +static double fwd(int op, double x, double g, double alpha, double limit) { + if (op == GGML_GLU_OP_SWIGLU_OAI) { + double xc = x < limit ? x : limit; + double gc = g > limit ? limit : (g < -limit ? -limit : g); + return (xc/(1.0+exp(alpha*(-xc)))) * (gc + 1.0); + } + return act(op, x, alpha, limit) * g; +} + + +// A TRANSPOSED grad must give the same answer as a contiguous one carrying the same values. +// +// ggml's autodiff produces non-contiguous grads routinely (the MUL_MAT backward passes +// ggml_transpose(grad) straight into ggml_out_prod), and no test-backend-ops case ever builds one -- +// so this bug is invisible to MODE_GRAD by construction. It was found in the S1-26/S1-27 MoE kernels +// by adversarial review, and GLU_BACK had it too: the same logical grad in two layouts disagreed by +// 1.30. Read src->nb[0]; never index a float*. +static int check_transposed_grad(void) { + const int NC = 8, NR = 4; + float out[2][2*8*4]; + + for (int pass = 0; pass < 2; ++pass) { + struct ggml_init_params ip = { 16*1024*1024, NULL, false }; + struct ggml_context * ctx = ggml_init(ip); + srand(9); + + struct ggml_tensor * a = ggml_new_tensor_2d(ctx, GGML_TYPE_F32, NC, NR); + struct ggml_tensor * b = ggml_new_tensor_2d(ctx, GGML_TYPE_F32, NC, NR); + for (int i = 0; i < NC*NR; ++i) { + ((float*)a->data)[i] = 2.f*(rand()/(float)RAND_MAX) - 1.f; + ((float*)b->data)[i] = 2.f*(rand()/(float)RAND_MAX) - 1.f; + } + float g[8*4]; + for (int i = 0; i < NC*NR; ++i) g[i] = 2.f*(rand()/(float)RAND_MAX) - 1.f; + + struct ggml_tensor * dy; + if (pass == 0) { + dy = ggml_new_tensor_2d(ctx, GGML_TYPE_F32, NC, NR); + for (int i = 0; i < NC*NR; ++i) ((float*)dy->data)[i] = g[i]; + } else { + struct ggml_tensor * base = ggml_new_tensor_2d(ctx, GGML_TYPE_F32, NR, NC); + for (int r = 0; r < NR; ++r) for (int k = 0; k < NC; ++k) ((float*)base->data)[k*NR+r] = g[r*NC+k]; + dy = ggml_transpose(ctx, base); // ne=[NC,NR], nb[0] = NR*4, NOT 4 + } + + struct ggml_tensor * o = ggml_glu_back(ctx, dy, a, b, GGML_GLU_OP_GEGLU, false, 0.f, 0.f); + struct ggml_cgraph * gf = ggml_new_graph(ctx); + ggml_build_forward_expand(gf, o); + ggml_graph_compute_with_ctx(ctx, gf, 2); + for (int i = 0; i < 2*NC*NR; ++i) out[pass][i] = ((float*)o->data)[i]; + ggml_free(ctx); + } + + double worst = 0; + for (int i = 0; i < 2*NC*NR; ++i) worst = fmax(worst, fabs(out[0][i] - out[1][i])); + printf(" %-12s transposed grad vs contiguous: max diff %.3e %s\n", + "STRIDES", worst, worst == 0.0 ? "identical" : "*** MISREADS A TRANSPOSED GRAD ***"); + return worst != 0.0; +} + +int main(void) { + const int ops[] = {GGML_GLU_OP_REGLU, GGML_GLU_OP_SWIGLU, GGML_GLU_OP_GEGLU, + GGML_GLU_OP_GEGLU_ERF, GGML_GLU_OP_GEGLU_QUICK, GGML_GLU_OP_SWIGLU_OAI}; + const char * nm[] = {"REGLU","SWIGLU","GEGLU","GEGLU_ERF","GEGLU_QUICK","SWIGLU_OAI"}; + const double alpha = 1.702, limit = 7.0; + const int N = 64, NR = 5; + int bad = 0; + + for (int oi = 0; oi < 6; ++oi) { + const int op = ops[oi]; + struct ggml_init_params ip = { 32*1024*1024, NULL, false }; + struct ggml_context * ctx = ggml_init(ip); + + struct ggml_tensor * a = ggml_new_tensor_2d(ctx, GGML_TYPE_F32, N, NR); // gate half + struct ggml_tensor * b = ggml_new_tensor_2d(ctx, GGML_TYPE_F32, N, NR); // linear half + struct ggml_tensor * dy = ggml_new_tensor_2d(ctx, GGML_TYPE_F32, N, NR); + srand(5 + oi); + for (int i = 0; i < N*NR; ++i) { + ((float*)a->data)[i] = 6.0f*(rand()/(float)RAND_MAX) - 3.0f; // spans kinks & tails + ((float*)b->data)[i] = 6.0f*(rand()/(float)RAND_MAX) - 3.0f; + ((float*)dy->data)[i] = 2.0f*(rand()/(float)RAND_MAX) - 1.0f; + } + + struct ggml_tensor * out = ggml_glu_back(ctx, dy, a, b, (enum ggml_glu_op) op, false, (float)alpha, (float)limit); + struct ggml_cgraph * gf = ggml_new_graph(ctx); + ggml_build_forward_expand(gf, out); + ggml_graph_compute_with_ctx(ctx, gf, 2); + + double worst_dx = 0, worst_dg = 0, scale_x = 0, scale_g = 0; + for (int r = 0; r < NR; ++r) { + const float * dxk = (const float*)out->data + r*2*N; // half 0 = d_a + const float * dgk = (const float*)out->data + r*2*N + N; // half 1 = d_b + for (int k = 0; k < N; ++k) { + double x = ((float*)a->data)[r*N+k], g = ((float*)b->data)[r*N+k], d = ((float*)dy->data)[r*N+k]; + // REGLU's kink and SWIGLU_OAI's clamps are non-differentiable points: skip a + // neighbourhood of them, since the analytic VJP uses a subgradient there by design. + const double h = 1e-5; + if (op == GGML_GLU_OP_REGLU && fabs(x) < 1e-3) continue; + if (op == GGML_GLU_OP_SWIGLU_OAI && (fabs(x-limit) < 1e-3 || fabs(fabs(g)-limit) < 1e-3)) continue; + + double fd_x = d * (fwd(op, x+h, g, alpha, limit) - fwd(op, x-h, g, alpha, limit)) / (2*h); + double fd_g = d * (fwd(op, x, g+h, alpha, limit) - fwd(op, x, g-h, alpha, limit)) / (2*h); + // Absolute error, normalized by the tensor's SCALE -- not by the element's own + // magnitude. silu' has a zero crossing at x ~ -1.278 and gelu' has one too, so a + // per-element relative error divides by ~0 there and explodes on an exact kernel. + // That is the same near-zero trap that makes MODE_GRAD's metric unusable; a + // reference must not repeat it. + if (fabs(fd_x - dxk[k]) > worst_dx) worst_dx = fabs(fd_x - dxk[k]); + if (fabs(fd_g - dgk[k]) > worst_dg) worst_dg = fabs(fd_g - dgk[k]); + if (fabs(fd_x) > scale_x) scale_x = fabs(fd_x); + if (fabs(fd_g) > scale_g) scale_g = fabs(fd_g); + } + } + double rx = worst_dx/(scale_x > 0 ? scale_x : 1), rg = worst_dg/(scale_g > 0 ? scale_g : 1); + int ok = rx < 1e-5 && rg < 1e-5; + if (!ok) bad++; + printf(" %-12s worst err / scale: d_gate %.2e d_linear %.2e %s\n", + nm[oi], rx, rg, ok ? "exact" : "*** WRONG ***"); + ggml_free(ctx); + } + bad += check_transposed_grad(); + + printf(bad ? "\n%d CHECKS FAILED\n" : "\nall six variants match ggml's own forwards, and strides are honoured\n", bad); + return bad != 0; +}