Branch and bound আর traveling salesman

১ · Backtracking-এর পরের ধাপ — “সবচেয়ে ভালো” উত্তরের খোঁজ

পর্ব ৫৬-এ backtracking দেখেছেন: সিদ্ধান্ত একটি একটি করে নিন, কোনো শর্ত ভাঙলেই ফিরে আসুন। N-queens বা subset sum-এ প্রশ্ন ছিল “একটি বৈধ উত্তর আছে কি?” কিন্তু অনেক সমস্যায় বৈধ উত্তর অগণিত, আর চাই সবচেয়ে ভালোটি — সবচেয়ে কম খরচ বা সবচেয়ে বেশি লাভ। এগুলো optimization problem.

ধরুন একজন কুরিয়ার ভ্যানচালক মিরপুরের অফিস থেকে বেরিয়ে চারটি এলাকায় পার্সেল দিয়ে আবার অফিসে ফিরবেন। প্রতিটি এলাকায় ঠিক একবার যাবেন, আর মোট সময় সবচেয়ে কম চান — এটাই traveling salesman problem (TSP). চারটি এলাকা হলে সম্ভাব্য ক্রম ৪! = ২৪টি, হাতে গোনা যায়। কিন্তু ১৯টি এলাকা হলে ১৯! ≈ ১.২ × ১০^১৭ — প্রতি সেকেন্ডে একশো কোটি ক্রম পরীক্ষা করলেও কয়েক বছর লাগবে।

সব ক্রম একে একে দেখা হলো brute force. Branch and bound-এর বুদ্ধি: অর্ধেক বানানো একটি রুট দেখে যদি প্রমাণ করা যায় “বাকিটা যেভাবেই সাজাই, মোট সময় ৯০ মিনিটের কম হবে না”, আর হাতে আগে থেকেই ৮০ মিনিটের একটি রুট থাকে — তাহলে এই অর্ধেক রুটের সব সম্প্রসারণ একসাথে বাদ।

সমাধান-গাছ আর তিন রকম node

Backtracking-এর মতোই সিদ্ধান্তগুলো একটি গাছ তৈরি করে — state-space tree. Root-এ এখনো কোনো সিদ্ধান্ত নেই; প্রতিটি স্তরে একটি সিদ্ধান্ত (“জিনিস ৩ নেব কি না”, “এরপর কোন শহর”); আর leaf মানে একটি সম্পূর্ণ সমাধান। একটি node-কে “খোলা” (expand) মানে তার সব সন্তান তৈরি করা — এটাই branch. বিশ্ববিদ্যালয়ের পাঠ্যবইয়ে (যেমন Horowitz–Sahni ধাঁচের বইয়ে) তিনটি নাম খুব চলে:

• Live node: তৈরি হয়েছে, কিন্তু এখনো খোলা হয়নি — অপেক্ষার তালিকায় আছে।
• E-node (expanding node): এই মুহূর্তে যার সন্তান তৈরি হচ্ছে।
• Dead node: খোলা শেষ, অথবা ছাঁটাই হয়ে বাদ পড়েছে।

Live node-গুলো কোন ডেটা স্ট্রাকচারে রাখবেন, সেটাই ঠিক করে খোঁজের ক্রম — queue, stack, নাকি priority queue. অংশ ৪-এ তিনটির তুলনা।

২ · Bound — আশাবাদী অনুমান আর ছাঁটাইয়ের নিয়ম

পুরো পদ্ধতিটি দাঁড়িয়ে আছে একটি সংখ্যার ওপর। ধরুন সমস্যাটি খরচ কমানোর (minimisation). একটি node x-এর bound b(x) হলো এমন একটি সংখ্যা, যাতে

অর্থাৎ b(x) একটি lower bound — x থেকে যতই ভালো সাজান, খরচ এর নিচে নামবে না। একে বলে আশাবাদী (optimistic) অনুমান: এটি সবসময় বাস্তবের চেয়ে ভালো বা সমান ছবি দেখায়, কখনো খারাপ নয়। লাভ বাড়ানোর (maximisation) সমস্যায় সবকিছু উল্টো — bound তখন upper bound, “এর বেশি লাভ অসম্ভব”.

খোঁজের মাঝপথে পাওয়া সেরা সম্পূর্ণ সমাধানের নাম incumbent, তার খরচ best. এবার ছাঁটাইয়ের (pruning) নিয়ম এক লাইনের:

একটি bound দুটি টানাপোড়েনের মাঝে থাকে। সবচেয়ে সস্তা bound “এ পর্যন্ত যা খরচ হয়েছে” — খরচ ঋণাত্মক না হলে বৈধ, কিন্তু প্রায় কিছুই ছাঁটে না। সবচেয়ে শক্ত bound x-এর নিচের আসল সর্বোত্তম খরচ — কিন্তু সেটি বের করা মানে মূল সমস্যাটাই সমাধান করা! ভালো bound মাঝখানে: দ্রুত হিসাব হয় (সাধারণত polynomial সময়ে), আর আসল মানের কাছাকাছি। প্রায় সব ভালো bound আসে relaxation থেকে: কঠিন একটি শর্ত ঢিলে করে সহজ সমস্যা বানানো, যার সর্বোত্তম মান মূল সমস্যার সর্বোত্তমের চেয়ে খারাপ হতে পারে না।

