Cấu trúc dữ liệu và Giải thuật Bài 131: Thuật toán Hungarian (Kuhn-Munkres / KM Algorithm) - Ghép cặp cực đại với trọng số cực đại trên đồ thị hai phía trong $O(V^3)$.

Trong lý thuyết đồ thị và tối ưu hóa tổ hợp, Bài toán phân công (Assignment Problem) là một bài toán kinh điển: Có $N$ công nhân và $N$ công việc, mỗi công nhân thực hiện một công việc cụ thể sẽ tạo ra giá trị lợi nhuận $w(i, j)$ (hoặc chi phí tương ứng). Cần ghép mỗi công nhân với đúng một công việc sao cho tổng lợi nhuận đạt giá trị lớn nhất.

Dưới góc nhìn đồ thị, đây chính là bài toán Ghép cặp hoàn hảo có trọng số cực đại trên đồ thị hai phía (Maximum Weight Perfect Bipartite Matching).

Nếu mô hình hóa bài toán này thành mạng luồng chi phí cực tiểu (Min-Cost Max-Flow - Bài 118), thuật toán SPFA hoặc Dijkstra kết hợp thế năng Johnson giải quyết được trong thời gian $O(V^2 E) = O(V^4)$ hoặc $O(V \cdot E \log V)$. Harold Kuhn (1955) và James Munkres (1957) đã phát triển Thuật toán Kuhn-Munkres (KM Algorithm / Hungarian Algorithm), dựa trên nền tảng phương pháp đối ngẫu (Primal-Dual Method) của quy hoạch tuyến tính, hạ độ phức tạp xuống mức tối ưu $O(V^3)$ trên đồ thị hai phía đầy đủ.

1. Bản chất toán học: Cặp biến đối ngẫu và Nhãn đỉnh khả thi

Xét đồ thị hai phía đầy đủ $G = (X \cup Y, E)$ với $\vert{}X\vert{} = \vert{}Y\vert{} = N$ và mỗi cạnh $(x, y)$ có trọng số $w(x, y)$.

Bài toán quy hoạch nguyên tìm ma trận ghép $c_{x, y} \in \{0, 1\}$ tối đa hóa tổng trọng số:

$$\max \sum_{x \in X} \sum_{y \in Y} w(x, y) \cdot c_{x, y}$$

thỏa mãn điều kiện mỗi đỉnh chỉ tham gia đúng một cạnh ghép:

$$\sum_{y \in Y} c_{x, y} = 1 \quad (\forall x \in X), \quad \sum_{x \in X} c_{x, y} = 1 \quad (\forall y \in Y)$$

Bài toán đối ngẫu (Dual Problem) và Nhãn đỉnh khả thi

Gán cho mỗi đỉnh $x \in X$ một nhãn $l(x)$ và mỗi đỉnh $y \in Y$ một nhãn $l(y)$. Hàm nhãn $l$ được gọi là khả thi (Feasible Labeling) nếu nó thỏa mãn:

$$l(x) + l(y) \ge w(x, y) \quad (\forall x \in X, y \in Y)$$

Tổng giá trị của bất kỳ ghép cặp hoàn hảo $M$ nào đều bị chặn trên bởi tổng các nhãn đỉnh:

$$\sum_{(x, y) \in M} w(x, y) \le \sum_{(x, y) \in M} (l(x) + l(y)) = \sum_{x \in X} l(x) + \sum_{y \in Y} l(y)$$

Đồ thị đẳng thức (Equality Subgraph)

Từ hàm nhãn khả thi $l$, ta trích xuất Đồ thị đẳng thức $G_l = (X \cup Y, E_l)$ chỉ chứa các cạnh đạt dấu đẳng thức:

$$E_l = \{ (x, y) \in E \mid l(x) + l(y) = w(x, y) \}$$

Định lý Kuhn-Munkres

Định lý: Nếu đồ thị đẳng thức $G_l$ tồn tại một ghép cặp hoàn hảo $M^*$, thì $M^*$ chính là ghép cặp có trọng số lớn nhất trên toàn bộ đồ thị gốc $G$.

