Cấu trúc dữ liệu và Giải thuật Bài 145: Thuật toán xấp xỉ Christofides-Serdyukov (3/2-Approximation for Metric TSP) - Cây khung và Ghép cặp trọng số cực tiểu.

Ở Bài 144, chúng ta đã tiếp cận bài toán Người bán hàng (TSP) bằng các thuật toán chính xác như Quy hoạch động Held-Karp ($O(N^2 2^N)$) và Phân nhánh chặn (Branch and Bound). Tuy nhiên, vì TSP thuộc lớp bài toán NP-hard, việc tìm nghiệm chính xác tuyệt đối trên các tập dữ liệu thực tế cỡ lớn ($N > 1.000$) đòi hỏi tài nguyên tính toán vượt quá khả năng xử lý trong thời gian thực.

Đối với TSP tổng quát, việc tìm kiếm một thuật toán xấp xỉ với bất kỳ hệ số hằng số nào cũng là điều bất khả thi (NP-hard to approximate within $2^{O(N)}$). Tuy nhiên, trong thực tế vật lý và hình học, ma trận chi phí thường tuân thủ Bất đẳng thức tam giác (Triangle Inequality): "Đi thẳng từ $u$ đến $v$ không bao giờ tốn kém hơn đi đường vòng qua $x$" ($w(u, v) \le w(u, x) + w(x, v)$). Biến thể này được gọi là Metric TSP.

Năm 1976, Nicos Christofides (và độc lập bởi Anatoliy Serdyukov vào năm 1978) đã công bố Thuật toán Christofides-Serdyukov. Bằng cách kết hợp Cây khung nhỏ nhất (MST) và Ghép cặp hoàn hảo trọng số cực tiểu (Minimum-Weight Perfect Matching - MWPM), thuật toán đảm bảo tìm ra hành trình có chi phí không vượt quá $1.5 \times \text{OPT}$ trong thời gian đa thức $O(N^3)$. Đây là một trong những cột mốc vĩ đại nhất của lý thuyết thuật toán xấp xỉ trong thế kỷ 20.

1. Cơ sở lý thuyết: Từ 2-Approximation đến 3/2-Approximation

Phương pháp 2-xấp xỉ cơ bản (Double-Tree Algorithm)

Cách tiếp cận ngây thơ nhất trên Metric TSP là:

1. Dựng Cây khung nhỏ nhất (MST) $T$ của đồ thị. Trọng số $w(T) \le \text{OPT}$ (vì bỏ 1 cạnh khỏi chu trình TSP tối ưu ta được 1 cây khung).

2. Nhân đôi mọi cạnh của $T$ thành đa đồ thị $2T$. Khi đó mọi đỉnh đều có bậc chẵn $\implies 2T$ là đồ thị Euler.

3. Đi theo chu trình Euler và "đi tắt" (shortcut) qua các đỉnh đã thăm để tạo chu trình Hamilton $C$.

4. Do bất đẳng thức tam giác, việc shortcutting không làm tăng chi phí:

   $$w(C) \le w(2T) = 2 \cdot w(T) \le 2 \cdot \text{OPT}$$

        Thuật toán Double-Tree (2-OPT):
        Nhân đôi TẤT CẢ các cạnh của MST:
        (0) ===== (1) ===== (2)    ===>    Tổng trọng số: 2 * w(T) <= 2 * OPT

Ý tưởng đột phá của Christofides-Serdyukov

Tại sao phải nhân đôi toàn bộ cây khung khi mục tiêu duy nhất chỉ là làm cho bậc của mọi đỉnh trở thành số chẵn?

• Trong một cây $T$, chỉ có một tập hợp các đỉnh con $O \subset V$ có bậc lẻ (Odd-degree vertices).

• Theo Bổ đề Bắt tay (Handshaking Lemma), số lượng đỉnh có bậc lẻ $\vert{}O\vert{}$ luôn luôn là một số chẵn.