Pseudocode — best-first branch and bound (minimisation)
BRANCH-AND-BOUND(root):
    best ← cost of a quick heuristic solution   (or ∞ if none)
    Q ← priority queue keyed by bound;  INSERT(Q, root)
    while Q is not empty:
        x ← EXTRACT-MIN(Q)
        if b(x) ≥ best: stop                   // every live node is ≥ best
        for each child y of x:                 // branch
            if y is a complete solution:
                if cost(y) < best: best ← cost(y); remember y
            else if b(y) < best: INSERT(Q, y)  // bound: keep y alive
            // otherwise y is pruned
    return best and the remembered solution

Q-কে সাধারণ queue বা stack বানালে বাকি সবকিছু একই থাকে, শুধু “stop”-এর জায়গায় “এই x বাদ দিয়ে পরেরটি নাও” (continue) লিখতে হয় — কারণ তখন Q-এর সামনের node-টিই সবচেয়ে ছোট bound-এর হবে, এমন কোনো নিশ্চয়তা নেই।

কেন উত্তর ঠিক — একটি invariant

শুরুতে root live, আর সব সমাধানই root-এর নিচে — সত্য। একটি node খুললে x-এর জায়গায় তার সন্তানেরা আসে; x-এর subtree-র প্রতিটি সমাধান ঠিক একটি সন্তানের subtree-তে, তাই কিছু হারায় না। y ছাঁটাই হলে b(y) ≥ best; y-এর নিচের প্রতিটি সমাধানের খরচ ≥ b(y) ≥ best, তাই সেখানে incumbent-এর চেয়ে কম খরচের কিছু নেই — যদি সেখানে কোনো সর্বোত্তম থাকত, তার খরচ best-এর সমান হতো, অর্থাৎ incumbent নিজেই সর্বোত্তম। দুই ক্ষেত্রেই invariant টেকে। শেষে Q খালি, অথবা সবচেয়ে ছোট bound-ও ≥ best — কোনো live node-এর নিচে আরও ভালো কিছু নেই, তাই incumbent সর্বোত্তম। গাছ সসীম, প্রতিটি node বড়জোর একবার খোলা হয়, তাই লুপ থামবেই। ∎

খেয়াল করুন, প্রমাণে লাগেনি যে bound “ভালো” — শুধু লেগেছে যে এটি বৈধ। দুর্বল bound উত্তর ভুল করে না, শুধু ধীর করে। কিন্তু অবৈধ bound — যা কোনো সমাধানের আসল খরচের চেয়ে বেশি বলে — নিঃশব্দে সর্বোত্তম সমাধান ছেঁটে ফেলতে পারে। Branch and bound-এর প্রায় সব বাগ এখানেই।

৩ · 0/1 knapsack — fractional knapsack দিয়ে upper bound

প্রথম পূর্ণ উদাহরণ 0/1 knapsack. একজন হকারের ঝুড়িতে ১০ কেজি ধরে। চারটি জিনিস, প্রতিটি হয় পুরো নেবেন, নয় নেবেন না:

জিনিসগুলো আগেই মূল্য/ওজন অনুপাতের বড় থেকে ছোট ক্রমে সাজানো — এটা জরুরি, একটু পরেই বুঝবেন। গাছের স্তর k-তে সিদ্ধান্ত: জিনিস k নেব (১) নাকি নেব না (০)। একটি node-এর নাম তার সিদ্ধান্তের সারি — যেমন [10] মানে “জিনিস ১ নিয়েছি, ২ নিইনি, বাকিগুলো এখনো ঠিক হয়নি”.

Bound: ভগ্নাংশ নেওয়া গেলে কত পেতাম?

এটি লাভ বাড়ানোর সমস্যা, তাই চাই upper bound. Relaxation: “পুরো নাও বা নিও না” শর্তটি ঢিলে করে বলি, জিনিসের ভগ্নাংশও নেওয়া যায় — চাল বা ডালের মতো, যতটুকু জায়গা আছে ততটুকু। এটাই fractional knapsack, আর পর্ব ৫০-এ দেখেছেন এটি greedy দিয়েই নিখুঁত সমাধান হয়: অনুপাতের ক্রমে পুরো জিনিস নিতে থাকুন, যেটি আর পুরো ধরে না তার যতটুকু ধরে ততটুকু নিন, ব্যস।

কেন এটি বৈধ upper bound: node x-এর নিচের প্রতিটি 0/1 সমাধান একই সাথে x-এর সিদ্ধান্ত মেনে চলা একটি fractional সমাধানও (০ আর ১-ও তো ভগ্নাংশ)। Fractional সমস্যার সর্বোত্তম মান তাই প্রতিটি 0/1 সমাধানের মানের ≥. আর greedy সেই সর্বোত্তম মানটিই দেয় (পর্ব ৫০-এর exchange argument). মূল্যগুলো পূর্ণসংখ্যা, তাই 0/1 উত্তরও পূর্ণসংখ্যা — bound-কে নিচের দিকে পূর্ণসংখ্যায় নামালেও (floor) বৈধ থাকে, আর একটু শক্ত হয়।

Root-এ হিসাব: জিনিস ১ পুরো (৪ কেজি, ৪০), বাকি জায়গা ৬ কেজি; জিনিস ২-এর ওজন ৭, পুরো ধরে না, তাই তার ৬/৭ অংশ: ৬ × ৪২/৭ = ৩৬. Upper bound = ৪০ + ৩৬ = ৭৬. কোনো 0/1 বাছাই ৭৬-এর বেশি দিতে পারবে না। এবার best-first খোঁজ: প্রতিবার সবচেয়ে বড় upper bound-এর live node খুলুন।

কোডে — best-first knapsack, DP দিয়ে যাচাই