Chứng minh: Vì mọi cạnh $(x, y) \in M^*$ đều thuộc $E_l$, ta có:

$$\sum_{(x, y) \in M^*} w(x, y) = \sum_{(x, y) \in M^*} (l(x) + l(y)) = \sum_{x \in X} l(x) + \sum_{y \in Y} l(y)$$

Tổng trọng số của $M^*$ chạm đúng chặn trên của bài toán đối ngẫu, do đó $M^*$ tối ưu tuyệt đối.

2. Cơ chế điều chỉnh nhãn đỉnh và Mảng chùng (Slack Array)

Thuật toán KM bắt đầu bằng một hàm nhãn khả thi khởi tạo:

$$l(x) = \max_{y \in Y} w(x, y) \quad (\forall x \in X), \quad l(y) = 0 \quad (\forall y \in Y)$$

Sau đó, thuật toán lần lượt tìm đường tăng luồng cho từng đỉnh $x \in X$ trên đồ thị đẳng thức $G_l$.

Điều chỉnh nhãn khi không tìm thấy đường tăng luồng

Nếu quá trình tìm kiếm (bằng BFS hoặc DFS) từ một đỉnh $u \in X$ bị tắc nghẽn, gọi $S \subseteq X$ và $T \subseteq Y$ là tập các đỉnh đã được duyệt qua trên cây xen kẽ. Vì không thể tăng luồng, ta cần nới lỏng nhãn để đưa thêm ít nhất một cạnh mới vào đồ thị đẳng thức $G_l$ mà vẫn bảo toàn tính khả thi.

Độ lệch tối thiểu giữa các cặp đỉnh $(x, y)$ với $x \in S$ và $y \notin T$ được xác định bởi:

$$\delta = \min_{x \in S, \, y \notin T} \{ l(x) + l(y) - w(x, y) \}$$

Cập nhật nhãn mới $l'$:

• Với mọi $x \in S$: $l'(x) = l(x) - \delta$

• Với mọi $y \in T$: $l'(y) = l(y) + \delta$

• Các đỉnh còn lại ($x \notin S, y \notin T$) giữ nguyên nhãn.

Hệ quả của phép cập nhật:

1. Cạnh $(x, y)$ với $x \in S, y \in T$: $(l(x) - \delta) + (l(y) + \delta) = l(x) + l(y)$, đẳng thức giữ nguyên.

2. Cạnh $(x, y)$ với $x \notin S, y \notin T$: không đổi, tính khả thi giữ nguyên.

3. Cạnh $(x, y)$ với $x \notin S, y \in T$: tổng nhãn tăng lên $\delta$, tính khả thi vẫn thỏa mãn.

4. Cạnh $(x, y)$ với $x \in S, y \notin T$: tổng nhãn giảm đi $\delta$. Cạnh đạt giá trị cực tiểu $\delta$ sẽ có $l'(x) + l'(y) = w(x, y)$, tức là được thêm mới vào đồ thị đẳng thức $G_l$, mở ra đường đi tiếp theo cho cây xen kẽ.

Kỹ thuật Mảng chùng (Slack Array) để đạt $O(V^3)$

Nếu mỗi lần tìm $\delta$ ta duyệt qua mọi cặp $(x, y) \in S \times (Y \setminus T)$, thao tác này tốn $O(V^2)$, dẫn đến tổng thời gian của thuật toán chạm mức $O(V^4)$.

Để hạ xuống $O(V^3)$, ta duy trì mảng slack[y] cho mỗi $y \notin T$:

$$\text{slack}[y] = \min_{x \in S} \{ l(x) + l(y) - w(x, y) \}$$

• Khi thêm một đỉnh $x$ mới vào $S$: cập nhật slack[y] = min(slack[y], l(x) + l(y) - w(x, y)) cho mọi $y \notin T$ trong $O(V)$.

• Khi tính $\delta$: chỉ cần duyệt qua mảng slack tìm $\min_{y \notin T} \text{slack}[y]$ trong $O(V)$.

• Khi trừ nhãn: cập nhật $\text{slack}[y] \gets \text{slack}[y] - \delta$ cho mọi $y \notin T$ trong $O(V)$.