• Thay vì nhân đôi toàn bộ $T$, ta chỉ cần tìm một Ghép cặp hoàn hảo có trọng số nhỏ nhất (Minimum-Weight Perfect Matching - MWPM) trên tập đỉnh $O$ rồi bổ sung các cạnh ghép đó vào $T$.

• Sau khi ghép, mỗi đỉnh trong $O$ được tăng bậc thêm đúng 1, trở thành bậc chẵn. Các đỉnh vốn đã có bậc chẵn trong $T$ không bị ảnh hưởng.

• Đa đồ thị mới trở thành đồ thị Euler với số lượng cạnh bổ sung ít hơn rất nhiều.

        Thuật toán Christofides (1.5-OPT):
        Chỉ thêm các cạnh ghép cặp hoàn hảo giữa các đỉnh BẬC LẺ:
        (u) [lẻ] ------------- (v) [lẻ]    ===>    Tổng trọng số: w(T) + w(M) <= 1.5 * OPT

2. Chứng minh Toán học: Hệ số Xấp xỉ $3/2$ ($\alpha = 1.5$)

Gọi $C^*$ là chu trình TSP tối ưu toàn cục với chi phí $\text{OPT} = w(C^*)$.

Bổ đề 1: Trọng số của Cây khung nhỏ nhất bị chặn bởi $\text{OPT}$

Xóa bỏ một cạnh bất kỳ khỏi chu trình tối ưu $C^*$, ta thu được một đường đi Hamilton phủ toàn bộ $N$ đỉnh. Đường đi này là một cây khung của $G$. Do $T$ là cây khung có trọng số nhỏ nhất:

$$w(T) \le w(C^* \setminus \{e\}) < w(C^*) = \text{OPT}$$

Bổ đề 2: Trọng số của Ghép cặp hoàn hảo MWPM trên $O$ bị chặn bởi $\frac{1}{2} \text{OPT}$

1. Xét chu trình tối ưu $C^*$. Rút gọn $C^*$ bằng cách đi tắt (shortcutting) chỉ qua các đỉnh thuộc tập $O$ theo đúng thứ tự xuất hiện của chúng trên $C^*$.

2. Ta thu được một chu trình $C_O$ chỉ đi qua các đỉnh của $O$. Do bất đẳng thức tam giác:

   $$w(C_O) \le w(C^*) = \text{OPT}$$

3. Vì $\vert{}O\vert{}$ là số chẵn, chu trình $C_O$ có thể phân rã thành hai tập ghép cặp hoàn hảo rời nhau về cạnh: $M_1$ (các cạnh ở vị trí lẻ) và $M_2$ (các cạnh ở vị trí chẵn).

   $$w(M_1) + w(M_2) = w(C_O) \le \text{OPT}$$

4. Do $M$ là ghép cặp hoàn hảo có trọng số nhỏ nhất trên $O$:

   $$w(M) \le \min(w(M_1), w(M_2)) \le \frac{w(M_1) + w(M_2)}{2} \le \frac{1}{2} \text{OPT}$$

Kết luận Định lý

Hợp của cây $T$ và tập ghép $M$ tạo thành một đa đồ thị Euler $H = T \cup M$ có tổng trọng số:

$$w(H) = w(T) + w(M) \le \text{OPT} + \frac{1}{2} \text{OPT} = \frac{3}{2} \text{OPT}$$

Khi chuyển chu trình Euler trên $H$ thành chu trình Hamilton $C$ bằng phép đi tắt, bất đẳng thức tam giác bảo toàn:

$$w(C) \le w(H) \le \frac{3}{2} \text{OPT} = 1.5 \times \text{OPT}$$

3. Quy trình 6 bước của Thuật toán Christofides-Serdyukov

1. Bước 1 (Minimum Spanning Tree): Xây dựng Cây khung nhỏ nhất $T$ của đồ thị đầy đủ $G$ bằng thuật toán Kruskal hoặc Prim.

2. Bước 2 (Identify Odd Vertices): Đếm bậc của các đỉnh trên cây $T$ và trích xuất tập các đỉnh có bậc lẻ $O = \{v \in V \mid \deg_T(v) \equiv 1 \pmod 2\}$.