প্রোগ্রামটি ঠিক ওপরের খোঁজটি চালায়, প্রতিটি সিদ্ধান্ত ছাপে, আর শেষে পর্ব ৫৩-এর knapsack DP দিয়ে উত্তর মেলায়। Priority queue-তে max-heap দরকার; C++ আর Python দুই জায়গাতেই min-heap-এ -ub রেখে কাজটি সারা হয়েছে, যাতে দুটির ক্রম হুবহু এক থাকে।

#include <algorithm>
#include <cstdio>
#include <initializer_list>
#include <queue>
#include <string>
#include <tuple>
#include <vector>
using namespace std;
const int W = 10;
const vector<pair<int, int>> items = {{4, 40}, {7, 42}, {5, 25}, {3, 12}};  // (weight, value)
int bound(size_t k, int w, int v) {            // fractional knapsack on items k.. (upper bound)
    for (; k < items.size(); k++) {
        auto [wi, vi] = items[k];
        if (w + wi > W) return v + (W - w) * vi / wi;   // floor: still valid for integers
        w += wi; v += vi; }
    return v;
}
int main() {
    int best = 0, expanded = 0, pruned = 0; string bestPath;
    using Node = tuple<int, string, int, int>;            // (-ub, path, weight, value)
    priority_queue<Node, vector<Node>, greater<Node>> pq;  // smallest -ub = largest ub first
    pq.push({-bound(0, 0, 0), "", 0, 0});
    while (!pq.empty()) {
        auto [nub, path, w, v] = pq.top(); pq.pop();
        if (-nub <= best) { printf("pop [%s] ub=%d <= best %d: pruned\n", path.c_str(), -nub, best); pruned++; continue; }
        expanded++;
        printf("expand [%s] w=%d v=%d ub=%d\n", path.c_str(), w, v, -nub);
        auto [wi, vi] = items[path.size()];
        for (auto [bit, nw, nv] : {tuple{'1', w + wi, v + vi}, tuple{'0', w, v}}) {
            string p = path + bit;
            if (nw > W) { printf("  [%s] infeasible\n", p.c_str()); continue; }
            int ub = bound(p.size(), nw, nv);
            if (nv > best) { best = nv; bestPath = p; }
            if (ub <= best) { printf("  [%s] v=%d ub=%d pruned\n", p.c_str(), nv, ub); pruned++; }
            else { printf("  [%s] v=%d ub=%d queued\n", p.c_str(), nv, ub); pq.push({-ub, p, nw, nv}); }
        }
    }
    vector<int> dp(W + 1, 0);                              // 0/1 knapsack DP to check
    for (auto [wi, vi] : items) for (int c = W; c >= wi; c--) dp[c] = max(dp[c], dp[c - wi] + vi);
    printf("B&B best = %d take [%s] | expanded %d, pruned %d\n", best, bestPath.c_str(), expanded, pruned);
    printf("DP best = %d | match: %s\n", dp[W], dp[W] == best ? "yes" : "no");
}

expand [] w=0 v=0 ub=76
  [1] v=40 ub=76 queued
  [0] v=0 ub=57 queued
expand [1] w=4 v=40 ub=76
  [11] infeasible
  [10] v=40 ub=69 queued
expand [10] w=4 v=40 ub=69
  [101] v=65 ub=69 queued
  [100] v=40 ub=52 pruned
expand [101] w=9 v=65 ub=69
  [1011] infeasible
  [1010] v=65 ub=65 pruned
pop [0] ub=57 <= best 65: pruned
B&B best = 65 take [101] | expanded 4, pruned 3
DP best = 65 | match: yes

import heapq
W = 10
items = [(4, 40), (7, 42), (5, 25), (3, 12)]    # (weight, value), best value/weight first

def bound(k, w, v):                 # fractional knapsack on items k.. (an upper bound)
    for wi, vi in items[k:]:
        if w + wi > W: return v + (W - w) * vi // wi   # floor: still valid for integers
        w, v = w + wi, v + vi
    return v

best, best_path, expanded, pruned = 0, "", 0, 0
heap = [(-bound(0, 0, 0), "", 0, 0)]             # (-ub, path, weight, value): max ub first
while heap:
    nub, path, w, v = heapq.heappop(heap)
    if -nub <= best:                              # cannot beat the incumbent
        print(f"pop [{path}] ub={-nub} <= best {best}: pruned"); pruned += 1; continue
    expanded += 1
    print(f"expand [{path}] w={w} v={v} ub={-nub}")
    wi, vi = items[len(path)]
    for bit, nw, nv in (("1", w + wi, v + vi), ("0", w, v)):
        if nw > W: print(f"  [{path + bit}] infeasible"); continue
        ub = bound(len(path) + 1, nw, nv)
        if nv > best: best, best_path = nv, path + bit
        if ub <= best: print(f"  [{path + bit}] v={nv} ub={ub} pruned"); pruned += 1
        else: print(f"  [{path + bit}] v={nv} ub={ub} queued"); heapq.heappush(heap, (-ub, path + bit, nw, nv))

dp = [0] * (W + 1)                                  # 0/1 knapsack DP to check
for wi, vi in items:
    for c in range(W, wi - 1, -1): dp[c] = max(dp[c], dp[c - wi] + vi)
print(f"B&B best = {best} take [{best_path}] | expanded {expanded}, pruned {pruned}")
print(f"DP best = {dp[W]} | match: {'yes' if dp[W] == best else 'no'}")

expand [] w=0 v=0 ub=76
  [1] v=40 ub=76 queued
  [0] v=0 ub=57 queued