Nhờ vậy, mỗi lần mở rộng cây xen kẽ hoặc cập nhật nhãn chỉ tiêu tốn $O(V)$, đưa chi phí tìm đường tăng luồng cho một đỉnh về $O(V^2)$, và tổng thời gian cho toàn bộ $N$ đỉnh là $O(V^3)$.

3. Hiện thực hóa Thuật toán Kuhn-Munkres bằng Golang

Chương trình Golang dưới đây cài đặt hoàn chỉnh thuật toán Kuhn-Munkres dạng lặp (Iterative / BFS-based) với mảng slack, hỗ trợ trọng số âm và đồ thị hai phía đầy đủ:

Go

package main

import (
	"fmt"
	"math"
)

const INF = math.MaxInt64 / 4

type Hungarian struct {
	n      int
	weight [][]int64
	lx, ly []int64
	matchY []int   // matchY[y] = x nghĩa là y ghép với x
	matchX []int   // matchX[x] = y
}

func NewHungarian(n int) *Hungarian {
	h := &Hungarian{
		n:      n,
		weight: make([][]int64, n),
		lx:     make([]int64, n),
		ly:     make([]int64, n),
		matchY: make([]int, n),
		matchX: make([]int, n),
	}
	for i := 0; i < n; i++ {
		h.weight[i] = make([]int64, n)
		h.matchY[i] = -1
		h.matchX[i] = -1
	}
	return h
}

func (h *Hungarian) SetWeight(u, v int, w int64) {
	h.weight[u][v] = w
}

// Khởi tạo nhãn đỉnh khả thi ban đầu
func (h *Hungarian) initLabels() {
	for i := 0; i < h.n; i++ {
		h.lx[i] = -INF
		h.ly[i] = 0
		for j := 0; j < h.n; j++ {
			if h.weight[i][j] > h.lx[i] {
				h.lx[i] = h.weight[i][j]
			}
		}
	}
}

// Tìm đường tăng luồng cho một đỉnh u trong O(V^2)
func (h *Hungarian) augment(root int) {
	slack := make([]int64, h.n)
	slackx := make([]int, h.n)
	prev := make([]int, h.n)
	visX := make([]bool, h.n)
	visY := make([]bool, h.n)

	for i := 0; i < h.n; i++ {
		slack[i] = INF
		prev[i] = -1
	}

	queue := []int{root}
	visX[root] = true

	for {
		for len(queue) > 0 {
			x := queue[0]
			queue = queue[1:]

			for y := 0; y < h.n; y++ {
				if !visY[y] {
					delta := h.lx[x] + h.ly[y] - h.weight[x][y]
					if delta < slack[y] {
						slack[y] = delta
						slackx[y] = x
					}
					if slack[y] == 0 {
						visY[y] = true
						prev[y] = slackx[y]
						if h.matchY[y] == -1 {
							// Tìm thấy đường tăng luồng đến đỉnh tự do y
							h.applyAugment(y, prev)
							return
						}
						// Đi tiếp qua cạnh ghép hiện tại của y
						matchedX := h.matchY[y]
						visX[matchedX] = true
						queue = append(queue, matchedX)
					}
				}
			}
		}

		// Không tìm thấy đường đi trong đồ thị đẳng thức hiện tại -> Tính delta
		var delta int64 = INF
		for y := 0; y < h.n; y++ {
			if !visY[y] && slack[y] < delta {
				delta = slack[y]
			}
		}

		// Cập nhật nhãn đỉnh
		for i := 0; i < h.n; i++ {
			if visX[i] {
				h.lx[i] -= delta
			}
			if visY[i] {
				h.ly[i] += delta
			} else {
				slack[i] -= delta
			}
		}

		// Các đỉnh y mới đạt slack[y] == 0 tham gia vào cây xen kẽ
		for y := 0; y < h.n; y++ {
			if !visY[y] && slack[y] == 0 {
				visY[y] = true
				prev[y] = slackx[y]
				if h.matchY[y] == -1 {
					h.applyAugment(y, prev)
					return
				}
				matchedX := h.matchY[y]
				visX[matchedX] = true
				queue = append(queue, matchedX)
			}
		}
	}
}

