I think the Bonami, Salvagnin and Tramontani paper actually gives a fair amount of insight into this, even if we cannot know every detail of the current proprietary implementation. The implementation described there is not simply the same callback algorithm running in C. CPLEX identifies the decomposition, presolves the complete model, decomposes the presolved model, runs a stabilized Benders loop on the LP relaxation, and then starts branch-and-cut with Benders cuts separated as lazy constraints. Full model presolve and native access to LP bases and cut management are difficult to reproduce when the master and worker are constructed separately.
Walter’s point about the root bound is therefore where I would start. If the callback implementation enters the tree with a weaker bound than CPLEX obtains after its initial Benders loop, the larger tree is not surprising. I would compare the final root bounds and the time and number of cuts needed to reach them before comparing the rest of the search.
There is also a paper that is especially relevant to this example. Fischetti, Ljubić and Sinnl describe a callback Benders implementation for essentially the same capacitated facility location problem in Benders decomposition without separability. It contains several refinements that seem more directly applicable than generic Benders advice.
One thing that is easy to miss in the posted formulation is that the aggregate capacity constraint already guarantees feasible recourse. Assignments are continuous and every customer can be served by every facility, so the worker should not generate feasibility cuts or Farkas rays. The relevant cut quality issue is the degeneracy of the transportation LP and the choice between alternative optimal dual solutions.
The CFL paper handles this by recomputing the optimality cut coefficients rather than simply accepting the reduced costs returned by the LP solver. It fixes the dual multipliers for the customer assignment constraints and obtains each facility coefficient from a continuous knapsack problem. These knapsacks are cheap to solve by sorting. This reduces the arbitrary coefficient choice caused by dual degeneracy and produced more stable cuts in their implementation. It seems a more targeted improvement for this problem than feasibility cut normalization.
They also keep the transportation worker persistent and solve it using dual simplex so that each solve reuses the previous basis. Their in-out procedure deliberately keeps successive separation points fairly close in the linear CFL case, making the basis warm start more effective. This is worth considering if the current in-out implementation chooses points without accounting for how much the worker changes between solves.
Their separation policy in the tree is also quite selective. They limit repeated separation at a node and stop calling the worker at fractional points when separation is taking too much time or is unlikely to prune the node. Integer candidates are still always checked through the lazy callback, so correctness is preserved. This directly addresses the problem described in the question, where the callback implementation generates a very large number of cuts.
Another useful idea is their restart procedure. After processing the root, they retain the useful callback cuts, add them to the master as ordinary constraints, and resolve the root before starting the final tree search. This allows presolve, variable fixing and CPLEX’s internal cuts to use the accumulated Benders information from the beginning. It also gives a practical reason to store callback cuts separately, as discussed earlier in the thread.
I would add y[i,j] <= x[i] as well. These inequalities are redundant when x is binary, but they strengthen the LP relaxation and give the generalized bound structure discussed in the CPLEX paper. In a persistent worker they can be handled as bound changes, provided their dual contribution is included when constructing the cut.
On multithreading, the limitation is not Benders itself. A native CPLEX callback can run during parallel MIP search if it is thread safe. The CFL implementation creates one persistent worker clone for each CPLEX thread, so simultaneous callback calls do not share a worker model or LP basis. A JuMP implementation would need the same architecture, together with CPLEX.jl support for callbacks arriving on multiple solver threads, thread local worker models and careful handling of shared cut state. CPLEX.jl currently documents its callbacks for single threaded use, so this cannot safely be reproduced merely by setting Threads > 1. It does not explain the single thread comparison here, but it does explain an architectural advantage of CPLEX’s internal implementation.
There are also a few benchmarking issues in the attached script. The seed argument is unused, the transportation cost matrix is transposed relative to its later indexing, and only cover, flow cover and MIR cuts are disabled rather than all other CPLEX cut families. These should be corrected before relying on the cut and node counts. CPXPARAM_Benders_Strategy = 1 is not itself a problem and is probably the fairest comparison when both methods use the same partition.
My reading is that the difference is likely to come from the combination of full model presolve, a stronger root phase, better treatment of degenerate optimality cuts, basis aware worker solves and more selective separation in the tree. It is not simply that CPLEX generates the same cuts faster.