expand [1] w=4 v=40 ub=76
  [11] infeasible
  [10] v=40 ub=69 queued
expand [10] w=4 v=40 ub=69
  [101] v=65 ub=69 queued
  [100] v=40 ub=52 pruned
expand [101] w=9 v=65 ub=69
  [1011] infeasible
  [1010] v=65 ub=65 pruned
pop [0] ub=57 <= best 65: pruned
B&B best = 65 take [101] | expanded 4, pruned 3
DP best = 65 | match: yes

জটিলতা: গাছে বড়জোর 2ⁿ⁺¹ − 1 টি node, প্রতিটির bound O(n) (prefix sum আর binary search দিয়ে O(log n)), heap-এর প্রতিটি কাজ O(log 2ⁿ) = O(n). তাই worst case সময় O(n · 2ⁿ), শুরুর সর্ট O(n log n). Best-first-এর queue-তে worst case-এ O(2ⁿ) টি node, প্রতিটিতে O(n) লম্বা path — জায়গা O(n · 2ⁿ); depth-first-এ মাত্র O(n). বাস্তবে fractional bound খুব শক্ত, তবে নিশ্চয়তা নেই — knapsack NP-hard.

৪ · খোঁজের তিন কৌশল — FIFO, depth-first আর best-first

একই bound, একই ছাঁটাইয়ের নিয়ম — শুধু live node রাখার পাত্র বদলালে খোঁজের চেহারা পুরো বদলে যায়। ট্রেনের টিকিটের লাইনের কথা ভাবুন: যে আগে এসেছে সে আগে (queue); যে সবশেষে এসেছে সে আগে (stack); নাকি যার প্রয়োজন সবচেয়ে জরুরি সে আগে (হাসপাতালের জরুরি বিভাগের মতো — priority queue).

অংশ ৫-এর পাঁচ শহরের TSP-তে (reduced cost matrix bound) তিনটি কৌশল চালিয়ে গোনা হয়েছে — কতটি node খোলা হলো, কতটি সন্তান তৈরি হলো, আর queue/stack সবচেয়ে বড় কত হলো। দুইভাবে: শুরুতে কোনো incumbent ছাড়া (best = ∞), আর nearest neighbour-এর দেওয়া ৩৬ খরচের রুট দিয়ে শুরু করে (অংশ ৯)। সম্পূর্ণ tour তৈরি হওয়ামাত্র incumbent আপডেট হয়।

তিনটি শিক্ষা। এক, best-first সবচেয়ে কম node খোলে। দুই, FIFO-র সবচেয়ে বড় শত্রু incumbent-এর অভাব — ভালো শুরুর উত্তরেই ২০ থেকে ৬. তিন, best-first-এ শুরুর উত্তর খোলা node কমায়নি, কিন্তু queue ৭ থেকে ৩ — খারাপ সন্তান তৈরি হওয়ামাত্র ছাঁটা পড়ে। বড় সমস্যায় এই মেমরির পার্থক্যই ঠিক করে প্রোগ্রাম শেষ হবে, নাকি মেমরি ফুরিয়ে থেমে যাবে।

Best-first কেন সবচেয়ে কম node খোলে

ধরুন bound monotone: সন্তানের bound কখনো মা-বাবার চেয়ে কম নয় (এই পর্বের সব bound-ই এমন)। সর্বোত্তম খরচ OPT. দাবি ১: b(x) < OPT এমন প্রতিটি node যেকোনো সঠিক branch and bound-কে খুলতেই হবে, যদি তার পূর্বপুরুষেরা খোলা হয় — কারণ best কখনো OPT-এর নিচে নামে না, তাই b(x) < OPT ≤ best, ছাঁটাইয়ের শর্ত খাটে না। দাবি ২: best-first কখনো b(x) > OPT এমন node খোলে না — সর্বোত্তম সমাধানের পথের প্রতিটি node-এর bound ≤ OPT, তাদের কেউ না কেউ সবসময় queue-তে থাকে আর আগে বের হয়, আর সর্বোত্তম leaf পাওয়ামাত্র best = OPT হয়ে যায়। ফলে best-first শুধু bound < OPT-এর node (আর কিছু bound = OPT-এর node) খোলে — যা না খুলে উপায় নেই। পর্ব ৪৩-এর A*-এর সর্বোত্তমতার যুক্তিও হুবহু এটাই: A* আসলে shortest path-এর ওপর best-first branch and bound, যেখানে admissible heuristic হলো lower bound.

৫ · TSP-এর reduced cost matrix — পাঠ্যবইয়ের ক্লাসিক bound

এবার কুরিয়ার ভ্যানের সমস্যায় ফিরি। পাঁচটি জায়গা, ১ হলো অফিস। Cᵢⱼ = i থেকে j যেতে মিনিট — ঢাকার একমুখী রাস্তা আর জ্যামের কারণে i → j আর j → i সমান নয় (asymmetric TSP). নিজের কাছে যাওয়া নিষেধ, তাই কর্ণে ∞. পর্ব ৫৫-এ এই সমস্যার Held–Karp DP দেখেছেন — O(n² 2ⁿ) সময়, কিন্তু O(n 2ⁿ) মেমরিও লাগে; পর্ব ৫৬-এ Hamiltonian cycle-এর backtracking. এখন branch and bound.

কেন সারি-কলাম বিয়োগ একটি lower bound