3. Bước 3 (Minimum-Weight Perfect Matching): Dựng đồ thị con cảm ứng trên tập đỉnh $O$ và tìm Ghép cặp hoàn hảo có trọng số cực tiểu $M$.

4. Bước 4 (Eulerian Multigraph): Hợp nhất các cạnh của $T$ và các cạnh của $M$ thành đa đồ thị $H$. Mọi đỉnh trong $H$ đều có bậc chẵn.

5. Bước 5 (Eulerian Tour): Tìm chu trình Euler trên $H$ bằng thuật toán Hierholzer.

6. Bước 6 (Shortcutting): Duyệt tuần tự chu trình Euler, bỏ qua các đỉnh đã xuất hiện trước đó, tạo thành chu trình Hamilton khép kín.

4. Hiện thực hóa Thuật toán bằng Golang

Đoạn mã Golang dưới đây cài đặt hoàn chỉnh thuật toán Christofides:

• Hỗ trợ tọa độ 2D Euclidean (tự động đảm bảo tính Metric).

• Thuật toán Kruskal với DSU để dựng MST.

• Thuật toán Quy hoạch động Bitmask để giải quyết chính xác bài toán Ghép cặp hoàn hảo cực tiểu (MWPM) trên tập các đỉnh bậc lẻ.

• Thuật toán Hierholzer tìm chu trình Euler và kỹ thuật Shortcutting.

Go

package main

import (
	"fmt"
	"math"
	"sort"
)

const INF = math.MaxInt64 / 4

type Point struct {
	x, y float64
}

func distance(a, b Point) int64 {
	dx := a.x - b.x
	dy := a.y - b.y
	return int64(math.Round(math.Sqrt(dx*dx + dy*dy)))
}

type Edge struct {
	u, v int
	w    int64
}

// Cấu trúc DSU cho Kruskal
type DSU struct {
	parent []int
}

func NewDSU(n int) *DSU {
	p := make([]int, n)
	for i := range p {
		p[i] = i
	}
	return &DSU{parent: p}
}

func (d *DSU) Find(i int) int {
	if d.parent[i] == i {
		return i
	}
	d.parent[i] = d.Find(d.parent[i])
	return d.parent[i]
}

func (d *DSU) Union(i, j int) bool {
	rootI := d.Find(i)
	rootJ := d.Find(j)
	if rootI != rootJ {
		d.parent[rootI] = rootJ
		return true
	}
	return false
}

type ChristofidesTSP struct {
	n      int
	points []Point
	dist   [][]int64
}

func NewChristofidesTSP(points []Point) *ChristofidesTSP {
	n := len(points)
	dist := make([][]int64, n)
	for i := 0; i < n; i++ {
		dist[i] = make([]int64, n)
		for j := 0; j < n; j++ {
			dist[i][j] = distance(points[i], points[j])
		}
	}
	return &ChristofidesTSP{
		n:      n,
		points: points,
		dist:   dist,
	}
}

// Bước 1: Dựng MST bằng Kruskal
func (c *ChristofidesTSP) computeMST() []Edge {
	var edges []Edge
	for i := 0; i < c.n; i++ {
		for j := i + 1; j < c.n; j++ {
			edges = append(edges, Edge{u: i, v: j, w: c.dist[i][j]})
		}
	}

	sort.Slice(edges, func(i, j int) bool {
		return edges[i].w < edges[j].w
	})

	dsu := NewDSU(c.n)
	var mstEdges []Edge

	for _, e := range edges {
		if dsu.Union(e.u, e.v) {
			mstEdges = append(mstEdges, e)
			if len(mstEdges) == c.n-1 {
				break
			}
		}
	}
	return mstEdges
}

// Bước 2: Tìm các đỉnh có bậc lẻ
func (c *ChristofidesTSP) findOddVertices(mst []Edge) []int {
	deg := make([]int, c.n)
	for _, e := range mst {
		deg[e.u]++
		deg[e.v]++
	}

	var odd []int
	for i := 0; i < c.n; i++ {
		if deg[i]%2 != 0 {
			odd = append(odd, i)
		}
	}
	return odd
}