// Lần ngược đường đi và đảo cạnh ghép
func (h *Hungarian) applyAugment(endY int, prev []int) {
	currY := endY
	for currY != -1 {
		x := prev[currY]
		nextY := h.matchX[x]
		h.matchY[currY] = x
		h.matchX[x] = currY
		currY = nextY
	}
}

// Solve giải bài toán Maximum Weight Perfect Bipartite Matching: O(V^3)
func (h *Hungarian) Solve() (int64, []int) {
	h.initLabels()

	for x := 0; x < h.n; x++ {
		h.augment(x)
	}

	var maxWeight int64 = 0
	for x := 0; x < h.n; x++ {
		y := h.matchX[x]
		maxWeight += h.weight[x][y]
	}

	return maxWeight, h.matchX
}

func main() {
	// Ma trận phân công 4 công nhân (X) cho 4 công việc (Y)
	// Trọng số thể hiện mức độ hiệu quả / lợi nhuận:
	//      Y0  Y1  Y2  Y3
	// X0 [  7,  5,  8,  2 ]
	// X1 [ 10,  9,  4,  5 ]
	// X2 [  3,  8, 12,  6 ]
	// X3 [  6,  3,  5,  7 ]
	n := 4
	matrix := [][]int64{
		{7, 5, 8, 2},
		{10, 9, 4, 5},
		{3, 8, 12, 6},
		{6, 3, 5, 7},
	}

	hungarian := NewHungarian(n)
	for i := 0; i < n; i++ {
		for j := 0; j < n; j++ {
			hungarian.SetWeight(i, j, matrix[i][j])
		}
	}

	totalWeight, matches := hungarian.Solve()

	fmt.Printf("Tổng trọng số cực đại đạt được: %d\n", totalWeight)
	fmt.Println("Phương án phân công tối ưu:")
	for worker, job := range matches {
		fmt.Printf("- Công nhân X%d ===> Công việc Y%d (Hiệu quả: %d)\n",
			worker, job, matrix[worker][job])
	}
}

4. Biến thể tìm Trọng số cực tiểu (Minimum Weight Matching)

Nếu bài toán yêu cầu tìm phân công có chi phí nhỏ nhất thay vì lợi nhuận lớn nhất, ta có hai cách chuyển đổi đơn giản:

1. Đổi dấu ma trận trọng số: Đặt $w'(x, y) = -w(x, y)$. Chạy thuật toán tìm Maximum Weight trên $w'$, sau đó đảo dấu kết quả tổng thu được:

   $$\min \sum w(x, y) = - \max \sum (-w(x, y))$$

2. Trừ vào hằng số lớn: Chọn $M = \max_{x, y} w(x, y)$, định nghĩa $w'(x, y) = M - w(x, y) \ge 0$. Khi đó:

   $$\sum w'(x, y) = N \cdot M - \sum w(x, y)$$

   Cực đại hóa $\sum w'$ đồng nghĩa với cực tiểu hóa $\sum w$.

5. So sánh hiệu năng giữa các thuật toán giải Assignment Problem

6. Ứng dụng thực tế

• Hệ thống theo dõi đa đối tượng (Multi-Object Tracking - MOT / DeepSORT): Trong thị giác máy tính và xe tự hành, thuật toán Hungarian được thực thi ở mỗi khung hình để ghép các bounding box vừa phát hiện (YOLO detections) với quỹ đạo dự đoán của bộ lọc Kalman (Kalman filter tracks) dựa trên chỉ số IoU hoặc khoảng cách Cosine.

• Hệ thống phân bổ cuốc xe công nghệ (Ride-Hailing Dispatch): Ghép cặp đồng thời $N$ tài xế với $N$ hành khách đặt xe trong cùng một khu vực địa lý để tối thiểu hóa thời gian chờ đón xe trung bình toàn hệ thống.

• Khớp lệnh thị trường đấu thầu (Combinatorial Auction & Ad-Exchange): Ghép các lượt hiển thị quảng cáo (ad impressions) với các nhà quảng cáo đặt giá thầu cao nhất thỏa mãn các ràng buộc ngân sách theo thời gian thực.