মূল পর্যবেক্ষণ: প্রতিটি tour প্রতিটি শহর থেকে ঠিক একবার বের হয়, আর প্রতিটি শহরে ঠিক একবার ঢোকে। অর্থাৎ ম্যাট্রিক্সের প্রতিটি সারি থেকে ঠিক একটি আর প্রতিটি কলাম থেকে ঠিক একটি ঘর ব্যবহার করে। এখন সারি i-এর প্রতিটি ঘর থেকে একই সংখ্যা r বিয়োগ করলে প্রতিটি tour-এর খরচ ঠিক r কমে — কারণ প্রতিটি tour সারি i-এর ঠিক একটি ঘর নেয়। কলামের বেলায়ও তাই। সারির minimum বিয়োগ করলে কোনো ঘর ঋণাত্মক হয় না।

তাই সব সারি আর কলাম থেকে মোট R বিয়োগ করে পাওয়া ম্যাট্রিক্সে (reduced cost matrix) প্রতিটি tour-এর খরচ = আসল খরচ − R, আর এটি ≥ ০ (সব ঘর অঋণাত্মক)। সুতরাং আসল খরচ ≥ R — প্রতিটি tour-এর জন্য। R-ই root-এর lower bound. নিচের স্টেপারে পুরো হিসাব, তারপর root-এর চারটি সন্তানের bound.

সন্তানের নিয়মটি কেন বৈধ: মা-বাবার reduced ম্যাট্রিক্স M-এ প্রতিটি tour-এর খরচ = মা-বাবার bound + M-এ tour-এর ঘরগুলোর যোগফল (ওপরের যুক্তিতে)। i → j ধরে নিলে tour-এর একটি ঘর M[i][j]; বাকি ঘরগুলো সারি i বা কলাম j-তে থাকতে পারে না, আর (j, ১)-ও নয় (শেষ ধাপ ছাড়া)। তাই ওগুলো ∞ করলে কোনো বৈধ tour হারায় না, আর নতুন ম্যাট্রিক্সকে আবার reduce করলে আরও r নিশ্চিত হয়। মোট: b(সন্তান) = b(মা-বাবা) + M[i][j] + r ≥ b(মা-বাবা) — monotone. Little, Murty, Sweeney ও Karel-এর ১৯৬৩ সালের মূল পদ্ধতি একটু অন্যভাবে branch করত — প্রতিটি ধাপে একটি edge “নাও” আর “নিও না” এই দুই ভাগে — কিন্তু bound-এর চিন্তাটি এই সারি-কলাম বিয়োগই।

পুরো খোঁজ — best-first, nearest neighbour দিয়ে শুরু

শুরুর incumbent হিসেবে একটি দ্রুত আনুমানিক রুট নিই: অফিস থেকে প্রতিবার সবচেয়ে কাছের না-যাওয়া জায়গায় যাওয়া — 1 → 4 (১০) → 5 (৩) → 2 (৪) → 3 (১৬) → 1 (৩), মোট ৩৬ (অংশ ৯-এ বিস্তারিত)। এবার গাছ:

জটিলতা: গভীরতা k-তে node সংখ্যা (n−1)(n−2)…(n−k), মোট ≤ e · (n−1)!. প্রতিটি node-এ ম্যাট্রিক্স কপি আর reduce O(n²), heap-এর কাজ O(log((n−1)!)) = O(n log n). তাই worst case সময় O(n² · (n−1)!) = O(n · n!) — brute force-এর মতোই খারাপ; উন্নতি শুধু bound যতটা ছাঁটে। জায়গা: best-first-এ প্রতিটি live node নিজের n × n ম্যাট্রিক্স রাখে, worst case O(n² · (n−1)!); depth-first-এ শুধু পথের n টি ম্যাট্রিক্স, O(n³).

৬ · কোডে TSP branch and bound — brute force দিয়ে যাচাই

প্রোগ্রামটি প্রথমে nearest neighbour দিয়ে incumbent বানায়, তারপর ওপরের best-first খোঁজ চালায়: priority queue-তে (bound, পথ, reduced ম্যাট্রিক্স)। শেষে ২৪টি tour-এর brute force দিয়ে উত্তর মেলায় — ছোট ইনপুটে এমন যাচাই প্রতিটি branch and bound কোডের সাথে থাকা উচিত।