// Bước 3: Tìm Ghép cặp hoàn hảo trọng số cực tiểu (MWPM) trên tập đỉnh lẻ bằng DP Bitmask
func (c *ChristofidesTSP) computeMWPM(odd []int) []Edge {
	k := len(odd)
	if k == 0 {
		return nil
	}

	numStates := 1 << k
	dp := make([]int64, numStates)
	partner := make([]int, numStates) // Lưu vết ghép cặp
	for i := range dp {
		dp[i] = INF
	}
	dp[0] = 0

	for mask := 0; mask < numStates; mask++ {
		if dp[mask] == INF {
			continue
		}

		// Tìm bit 0 đầu tiên chưa được ghép
		first := -1
		for i := 0; i < k; i++ {
			if (mask & (1 << i)) == 0 {
				first = i
				break
			}
		}
		if first == -1 {
			continue
		}

		// Thử ghép 'first' với một đỉnh 'j' khác chưa được ghép
		for j := first + 1; j < k; j++ {
			if (mask & (1 << j)) == 0 {
				nextMask := mask | (1 << first) | (1 << j)
				cost := dp[mask] + c.dist[odd[first]][odd[j]]
				if cost < dp[nextMask] {
					dp[nextMask] = cost
					partner[nextMask] = (first << 8) | j
				}
			}
		}
	}

	// Lần vết khôi phục các cặp ghép
	var matchedEdges []Edge
	currMask := numStates - 1
	for currMask > 0 {
		val := partner[currMask]
		first := val >> 8
		j := val & 0xFF
		matchedEdges = append(matchedEdges, Edge{
			u: odd[first],
			v: odd[j],
			w: c.dist[odd[first]][odd[j]],
		})
		currMask ^= (1 << first) | (1 << j)
	}

	return matchedEdges
}

// Bước 4 & 5: Dựng chu trình Euler trên đa đồ thị bằng thuật toán Hierholzer
func (c *ChristofidesTSP) eulerianTour(edges []Edge) []int {
	type AdjEdge struct {
		to int
		id int
	}
	adj := make([][]AdjEdge, c.n)
	for id, e := range edges {
		adj[e.u] = append(adj[e.u], AdjEdge{to: e.v, id: id})
		adj[e.v] = append(adj[e.v], AdjEdge{to: e.u, id: id})
	}

	used := make([]bool, len(edges))
	var tour []int
	var stack []int
	stack = append(stack, edges[0].u)

	for len(stack) > 0 {
		u := stack[len(stack)-1]
		foundEdge := false

		for len(adj[u]) > 0 {
			top := adj[u][len(adj[u])-1]
			adj[u] = adj[u][:len(adj[u])-1]

			if !used[top.id] {
				used[top.id] = true
				stack = append(stack, top.to)
				foundEdge = true
				break
			}
		}

		if !foundEdge {
			tour = append(tour, u)
			stack = stack[:len(stack)-1]
		}
	}

	return tour
}

// Bước 6: Shortcutting chuyển chu trình Euler thành chu trình Hamilton
func (c *ChristofidesTSP) shortcut(eulerTour []int) ([]int, int64) {
	visited := make([]bool, c.n)
	var hamiltonTour []int

	for _, u := range eulerTour {
		if !visited[u] {
			visited[u] = true
			hamiltonTour = append(hamiltonTour, u)
		}
	}
	// Khép kín chu trình
	hamiltonTour = append(hamiltonTour, hamiltonTour[0])

	var totalCost int64 = 0
	for i := 0; i < len(hamiltonTour)-1; i++ {
		totalCost += c.dist[hamiltonTour[i]][hamiltonTour[i+1]]
	}

	return hamiltonTour, totalCost
}

