--- a/nnCostFunction.m
+++ b/nnCostFunction.m
@@ -42,19 +42,15 @@
       
       ## The regularisation term has to exclude the first column of the Thetas,
       ## because we don't regularise the bias nodes.
-      + lambda*(sum (Theta1(:, 2:end)(:).^2)                       \
-                + sum (Theta2(:, 2:end)(:).^2))/(2*m);
+      + lambda*(sumsq (Theta1(:, 2:end)(:))                         \
+                + sumsq (Theta2(:, 2:end)(:)))/(2*m);
 
   ## Backprop
   delta3 = a3 - y_idx;
   delta2 = (delta3*Theta2)(:, 2:end) .* sigmoidGradient (z2);
 
-  Theta2_grad = sum (bsxfun (@times, permute (delta3, [2, 3, 1]),
-                             permute ([one_vec, a2], [3, 2, 1])),
-                     3)/m;
-  Theta1_grad = sum (bsxfun (@times, permute (delta2, [2, 3, 1]),
-                             permute ([one_vec, a1], [3, 2, 1])),
-                     3)/m;
+  Theta2_grad = delta3' * [one_vec, a2] / m;
+  Theta1_grad = delta2' * [one_vec, a1] / m;
 
   ## Add regularisation terms
   Theta2_grad(:, 2:end) += Theta2(:, 2:end)*lambda/m;