diff --git a/apps/bgu/bgu_generator.cpp b/apps/bgu/bgu_generator.cpp
index 653e5ed87c5a..e375241c56d4 100644
--- a/apps/bgu/bgu_generator.cpp
+++ b/apps/bgu/bgu_generator.cpp
@@ -395,39 +395,20 @@ class BGU : public Generator<BGU> {
             b(2, 2) = blurx(x, y, z, 20);
             b(3, 2) = blurx(x, y, z, 21);
 
-            // Regularize by pushing the solution towards the average gain
-            // in this cell = (average output luma + eps) / (average input luma + eps).
-            const float lambda = 1e-6f;
-            const float epsilon = 1e-6f;
-
-            // The bottom right entry of A is a count of the number of
-            // constraints affecting this cell.
-            Expr N = A(3, 3);
-
-            // The last row of each matrix is the sum of input and output
-            // RGB values for the pixels affecting this cell. Instead of
-            // dividing them by N+1 to get averages, we'll multiply
-            // epsilon by N+1. This saves two divisions.
-            Expr output_luma = b(3, 0) + 2 * b(3, 1) + b(3, 2) + epsilon * (N + 1);
-            Expr input_luma = A(3, 0) + 2 * A(3, 1) + A(3, 2) + epsilon * (N + 1);
-            Expr gain = output_luma / input_luma;
-
-            // Add lambda and lambda*gain to the diagonal of the
-            // matrices. The matrices are sums/moments rather than
-            // means/covariances, so just like above we need to multiply
-            // lambda by N+1 so that it's equivalent to adding a constant
-            // to the diagonal of a covariance matrix. Otherwise it does
-            // nothing in cells with lots of linearly-dependent
-            // constraints.
-            Expr weighted_lambda = lambda * (N + 1);
-            A(0, 0) += weighted_lambda;
-            A(1, 1) += weighted_lambda;
-            A(2, 2) += weighted_lambda;
-            A(3, 3) += weighted_lambda;
-
-            b(0, 0) += weighted_lambda * gain;
-            b(1, 1) += weighted_lambda * gain;
-            b(2, 2) += weighted_lambda * gain;
+            // Regularize it with 1/10th of a sample that pulls the result towards the identity function.
+            // Regions in the grid that are well populated will have way more than 1/10th of a sample.
+            // The original paper on BGU had a more complex regularization scheme, but the regularization
+            // logic was backwards: when a cell has fewer samples, it got less regularization; and when
+            // the cell has a lot of samples, it got regularized a lot.
+            const float lambda = 1e-1f;
+            A(0, 0) += lambda;
+            A(1, 1) += lambda;
+            A(2, 2) += lambda;
+            A(3, 3) += lambda;
+
+            b(0, 0) += lambda;
+            b(1, 1) += lambda;
+            b(2, 2) += lambda;
 
             // Now solve Ax = b
             Matrix<3, 4> result = transpose(solve_symmetric(A, b, line, x, using_autoscheduler(), get_target()));