#include <algorithm>
#include <cstdio>
#include <queue>
#include <string>
#include <tuple>
#include <vector>
using namespace std;
using Mat = vector<vector<int>>; const int INF = 1e9, n = 5;
const Mat C = {{INF, 20, 30, 10, 11}, {15, INF, 16, 4, 2}, {3, 5, INF, 2, 4}, {19, 6, 18, INF, 3}, {16, 4, 7, 16, INF}};
int cost(const vector<int>& t) { int s = 0; for (int k = 0; k < n; k++) s += C[t[k]][t[(k + 1) % n]]; return s; }
string fmt(const vector<int>& p) { string s; for (int c : p) s += (s.empty() ? "" : "-") + to_string(c + 1); return s; }
int reduce(Mat& M) {                              // subtract row minima, then column minima
    int total = 0;
    for (auto& row : M) { int m = ranges::min(row); if (m > 0 && m < INF) { total += m; for (int& x : row) if (x < INF) x -= m; } }
    for (int j = 0; j < n; j++) { int m = INF; for (auto& row : M) m = min(m, row[j]);
        if (m > 0 && m < INF) { total += m; for (auto& row : M) if (row[j] < INF) row[j] -= m; } }
    return total;
}
int main() {
    vector<int> t = {0};                          // nearest neighbour gives the first incumbent
    while ((int)t.size() < n) { int nx = -1;
        for (int j = 0; j < n; j++) if (ranges::find(t, j) == t.end() && (nx < 0 || C[t.back()][j] < C[t.back()][nx])) nx = j;
        t.push_back(nx); }
    int best = cost(t), expanded = 0; vector<int> tour = t;
    printf("nearest neighbour %s-1 cost %d\n", fmt(t).c_str(), best);
    using Node = tuple<int, vector<int>, Mat>;    // (bound, path, reduced matrix), min-heap
    priority_queue<Node, vector<Node>, greater<Node>> pq; Mat M0 = C; pq.push({reduce(M0), {0}, M0});
    while (!pq.empty() && get<0>(pq.top()) < best) {   // smallest bound >= best: stop
        auto [b, path, M] = pq.top(); pq.pop();
        expanded++; printf("expand %s bound %d\n", fmt(path).c_str(), b);
        for (int i = path.back(), j = 0; j < n; j++) if (M[i][j] < INF) {
            Mat M2 = M; vector<int> p = path; p.push_back(j);
            for (int k = 0; k < n; k++) M2[i][k] = M2[k][j] = INF;
            if ((int)p.size() < n) M2[j][0] = INF;        // no early return to city 1
            int nb = b + M[i][j] + reduce(M2);
            if (nb >= best) printf("  %s bound %d pruned\n", fmt(p).c_str(), nb);
            else if ((int)p.size() == n) { best = nb; tour = p; printf("  %s-1 complete, new best %d\n", fmt(p).c_str(), nb); }
            else { printf("  %s bound %d queued\n", fmt(p).c_str(), nb); pq.push({nb, p, M2}); }
        }
    }
    printf("best %s-1 cost %d | expanded %d, left in queue %zu\n", fmt(tour).c_str(), best, expanded, pq.size());
    vector<int> perm = {0, 1, 2, 3, 4}; int bf = INF;   // brute force: city 1 fixed, 4! orders
    do bf = min(bf, cost(perm)); while (next_permutation(perm.begin() + 1, perm.end()));
    printf("brute force over all 24 tours: %d | match: %s\n", bf, bf == best ? "yes" : "no");
}

nearest neighbour 1-4-5-2-3-1 cost 36
expand 1 bound 25
  1-2 bound 35 queued
  1-3 bound 53 pruned
  1-4 bound 25 queued
  1-5 bound 31 queued
expand 1-4 bound 25
  1-4-2 bound 28 queued
  1-4-3 bound 50 pruned
  1-4-5 bound 36 pruned
expand 1-4-2 bound 28
  1-4-2-3 bound 52 pruned
  1-4-2-5 bound 28 queued
expand 1-4-2-5 bound 28
  1-4-2-5-3-1 complete, new best 28
best 1-4-2-5-3-1 cost 28 | expanded 4, left in queue 2
brute force over all 24 tours: 28 | match: yes

import heapq, itertools
INF = 10**9
C = [[INF, 20, 30, 10, 11], [15, INF, 16, 4, 2], [3, 5, INF, 2, 4],
     [19, 6, 18, INF, 3], [16, 4, 7, 16, INF]]
n = len(C)
cost = lambda t: sum(C[t[k]][t[(k + 1) % n]] for k in range(n))
fmt = lambda p: "-".join(str(c + 1) for c in p)

def reduce(M):                                  # subtract row minima, then column minima
    total = 0
    for row in M:
        m = min(row)
        if 0 < m < INF: total += m; row[:] = [x - m if x < INF else INF for x in row]
    for j in range(n):
        m = min(row[j] for row in M)
        if 0 < m < INF:
            total += m
            for row in M: row[j] = row[j] - m if row[j] < INF else INF
    return total

t = [0]                                         # nearest neighbour gives the first incumbent
while len(t) < n: t.append(min((j for j in range(n) if j not in t), key=lambda j: C[t[-1]][j]))
best, tour = cost(t), t
print(f"nearest neighbour {fmt(t)}-1 cost {best}")
M = [row[:] for row in C]
heap, expanded = [(reduce(M), [0], M)], 0      # (bound, path, reduced matrix)
while heap and heap[0][0] < best:               # smallest bound >= best: stop
    b, path, M = heapq.heappop(heap)
    expanded += 1; print(f"expand {fmt(path)} bound {b}")
    i = path[-1]
    for j in (j for j in range(n) if M[i][j] < INF):
        M2 = [row[:] for row in M]
        for k in range(n): M2[i][k] = M2[k][j] = INF
        if len(path) + 1 < n: M2[j][0] = INF   # no early return to city 1
        nb, p = b + M[i][j] + reduce(M2), path + [j]
        if nb >= best: print(f"  {fmt(p)} bound {nb} pruned")
        elif len(p) == n: best, tour = nb, p; print(f"  {fmt(p)}-1 complete, new best {nb}")
        else: print(f"  {fmt(p)} bound {nb} queued"); heapq.heappush(heap, (nb, p, M2))
print(f"best {fmt(tour)}-1 cost {best} | expanded {expanded}, left in queue {len(heap)}")
bf = min(cost((0,) + p) for p in itertools.permutations(range(1, n)))
print(f"brute force over all 24 tours: {bf} | match: {'yes' if bf == best else 'no'}")

nearest neighbour 1-4-5-2-3-1 cost 36
expand 1 bound 25
  1-2 bound 35 queued
  1-3 bound 53 pruned
  1-4 bound 25 queued
  1-5 bound 31 queued
