# Some question on column generation, branch and price, Coluna

**URL:** <https://discourse.julialang.org/t/some-question-on-column-generation-branch-and-price-coluna/131084>\
**Category:** Optimization (Mathematical)\
**Tags:** question, dantzig-wolfe\
**Created:** [July 27, 2025, 2:33am UTC](https://discourse.julialang.org/t/some-question-on-column-generation-branch-and-price-coluna/131084 "2025-07-27T02:33:06Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [July 27, 2025, 2:33am UTC](https://discourse.julialang.org/t/some-question-on-column-generation-branch-and-price-coluna/131084/1 "2025-07-27T02:33:06Z")

</div>

Since I’m studying this topic recently. I wonder if there is anyone who is interested and can hence discuss with me.

#### Benders

 ![image](https://global.discourse-cdn.com/julialang/original/3X/8/f/8f7cd8216c0b424ddac000d4695a5e245f10d9af.jpeg)  
~~I don’t think this structure is going to be handled by Benders decomposition.  
For this structure, a more appropriate method is still Dantzig-Wolfe decomposition, if you add copy variables and copy constraints and then relax all complicating constraints to achieve blockwise decomposition.~~

#### Diving heuristic

As it is known that branch-and-bound becomes slow and hard, in practical applications we may not have enough time to exactly conduct a branch-and-price algorithm. Therefore some heuristics are entailed. I’m interested in this diving heuristic marked below

 ![image](https://global.discourse-cdn.com/julialang/original/3X/e/8/e86cc7b7c35d98a8cd31b611e48db07eea10055f.jpeg)  
Could anyone offer some more details about this?  
PS I wrote a simple prototype myself

> **code**
>
> ```julia
> import JuMP, Gurobi; GRB_ENV = Gurobi.Env();
> import JuMP.value as ı
> import LinearAlgebra.⋅ as ⋅
> import LinearAlgebra.norm as norm
> import Random
> 
> function is_binary_bold(x) return abs(.5 - x) > (.5 - 5e-12) end; # a bold version of `is_binary`
> function prev(t, d, T)
> n = t - d
> return n < 1 ? T+n : n
> end; function gen_a_D_house(ID, DˈUB) return (
> ID = ID,
> D = rand(0:DˈUB, T), # base demand
> P_BUS = rand(12:24), # unidirectional due to no renewables
> P_EV = [2, rand(3:7)], # MIN and MAX charging power
> E_EV = rand(10:39),
> DW = rand(1:4, rand(2:5)), # cycle vector   
> ) end; function gen_a_G_house(ID) return (
> ID = ID,
> G = rand(0:15, T), # 🌞 renewables generation
> D = rand(2:9, T), # base demand
> P_BUS = rand(12:14), # bidirectional due to 🌞
> P_EV = [2, rand(3:10)], # MIN and MAX charging power
> E_EV = rand(10:39),
> DW = rand(1:7, rand(2:7)), # cycle vector
> ) end; function find_a_frac_block(χ)
> for j = Random.shuffle(1:J)
> ks = haskey(χ[j], :bLent) ? [:bLent, :bEV, :bDW] : [:bEV, :bDW]
> for k = ks, x = χ[j][k]
> is_binary_bold(ı(x)) || return j # [early return]
> end
> end; return 0 # all solutions are now agreeably integer
> end; function solve_restricted_primal(B; mip = false)
> model = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); # restricted primal
> if mip
> JuMP.set_attribute(model, "MIPGap", 0)
> JuMP.set_attribute(model, "MIPGapAbs", 0)
> else
> JuMP.set_silent(model) # This is an LP, used to generate fractional trial solutions
> end
> JuMP.@variable(model, λ[j = 1:J, v = 1:length(B[j])] ≥ 0); JuMP.@constraint(model, [j = 1:J], sum(λ[j, :]) == 1);
> χ = NamedTuple[(; (k => similar(v, JuMP.VariableRef) for (k, v) = pairs(first(keys(B[j]))))...) for j = 1:J]; # The "mean" vector
> for j = 1:J, k = keys(χ[j]), i = eachindex(χ[j][k])
> χ[j][k][i] = scalar_mean = JuMP.@variable(model)
> JuMP.@constraint(model, scalar_mean == sum(λ[j, v]nt[k][i] for (v, nt) = enumerate(keys(B[j]))))
> mip || continue
> k ∈ [:bLent, :bEV, :bDW] && JuMP.set_integer(scalar_mean)
> end
> JuMP.@expression(model, pAgr[t = 1:T], sum(haskey(χ[j], :pBus) ? χ[j][:pBus][t] : χ[j][:pBus1][t]+χ[j][:pBus2][t] for j = 1:J));
> JuMP.@constraint(model, β[t = 1:T], P_AGR ≥ pAgr[t]); # 🟡 linking constr
> # ⚠️ we do not explicitly enforce the "no reverse flow" constraint, but will test it post-optimizing
> JuMP.@objective(model, Min, MarketPrice ⋅ pAgr);
> JuMP.optimize!(model); JuMP.assert_is_solved_and_feasible(model; allow_local = false)
> return χ
> end; function initialize_inn(D)
> inn = Vector{JuMP.Model}(undef, J); # the intrablock MIP subproblems
> for (j, block) = enumerate(D)
> model = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)) # MIP of the j-th block
> JuMP.set_silent(model)
> JuMP.set_attribute(model, "Threads", 1)
> JuMP.set_attribute(model, "TimeLimit", 3)
> JuMP.set_attribute(model, "MIPGap", 0)
> JuMP.set_attribute(model, "MIPGapAbs", 0)
> if length(block) == 2 # [pBus1; pBus2; pLent; vec(pEV); bLent; vec(bEV); vec(bDW)]
> h1, h2 = block
> JuMP.@variable(model, bLent[1:T], Bin) # 0 is lending
> JuMP.@variable(model, pLent[1:T])
> JuMP.@constraint(model, [t = 1:T], pLent[t] ≤ (1 - bLent[t])h1[:P_EV][end])
> JuMP.@constraint(model, [t = 1:T], (1 - bLent[t])h1[:P_EV][begin] ≤ pLent[t])
> JuMP.@variable(model, bEV[1:T, 1:2], Bin)
> JuMP.@variable(model, pEV[1:T, 1:2])
> JuMP.@constraint(model, [t = 1:T, h = 1:2], pEV[t, h] ≤ bEV[t, h]block[h][:P_EV][end])
> JuMP.@constraint(model, [t = 1:T, h = 1:2], bEV[t, h]block[h][:P_EV][begin] ≤ pEV[t, h])
> JuMP.@constraint(model, [t = 1:T, h = 1:2], bEV[t, h] ≤ bLent[t])
> JuMP.@constraint(model, sum( view(pEV, :, 1)) ≥ h1[:E_EV])
> JuMP.@constraint(model, sum(pLent + view(pEV, :, 2)) ≥ h2[:E_EV])
> JuMP.@variable(model, bDW[1:T, 1:2], Bin)
> JuMP.@expression(model, pDW[t = 1:T, h = 1:2], sum(bDW[prev(t, φ-1, T), h]P for (φ, P) = enumerate(block[h][:DW])))
> JuMP.@constraint(model, [h = 1:2], sum(view(bDW, :, h)) ≥ true)
> JuMP.@variable(model, -h1[:P_BUS] ≤ pBus1[1:T] ≤ h1[:P_BUS]) # 🟡
> JuMP.@variable(model, false ≤ pBus2[1:T] ≤ h2[:P_BUS]) # 🟡
> JuMP.@constraint(model, pBus1 + h1[:G] ≥ pLent + h1[:D] + view(pEV, :, 1) + view(pDW, :, 1))
> JuMP.@constraint(model, pBus2 ≥ h2[:D] + view(pEV, :, 2) + view(pDW, :, 2))
> JuMP.@expression(model, prim_obj, MarketPrice ⋅ (pBus1 + pBus2)) # the additional penalty is merely β-coefficiented
> elseif length(block) == 1 # [pBus; pEV; bEV; bDW]
> h = block[1]
> JuMP.@variable(model, bEV[1:T], Bin)
> JuMP.@variable(model, pEV[1:T])
> JuMP.@constraint(model, [t = 1:T], pEV[t] ≤ bEV[t]h[:P_EV][end])
> JuMP.@constraint(model, [t = 1:T], bEV[t]h[:P_EV][begin] ≤ pEV[t])
> JuMP.@constraint(model, sum(pEV) ≥ h[:E_EV])
> JuMP.@variable(model, bDW[1:T], Bin)
> JuMP.@expression(model, pDW[t = 1:T], sum(bDW[prev(t, φ-1, T)]P for (φ, P) = enumerate(h[:DW])))
> JuMP.@constraint(model, sum(bDW) ≥ true)
> JuMP.@variable(model, false ≤ pBus[1:T] ≤ h[:P_BUS]) # 🟡
> JuMP.@constraint(model, pBus ≥ h[:D] + pEV + pDW)
> JuMP.@expression(model, prim_obj, MarketPrice ⋅ pBus) # the additional penalty is merely β-coefficiented
> end; inn[j] = model # save
> end # not include: penalty objective, non-block linking constr
> return inn
> end; function initialize_out(an_UB) # an_UB ⚠️ should be valid
> model = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)) # outer in DantzigWolfe decomposition
> JuMP.set_silent(model)
> JuMP.@variable(model, β[1:T] ≥ 0) # multiplier vector of the 🟡 linking constr
> JuMP.@variable(model, θ[1:J]) # This is a concise formulation---we omit the `pi` variable
> JuMP.@expression(model, common, -sum(β)P_AGR)
> JuMP.@expression(model, out_obj_tbMax, common + sum(θ))
> JuMP.@objective(model, Max, out_obj_tbMax)
> JuMP.@constraint(model, out_obj_tbMax ≤ an_UB)
> return model, θ, β
> end; function reset_mj_obj!(mj, Β)
> ( haskey(mj, :pBus) ? 
> JuMP.@objective(mj, Min, mj[:prim_obj] + Β ⋅ mj[:pBus]) :
> JuMP.@objective(mj, Min, mj[:prim_obj] + Β ⋅ (mj[:pBus1] + mj[:pBus2])) )
> end; function build_con(j, mj, θ, β)
> haskey(mj, :pBus) && return JuMP.@build_constraint(θ[j] ≤ ı(mj[:prim_obj]) + β ⋅ ı.(mj[:pBus]))
> return JuMP.@build_constraint(θ[j] ≤ ı(mj[:prim_obj]) + β ⋅ (ı.(mj[:pBus1]) + ı.(mj[:pBus2])))
> end; function get_full_vec(nt)
> haskey(nt, :pBus) && return [nt[:pBus]; nt[:pEV]; nt[:bEV]; nt[:bDW]]
> return [nt[:pBus1]; nt[:pBus2]; nt[:pLent]; vec(nt[:pEV]); nt[:bLent]; vec(nt[:bEV]); vec(nt[:bDW])]
> end; function get_bool_part_vec(nt)
> haskey(nt, :bLent) && return [nt[:bLent]; vec(nt[:bEV]); vec(nt[:bDW])]
> return [nt[:bEV]; nt[:bDW]]
> end; function get_full_num_tuple(nt)
> haskey(nt, :pBus) && return (
> pBus = ı.(nt[:pBus]),
> pEV = ı.(nt[:pEV]),
> bEV = round.(Bool, ı.(nt[:bEV])),
> bDW = round.(Bool, ı.(nt[:bDW])) )
> return (
> pBus1 = ı.(nt[:pBus1]),
> pBus2 = ı.(nt[:pBus2]),
> pLent = ı.(nt[:pLent]),
> pEV = ı.(nt[:pEV]),
> bLent = round.(Bool, ı.(nt[:bLent])),
> bEV = round.(Bool, ı.(nt[:bEV])),
> bDW = round.(Bool, ı.(nt[:bDW]))
> )
> end; function parallel_CG!(B, model, θ, β)
> for ite = 1:typemax(Int)
> JuMP.optimize!(model); JuMP.assert_is_solved_and_feasible(model; allow_local = false, dual = true)
> @info "ite = $ite, bound is $(JuMP.objective_bound(model))"
> Β = ı.(β) # ✂️
> out_upd_vec .= false
> Threads.@threads for j = 1:J
> Θ_j = ı(θ[j]); mj = inn[j] # ✂️
> reset_mj_obj!(mj, Β)
> JuMP.optimize!(mj); JuMP.assert_is_solved_and_feasible(mj; allow_local = false)
> Δ = Θ_j - JuMP.objective_value(mj)
> if Δ > 1e-11
> @info "Δ = $Δ, at j = $j"
> local intFullvertex = get_full_num_tuple(mj)
> if !haskey(B[j], intFullvertex)
> local con = build_con(j, mj, θ, β)
> local cut = @lock my_lock JuMP.add_constraint(model, con)
> B[j][intFullvertex] = cut
> out_upd_vec[j] = true
> end
> end
> end
> any(out_upd_vec) && continue
> break
> end
> end; function condense!(B, model) # ⚠️🔴 This post-procedure could be harmful
> # e.g. lose the chance to recover a integer-feasible primal solution
> # or, even if you can recover one, its quality is inferior than the one
> # recovered were this post-procedure unexecuted
> for j = 1:J, (vertex, cut) = collect(B[j]) # 🚰 condensing
> JuMP.dual(cut) == 0 || continue
> pop!(B[j], vertex)
> JuMP.delete(model, cut) # delete one-by-one is already fast, so we don't bother reoptimizing
> JuMP.optimize!(model); JuMP.assert_is_solved_and_feasible(model; allow_local = false, dual = true)
> end; @info "Lagrangian bound is $(JuMP.objective_bound(model)) after condensing"
> end; function restrict_a_frac_block!(B, model, inn)
> χ = solve_restricted_primal(B)
> j = find_a_frac_block(χ)
> if j == 0
> @info "▶▶▶ Now all fractional solutions has become integer"
> return true # you need to stop
> else
> @info "▶ eliminating fractionality in block $j"
> end
> mj = inn[j] # the block that we want to make 🌳 branching decision
> ⅁ = ı.(get_full_vec(χ[j])) # the concatenated fractional trial solution 
> v = get_full_vec(mj) # create a vec of decisions aligning with it
> o = JuMP.@variable(mj) # temporary epi-variable
> c1 = JuMP.@constraint(mj, v - ⅁ .≤ o) 
> c2 = JuMP.@constraint(mj, ⅁ - v .≤ o) # vector of constrs
> JuMP.@objective(mj, Min, o) # set-and-solve routine
> JuMP.optimize!(mj); JuMP.assert_is_solved_and_feasible(mj; allow_local = false) # the small MIP
> iv = get_full_num_tuple(mj) # integer-feasible vertex in NamedTuple format
> JuMP.delete.(mj, c1);
> JuMP.delete.(mj, c2);
> JuMP.delete(mj, o);
> iA = get_bool_part_vec(iv); # the Bool Anchor vector
> JuMP.fix.(get_bool_part_vec(mj), iA; force = true) # 🌳 The branching (= fixing) decision
> for (vertex, cut) = collect(B[j])
> get_bool_part_vec(vertex) == iA && continue # keep the compatible ones previously gened
> pop!(B[j], vertex)
> JuMP.delete(model, cut)
> end
> return false # don't need to stop
> end; 
> 
> begin # Data acquisition (D)
> seed = 6787495924352825842 # [interesting case] IP_OPT_VALUE > LAG_DUAL_BOUND, but diving can recover IP_OPT_VALUE!
> # seed = rand(Int);
> Random.seed!(seed);
> T = 24; # global
> MarketPrice = [2, 1, 1, 1, 1, 1, 3, 5, 6, 7, 8, 9, 9, 8, 8, 9, 10, 11, 9, 6, 5, 3, 2, 2];
> P_AGR = 47; # pAgr should be in [0, P_AGR]
> D = Vector{Vector{NamedTuple}}(undef, 0); # the global data
> ID = 0; # each house has a unique ID
> while ID < 4 # build house pairs
> ID += 1; h1 = gen_a_G_house(ID)
> ID += 1; h2 = gen_a_D_house(ID, 10)
> push!(D, NamedTuple[h1, h2])
> end
> while ID < 6 # build house singles
> ID += 1; h = gen_a_D_house(ID, 10)
> push!(D, NamedTuple[h])
> end
> J = length(D)
> @info "a community with $J blocks and $ID houses have been built"
> end;
> 
> if true # 🟪 centralized formulation: use primitive device constraints involving private data
> Χ = Vector{NamedTuple}(undef, J) # record the full standard optimal solution
> model = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)) # monolithic
> JuMP.set_attribute(model, "MIPGap", 0)
> JuMP.set_attribute(model, "MIPGapAbs", 0)
> for (j, block) = enumerate(D) # add intrablock constraints
> if length(block) == 2
> h1, h2 = block # decode
> # EV's part
> bLent = JuMP.@variable(model, [1:T], Bin) # 0 is lending
> pLent = JuMP.@variable(model, [1:T])
> JuMP.@constraint(model, [t = 1:T], pLent[t] ≤ (1 - bLent[t])h1[:P_EV][end])
> JuMP.@constraint(model, [t = 1:T], (1 - bLent[t])h1[:P_EV][begin] ≤ pLent[t])
> bEV = JuMP.@variable(model, [1:T, 1:2], Bin)
> pEV = JuMP.@variable(model, [1:T, 1:2])
> JuMP.@constraint(model, [t = 1:T, h = 1:2], pEV[t, h] ≤ bEV[t, h]block[h][:P_EV][end])
> JuMP.@constraint(model, [t = 1:T, h = 1:2], bEV[t, h]block[h][:P_EV][begin] ≤ pEV[t, h])
> JuMP.@constraint(model, [t = 1:T, h = 1:2], bEV[t, h] ≤ bLent[t])
> JuMP.@constraint(model, sum( view(pEV, :, 1)) ≥ h1[:E_EV])
> JuMP.@constraint(model, sum(pLent + view(pEV, :, 2)) ≥ h2[:E_EV])
> # DW's part
> bDW = JuMP.@variable(model, [1:T, 1:2], Bin)
> pDW = JuMP.@expression(model, [t = 1:T, h = 1:2], sum(bDW[prev(t, φ-1, T), h]P for (φ, P) = enumerate(block[h][:DW])))
> JuMP.@constraint(model, [h = 1:2], sum(view(bDW, :, h)) ≥ true)
> # household balance
> pBus1 = JuMP.@variable(model, [1:T], lower_bound = -h1[:P_BUS], upper_bound = h1[:P_BUS])
> pBus2 = JuMP.@variable(model, [1:T], lower_bound = false, upper_bound = h2[:P_BUS])
> JuMP.@constraint(model, pBus1 + h1[:G] ≥ pLent + h1[:D] + view(pEV, :, 1) + view(pDW, :, 1))
> JuMP.@constraint(model, pBus2 ≥ h2[:D] + view(pEV, :, 2) + view(pDW, :, 2))
> Χ[j] = (
> pBus1 = pBus1,
> pBus2 = pBus2,
> pLent = pLent,
> pEV = pEV,
> bLent = bLent,
> bEV = bEV,
> bDW = bDW
> )
> elseif length(block) == 1
> h = block[1] # decode
> # EV's part
> bEV = JuMP.@variable(model, [1:T], Bin)
> pEV = JuMP.@variable(model, [1:T])
> JuMP.@constraint(model, [t = 1:T], pEV[t] ≤ bEV[t]h[:P_EV][end])
> JuMP.@constraint(model, [t = 1:T], bEV[t]h[:P_EV][begin] ≤ pEV[t])
> JuMP.@constraint(model, sum(pEV) ≥ h[:E_EV])
> # DW's part
> bDW = JuMP.@variable(model, [1:T], Bin)
> pDW = JuMP.@expression(model, [t = 1:T], sum(bDW[prev(t, φ-1, T)]P for (φ, P) = enumerate(h[:DW])))
> JuMP.@constraint(model, sum(bDW) ≥ true)
> # household balance
> pBus = JuMP.@variable(model, [1:T], lower_bound = false, upper_bound = h[:P_BUS])
> JuMP.@constraint(model, pBus ≥ h[:D] + pEV + pDW)
> Χ[j] = (
> pBus = pBus,
> pEV = pEV,
> bEV = bEV,
> bDW = bDW
> )
> end
> end # we opt not to create an explicit public variable
> JuMP.@expression(model, pAgr, sum(haskey(Χ[j], :pBus) ? Χ[j][:pBus] : Χ[j][:pBus1]+Χ[j][:pBus2] for j = 1:J));
> JuMP.@constraint(model, pAgr .≤ P_AGR); # 🟡 linking constr
> # ⚠️ we do not explicitly enforce the "no reverse flow" constraint, but will test it post-optimize
> JuMP.@objective(model, Min, MarketPrice ⋅ pAgr);
> JuMP.optimize!(model); JuMP.assert_is_solved_and_feasible(model; allow_local = false)
> all(ı.(pAgr) .> 0) || error("This instance is not proper, please enlarge load")
> an_UB = JuMP.objective_value(model) + 500
> end; # Note: use this standard procedure to pick seed leading to proper tests
> 
> inn = initialize_inn(D); 
> my_lock = Threads.ReentrantLock();
> model, θ, β = initialize_out(an_UB); # a concise outer model
> B = [Dict{NamedTuple, JuMP.ConstraintRef}() for j = 1:J]; # a vec of Bijections::NamedTuple
> out_upd_vec = falses(J);
> parallel_CG!(B, model, θ, β); # ⚠️🔴 [deprecate] condense!(B, model)
> # Note: at this line we had derived a typically tight Lagrangian bound
> 
> # Note: The following 2 methods are both primal heuristics, hence cannot guarantee optimality, even feasibility!
> solve_restricted_primal(B; mip = true); # 1️⃣ This is a handy integer-feasible solution (typically good-quality)
> while true # 2️⃣ the diving heuristic to recover another integer-feasible solution (dive only once here)
> restrict_a_frac_block!(B, model, inn) && break
> parallel_CG!(B, model, θ, β);
> end
> 
> ```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [July 27, 2025, 3:18am UTC](https://discourse.julialang.org/t/some-question-on-column-generation-branch-and-price-coluna/131084/2 "2025-07-27T03:18:53Z")

</div>

The picture on the slide is correct. I suggest you look up some standard textbook about Benders.

You should also seach for “diving” in the context of solving a MIP and the branch and bound tree.

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [July 27, 2025, 3:55am UTC](https://discourse.julialang.org/t/some-question-on-column-generation-branch-and-price-coluna/131084/3 "2025-07-27T03:55:10Z")

</div>

> [@odow](#):
>
> The picture on the slide is correct

You’re right. I misunderstood.

So Dual decomposition and Benders decomposition has their own fields of application.

Suppose we have a large number of blocks (say, `N`), and we have a small number of non-block shared variables (say, `n`). e.g. `n = 3` while `N = 10000`. If the `n` variables are referred to in the constraints within all `N` blocks, then Benders decomposition is the fitting method. On the other hand, if a small number (`n`) of non-block constraints must be satisfied by the variables from within all `N` blocks, then Dual decomposition is the fitting method.

* * *

I noticed that in the slide above from JuMP-dev 2020, nested decomposition is mentioned. Personally I don’t think that is realistic, since the current Master-Subproblem two-layer structure is already very hard for large scale problems so that we aren’t prepared to handle the case Master-Middle-Subproblem three-layer nested formulation.

Instead, I think the following **pairwise coupling** case might worth more attention, this case is supposed to be a bit (or a lot?) harder than the _single point of coupling_ cases discussed above.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/6/e/6e93b22bc97f079859b8a8285e6a54b670c2b956.png)  
(this image itself is ready to be solved by Dual decomposition. If you lean your head towards your left shoulder by 90 degrees, then the problem corresponding to the image you see is ready to be solved by Benders decomposition. The complexity can be tuned up by lengthening the orange blocks, which indicates denser coupling.)

The reason that **pairwise coupling** is harder is:  
The coupling size is O(N), compared to O(1) in the _single point of coupling_ case, where N represents the number of blocks.