// Solve thực thi toàn bộ pipeline của Christofides
func (c *ChristofidesTSP) Solve() ([]int, int64) {
	mst := c.computeMST()
	odd := c.findOddVertices(mst)
	matching := c.computeMWPM(odd)

	multigraph := append(mst, matching...)
	euler := c.eulerianTour(multigraph)
	hamilton, cost := c.shortcut(euler)

	return hamilton, cost
}

func main() {
	// Khởi tạo tọa độ 2D của 8 thành phố
	points := []Point{
		{x: 0, y: 0},
		{x: 1, y: 5},
		{x: 5, y: 6},
		{x: 8, y: 4},
		{x: 7, y: 1},
		{x: 4, y: 0},
		{x: 3, y: 3},
		{x: 6, y: 3},
	}

	solver := NewChristofidesTSP(points)
	tour, cost := solver.Solve()

	fmt.Printf("=== KẾT QUẢ THUẬT TOÁN CHRISTOFIDES-SERDYUKOV (N = %d) ===\n", len(points))
	fmt.Printf("Chi phí hành trình xấp xỉ (<= 1.5 * OPT): %d\n\n", cost)

	fmt.Print("Chu trình Hamilton thu được: ")
	for i, city := range tour {
		if i > 0 {
			fmt.Print(" -> ")
		}
		fmt.Printf("P%d", city)
	}
	fmt.Println()

	fmt.Println("\nChi tiết tọa độ các điểm trên hành trình:")
	for _, city := range tour {
		fmt.Printf("- P%d: (%.1f, %.1f)\n", city, points[city].x, points[city].y)
	}
}

5. Phân tích Độ phức tạp Thuật toán

• Nút thắt thời gian: Khâu tốn kém nhất là tính toán Ghép cặp hoàn hảo trọng số cực tiểu (MWPM) trên tập đỉnh lẻ. Khi sử dụng thuật toán Edmonds' Blossom (Bài 130), toàn bộ thuật toán chạy trong thời gian đa thức chặt chẽ $O(N^3)$.

6. Kỷ lục Lịch sử: 44 Năm Bất bại và Đột phá Năm 2020

• Kỷ lục 44 năm: Kể từ năm 1976, hệ số $3/2 = 1.5$ của Christofides đứng vững như một "bức tường thành" không thể xô đổ trong lý thuyết tối ưu hóa tổ hợp. Hàng nghìn công trình nghiên cứu đều thất bại trong việc cải thiện hệ số này dù chỉ là $\epsilon$ nhỏ nhất cho Metric TSP tổng quát.

• Cột mốc lịch sử 2020: Anna Karlin, Nathan Klein và Shayan Oveis Gharan (Đại học Washington) đã công bố công trình chấn động, chứng minh tồn tại thuật toán xấp xỉ đạt hệ số:

  $$\alpha = 1.5 - 10^{-36}$$

  dựa trên phương pháp lấy mẫu cây khung ngẫu nhiên phân phối đều (Randomized Tree Sampling via Strongly Rayleigh Measures). Dù giá trị cải thiện rất nhỏ, đây là bước nhảy lý thuyết phá vỡ giới hạn tưởng chừng bất biến của khoa học máy tính.

7. Ứng dụng Thực tế

• Quy hoạch Định tuyến Xe tự hành & Giao hàng (Vehicle Routing Problem - VRP): Được sử dụng như thuật toán sinh lời giải khởi tạo (Initial Solution Generation) chất lượng cao trước khi đưa vào các bộ tối ưu cục bộ nâng cao như 2-Opt, 3-Opt hoặc Large Neighborhood Search (LNS).

• Vẽ đường viền công nghiệp bằng Laser & Cắt tia nước (Laser / Waterjet Path Planning): Giảm thiểu thời gian đầu cắt di chuyển không tải giữa các đường viền khép kín độc lập trên tấm kim loại tấm lớn.

• Tối ưu hóa Robot thám hiểm & Quét bề mặt 3D: Lập trình lộ trình di chuyển camera qua tập hợp các góc quan sát cố định sao cho hành trình khép kín ngắn nhất, tối ưu hóa thời lượng pin của thiết bị không người lái (Drone).