expand 1-4 bound 25
  1-4-2 bound 28 queued
  1-4-3 bound 50 pruned
  1-4-5 bound 36 pruned
expand 1-4-2 bound 28
  1-4-2-3 bound 52 pruned
  1-4-2-5 bound 28 queued
expand 1-4-2-5 bound 28
  1-4-2-5-3-1 complete, new best 28
best 1-4-2-5-3-1 cost 28 | expanded 4, left in queue 2
brute force over all 24 tours: 28 | match: yes

৭ · আরও bound — দুই সস্তা edge, MST আর 1-tree

Reduced matrix ছাড়াও TSP-এর জন্য বেশ কিছু bound আছে, বিশেষত symmetric TSP-এ (i → j আর j → i-এর খরচ সমান — যেমন সোজা দূরত্ব)। একটি পাঁচ শহরের উদাহরণ নিই — a, b, c, d, e:

সর্বোত্তম tour a–b–d–e–c–a: ৩ + ৭ + ৩ + ২ + ১ = ১৬ (১২টি আলাদা tour-এর brute force দিয়ে যাচাই করা)। এবার চারটি bound তুলনা করি।

দুই সস্তা edge-এর যোগফলের অর্ধেক

Tour-এ প্রতিটি শহর ঠিক দুটি edge-এ লাগানো — একটিতে ঢোকে, একটিতে বের হয়। প্রতিটি শহরের জন্য তার দুটি tour-edge যোগ করলে প্রতিটি edge দুবার গোনা হয় (দুই প্রান্তে), তাই যোগফল = ২ × tour. আবার প্রতিটি শহরের দুটি tour-edge-এর যোগফল ≥ সেই শহরের দুটি সবচেয়ে সস্তা edge-এর যোগফল। তাই

এখানে: a (১, ৩), b (৩, ৬), c (১, ২), d (৩, ৪), e (২, ৩) — মোট ২৮, bound ১৪. দূরত্ব পূর্ণসংখ্যা হলে tour-ও পূর্ণসংখ্যা, তাই ceil বৈধ। Branch করার সময় পথে নেওয়া edge-গুলোকে সংশ্লিষ্ট শহরের “দুই edge”-এর মধ্যে জোর করে ধরতে হয়, তাতে গাছের যত গভীরে যাবেন, bound তত শক্ত হয়। আগে থেকে প্রতিটি শহরের দুই minimum রাখলে হিসাব O(n).

MST আর 1-tree

পর্ব ৪৫-এর minimum spanning tree মনে করুন। একটি tour থেকে যেকোনো একটি edge সরালে থাকে সব শহর ছোঁয়া একটি পথ — সেটি একটি spanning tree. তাই tour > tour − (একটি edge) ≥ MST. এখানে MST = a–c ১, c–e ২, a–b ৩, d–e ৩ = ৯ — বৈধ, কিন্তু দুর্বল।

উন্নতি — 1-tree: একটি শহর (ধরুন a) আলাদা করে রাখুন। Tour থেকে a আর তার দুটি edge সরালে বাকি চার শহরে একটি পথ থাকে, যা ওই চার শহরের একটি spanning tree. তাই tour = (a-র দুটি edge) + (বাকিদের একটি spanning tree) ≥ (a-র দুটি সস্তা edge) + (বাকিদের MST). এখানে {b, c, d, e}-এর MST = c–e ২, d–e ৩, b–c ৬ = ১১, আর a-র দুই সস্তা edge ১ + ৩ = ৪; 1-tree bound = ১৫. কোন শহরকে আলাদা রাখবেন তা বদলে সবচেয়ে বড়টি নিতে পারেন — এখানে a বা b নিলে ১৫, অন্যগুলোয় ১৩-১৪.

৮ · Job assignment — প্রতিটি সারির সবচেয়ে সস্তা কাজ

একটি গার্মেন্টস লাইনে চারজন অপারেটর — a, b, c, d — আর চারটি কাজ: ১ কাটিং, ২ সেলাই, ৩ বোতাম, ৪ প্যাকিং. খরচের ম্যাট্রিক্সে সারি = অপারেটর, কলাম = কাজ, মান = মিনিট। প্রত্যেককে ঠিক একটি কাজ, প্রতিটি কাজ ঠিক একজনকে; মোট মিনিট সবচেয়ে কম চাই — এটাই assignment problem.

Bound: যাঁদের কাজ ঠিক হয়েছে তাঁদের খরচ + বাকি প্রত্যেক অপারেটরের, এখনো খালি কাজগুলোর মধ্যে সবচেয়ে সস্তাটি। বৈধ কেন: যেকোনো পূর্ণ বণ্টনে বাকি প্রত্যেক অপারেটর কোনো একটি খালি কাজ পাবেনই, আর তার খরচ সেই সারির minimum-এর ≥. শিথিল করা শর্ত হলো “দুজন একই কাজ পাবেন না” — bound দুজনকে একই সস্তা কাজে বসাতে দেয়। প্রতিটি স্তরে একজন অপারেটরের কাজ ঠিক হয়: স্তর ১-এ a, স্তর ২-এ b, এভাবে।

৯ · ভালো শুরু, বাস্তবের solver আর একটু ইতিহাস

শুরুর উত্তর কেন এত জরুরি

