Implementation of tensor product formula

This might help you get an idea:

n = 4;
BlockRandom[SeedRandom[42, Method -> "Legacy"]; (* for reproducibility *)
            cpts = Table[{i, j, RandomReal[{-1, 1}]}, {i, n + 1}, {j, n + 1}]];

GraphicsRow[{ParametricPlot3D[BezierFunction[cpts][u, v], {u, 0, 1}, {v, 0, 1},
                              Evaluated -> True], (* built-in function *)
             ParametricPlot3D[Evaluate[Fold[#2.#1 &, cpts, (* using dot-products *)
                                            {BernsteinBasis[n, Range[0, n], u],
                                             BernsteinBasis[n, Range[0, n], v]}]],
                              {u, 0, 1}, {v, 0, 1}]}]

You can of course replace BernsteinBasis[] with your own Bernstein[]; no need for Hold[] trickery!

Update

Fully compiling the code to C makes it as fast as the built-in:

cBernstein = 
 Compile @@ (Hold[{n, {j, _Real, 1}, u}, Table[expr u^i (1 - u)^(n - i), {i, j}], 
     CompilationTarget -> C] /. expr -> FunctionExpand[Binomial[n, i]])

BezierSurface4 = 
  With[{cBernstein = cBernstein}, 
   With[{AllBasis = Function[{deg, u0}, cBernstein[deg, Range[0, deg], u0]]}, 
    Compile[{{pts, _Real, 3}, u, v}, Module[{m, n}, {m, n} = Dimensions[pts, 2];
      
      With[{row = AllBasis[m - 1, u], col = AllBasis[n - 1, v]}, 
       row.Transpose[pts, {1, 3, 2}].col]], 
     RuntimeOptions -> {"EvaluateSymbolically" -> False}, 
     CompilationOptions -> {"InlineCompiledFunctions" -> True}, 
     CompilationTarget -> C]]];

f = BezierFunction[cpts];
ParametricPlot3D[f[u, v], {u, 0, 1}, {v, 0, 1}]; // AbsoluteTiming
(* {0.018014, Null} *)

ParametricPlot3D[BezierSurface4[cpts, u, v], {u, 0, 1}, {v, 0, 1}]; // AbsoluteTiming

(* {0.018006, Null} *)

---

Old Answer

Compile and Transpose and Listable will help:

cBernstein = 
 Compile[{n, {i, _Real, 1}, u}, 
  Evaluate[FunctionExpand[Binomial[n, i]] u^i (1 - u)^(n - i)]]

BezierSurface3[pts_, u_?NumericQ, v_?NumericQ] := 
 Module[{m, n, AllBasis}, {m, n} = Dimensions[pts, 2];
  AllBasis = Function[{deg, u0}, cBernstein[deg, Range[0, deg], u0]];
  With[{row = AllBasis[m - 1, u], col = AllBasis[n - 1, v]}, 
   row.Transpose[pts, {1, 3, 2}].col]]

ParametricPlot3D[BezierSurface1[cpts, u, v], {u, 0, 1}, {v, 0, 1}]; // AbsoluteTiming
(* {0.219170, Null} *)
ParametricPlot3D[BezierSurface3[cpts, u, v], {u, 0, 1}, {v, 0, 1}]; // AbsoluteTiming
(* {0.098571, Null} *)