ছাঁটাইয়ের নিয়মের দুই দিক: bound যত শক্ত, আর incumbent যত ভালো, তত বেশি ছাঁটা যায়। দ্বিতীয়টি সস্তায় পাওয়ার উপায় একটি দ্রুত heuristic. TSP-এর সবচেয়ে সরলটি nearest neighbour: যেখানে আছেন সেখান থেকে সবচেয়ে কাছের না-যাওয়া শহরে যান, O(n²) সময়ে। আমাদের ম্যাট্রিক্সে এটি দিয়েছে ৩৬, সর্বোত্তমের চেয়ে প্রায় ২৯% বেশি — তবু অংশ ৪-এর টেবিলে দেখেছেন, এটুকুতেই FIFO-র খোলা node ২০ থেকে ৬, আর সর্বোচ্চ queue ২৪ থেকে ৩.

Nearest neighbour-এর গুণগত নিশ্চয়তা নেই — triangle inequality মানা ইনপুটেও এটি সর্বোত্তমের log n-এর ক্রমের গুণ খারাপ হতে পারে (Rosenkrantz, Stearns ও Lewis, ১৯৭৭)। তাই এর পরে চালানো হয় local search, যেমন 2-opt: দুটি edge সরিয়ে উল্টো করে জোড়া দিলে tour ছোট হলে বদলান। কয়েক মিলিসেকেন্ডের এই উন্নতি branch and bound-এর ঘণ্টা বাঁচাতে পারে।

বাস্তবের solver — LP relaxation আর branch and cut

শিল্পে branch and bound-এর সবচেয়ে বড় ব্যবহার integer programming-এ: এয়ারলাইনের ক্রু-রোস্টার, বাস-ট্রেনের সময়সূচি, কারখানার উৎপাদন পরিকল্পনা, বিদ্যুৎকেন্দ্র কখন চালু হবে। সমস্যাগুলো লেখা হয় linear শর্ত আর পূর্ণসংখ্যা চলক দিয়ে। Bound আসে LP relaxation থেকে: “পূর্ণসংখ্যা হতে হবে” শর্ত তুলে দিয়ে সাধারণ linear program, যা simplex বা interior-point পদ্ধতিতে দ্রুত সমাধান হয়। LP-র উত্তরে কোনো চলক x = 2.6 এলে branch: একটি শাখায় x ≤ 2, অন্যটিতে x ≥ 3 — Land ও Doig-এর মূল পদ্ধতির আধুনিক রূপ। (Fractional knapsack আসলে knapsack-এর LP relaxation-ই।)

আধুনিক solver (Gurobi, CPLEX, SCIP, HiGHS ইত্যাদি) প্রতিটি node-এ LP-তে নতুন বৈধ শর্ত — cutting plane — যোগ করে, যা প্রতিটি পূর্ণসংখ্যা সমাধান মানে কিন্তু বর্তমান ভগ্নাংশ উত্তরটিকে কেটে বাদ দেয়; bound শক্ত হয়। এই মিশ্রণের নাম branch and cut. Solver চলার সময় incumbent আর সেরা bound-এর ফাঁক (gap) দেখায়, আর “gap ১%-এর নিচে নামলেই থামো” বলা যায় — প্রমাণিত মানের উত্তর, অনেক কম সময়ে।

TSP-এর চূড়ান্ত উদাহরণ Concorde — David Applegate, Robert Bixby, Vašek Chvátal ও William Cook-এর branch-and-cut প্রোগ্রাম। এটি ২০০৪ সালে সুইডেনের ২৪,৯৭৮টি শহরের আর ২০০৬ সালে একটি চিপ-নকশা থেকে আসা ৮৫,৯০০ বিন্দুর TSP প্রমাণসহ সর্বোত্তম সমাধান করেছে। সম্ভব হয়েছে শক্ত bound (Held–Karp-ধাঁচের relaxation আর হাজারো cut), ভালো শুরুর tour (Lin–Kernighan local search) আর চতুর branching-এর মিলনে।

১০ · AI-এর যুগে branch and bound — bound-এর বৈধতাই আসল প্রশ্ন

“TSP-এর branch and bound লিখে দাও” — AI চটপট কোড দেবে, আর ৪-৫ শহরের উদাহরণে উত্তরও প্রায়ই ঠিক আসবে। কিন্তু branch and bound-এর বাগ ভুল উত্তর দেয় মাঝেমধ্যে: একটি অবৈধ bound বেশিরভাগ ইনপুটে ক্ষতি করে না, কোনো এক ইনপুটে সর্বোত্তম শাখাটি ছেঁটে দেয়। আর জটিলতার দাবিতে AI প্রায়ই worst case আড়াল করে।

এই বিষয়ে AI সাধারণত যেসব ভুল করে

Branch and bound কোড যাচাইয়ের চেকলিস্ট

• Bound-এর বৈধতার একটি লিখিত যুক্তি আছে (relaxation কোনটি, কেন সেটি প্রতিটি সমাধানের চেয়ে আশাবাদী)
• একটি ছোট ইনপুটে গাছের প্রথম ৩-৪টি E-node হাতে trace করে কোডের ছাপা লাইনের সাথে মিলিয়েছি
• n ≤ 8-এর শত শত এলোমেলো ইনপুটে brute force বা DP-র সাথে উত্তর মিলিয়েছি — stress testing
• Edge case: n = 1 বা ২, সব খরচ সমান, ধারণক্ষমতা ০, কোনো জিনিসই ধরে না, একাধিক সর্বোত্তম উত্তর
• Minimisation না maximisation — ছাঁটাইয়ের অসমতার দিক আর heap-এর দিক (min বা max) দুটোই মিলিয়ে দেখেছি
• AI জটিলতা বললে জিজ্ঞেস করেছি “worst case-এ কতটি node?” — উত্তর exponential না হলে সন্দেহ করেছি