Cấu trúc dữ liệu và Giải thuật Bài 146: Cắt cực đại xấp xỉ (Max-Cut - Thuật toán Goemans-Williamson) - Phá vỡ giới hạn NP-hard bằng Quy hoạch bán xác định (SDP Relaxation) và Siêu phẳng ngẫu nhiên.

Trong lý thuyết đồ thị và tối ưu hóa tổ hợp, bài toán Lát cắt cực tiểu (Min-Cut) có thể giải quyết chính xác trong thời gian đa thức bằng các thuật toán mạng luồng hoặc co cụm tất định như Stoer-Wagner (Bài 137). Tuy nhiên, bài toán đối ngẫu của nó — Lát cắt cực đại (Maximum Cut / Max-Cut) — lại thuộc danh sách 21 bài toán NP-complete kinh điển của Richard Karp (1972).

Một thuật toán ngẫu nhiên hóa ngây thơ (gán ngẫu nhiên mỗi đỉnh vào tập $S$ hoặc $V \setminus S$ với xác suất $1/2$) hoặc thuật toán tham lam tìm kiếm cục bộ chỉ có thể đạt được tỷ lệ xấp xỉ tầm thường là $\alpha = 0.5$. Suốt hơn hai thập kỷ, mốc $0.5$ được coi là giới hạn không thể vượt qua đối với các phương pháp giải đa thức.

Mãi đến năm 1995, Michel Goemans và David P. Williamson đã tạo nên một bước ngoặt lịch sử: sử dụng Quy hoạch bán xác định (Semidefinite Programming - SDP Relaxation) để nới lỏng không gian biến nguyên thành các vector đơn vị nhiều chiều, sau đó kết hợp kỹ thuật Làm tròn bằng Siêu phẳng ngẫu nhiên (Random Hyperplane Rounding) để đạt tỷ lệ xấp xỉ kỷ lục:

$$\alpha_{\text{GW}} = \min_{0 < \theta \le \pi} \frac{2}{\pi} \frac{\theta}{1 - \cos \theta} \approx 0.878567$$

Đặc biệt, theo Giả thuyết Trò chơi Độc nhất (Unique Games Conjecture - Khot et al., 2007), $0.878567$ chính là giới hạn xấp xỉ tối thượng: không thể tồn tại bất kỳ thuật toán thời gian đa thức nào đạt tỷ lệ tốt hơn trừ khi $P = NP$.

1. Định nghĩa Toán học của Bài toán Max-Cut

Cho đồ thị vô hướng $G = (V, E)$ gồm $n = \vert{}V\vert{}$ đỉnh với trọng số các cạnh không âm $w_{i, j} \ge 0$ (nếu không có cạnh giữa $i$ và $j$, $w_{i, j} = 0$).

Một lát cắt là sự phân hoạch tập đỉnh thành hai tập rời nhau $S$ và $\bar{S} = V \setminus S$. Trọng số của lát cắt là tổng trọng số các cạnh có một đầu mút thuộc $S$ và đầu mút kia thuộc $\bar{S}$:

$$w(S, \bar{S}) = \sum_{i \in S, \, j \in \bar{S}} w_{i, j}$$

Mục tiêu là tìm tập $S \subset V$ sao cho $w(S, \bar{S})$ đạt giá trị cực đại toàn cục.

Mô hình hóa dưới dạng Quy hoạch Số nguyên Toàn phương (Integer Quadratic Program - IQP)

Gán cho mỗi đỉnh $i \in V$ một biến nhị phân $x_i \in \{-1, +1\}$:

• $x_i = +1$ nếu $i \in S$.

• $x_i = -1$ nếu $i \in \bar{S}$.

Khi đó:

• Nếu $i$ và $j$ nằm cùng một phía: $x_i x_j = (+1)(+1) = (-1)(-1) = 1 \implies \frac{1 - x_i x_j}{2} = 0$.

• Nếu $i$ và $j$ nằm ở hai phía khác nhau (cạnh bị cắt): $x_i x_j = -1 \implies \frac{1 - x_i x_j}{2} = 1$.

Bài toán Max-Cut được viết lại chính xác dưới dạng:

$$\max \frac{1}{2} \sum_{i < j} w_{i, j} (1 - x_i x_j) \quad \text{sao cho} \quad x_i \in \{-1, +1\} \quad (\forall i \in V)$$

2. Nới lỏng Quy hoạch Bán xác định (SDP Relaxation)

Rào cản khiến bài toán IQP trên thuộc lớp NP-hard là không gian tìm kiếm bị cô lập trên trục số thực 1 chiều: $x_i \in \{-1, +1\}$.

Bước nhảy không gian: Từ vô hướng 1D sang Vector $n$ chiều

Goemans và Williamson nới lỏng bài toán bằng cách thay thế mỗi biến vô hướng $x_i \in \{-1, +1\}$ bằng một vector đơn vị $v_i \in \mathbb{R}^n$ nằm trên mặt cầu đơn vị $S^{n-1}$ ($\Vert{}v_i\Vert{} = 1$). Tích hai biến $x_i x_j$ được thay thế bằng tích vô hướng $v_i \cdot v_j$:

$$\max \frac{1}{2} \sum_{i < j} w_{i, j} (1 - v_i \cdot v_j) \quad \text{sao cho} \quad \Vert{}v_i\Vert{}^2 = v_i \cdot v_i = 1 \quad (\forall i \in V)$$

   1D (Biến nguyên NP-hard):        n-Chiều (Nới lỏng SDP - Đa thức):
        -1              +1                         v2
   <-----o---------------o----->                 /   \
        x1              x2                     v1 --- v3  (Nằm trên mặt cầu S^{n-1})

Dạng Ma trận Bán xác định (Semidefinite Matrix Formulation)

Đặt ma trận Gram $Y \in \mathbb{R}^{n \times n}$ với các phần tử $Y_{i, j} = v_i \cdot v_j$. Ma trận $Y$ là đối xứng và bán xác định dương ($Y \succeq 0$). Bài toán trở thành:

$$\max \frac{1}{2} \sum_{i < j} w_{i, j} (1 - Y_{i, j})$$

$$\text{sao cho} \quad Y_{i, i} = 1 \quad (\forall i = 1, \dots, n), \quad Y \succeq 0$$

Vì hàm mục tiêu là tuyến tính theo $Y$ và tập các ma trận bán xác định dương tạo thành một nón lồi (Convex Cone), bài toán này có thể giải đến độ chính xác tùy ý $\epsilon$ trong thời gian đa thức bằng Phương pháp Điểm trong (Interior Point Methods) hoặc Phân rã Burer-Monteiro.

Gọi giá trị tối ưu của bài toán SDP này là $\text{SDP}_{\text{opt}}$. Do không gian vector rộng hơn không gian nhị phân $\{-1, +1\}$, ta luôn có:

$$\text{SDP}_{\text{opt}} \ge \text{OPT}$$

3. Kỹ thuật Làm tròn bằng Siêu phẳng Ngẫu nhiên (Random Hyperplane Rounding)

Sau khi giải bài toán SDP, ta thu được tập hợp các vector đơn vị tối ưu $\{v_1, v_2, \dots, v_n\}$ phân bố trên mặt cầu đơn vị $\mathbb{R}^n$ (thông qua phép phân tích Cholesky $Y = V^T V$). Làm thế nào để ánh xạ các vector này trở lại các nhãn nhị phân $\{-1, +1\}$?

Cơ chế Siêu phẳng đi qua gốc tọa độ

1. Sinh ngẫu nhiên một vector $r \in \mathbb{R}^n$ từ phân phối chuẩn nhiều chiều $\mathcal{N}(0, I_n)$ và chuẩn hóa về độ dài $1$. Vector $r$ đóng vai trò là vector pháp tuyến của một siêu phẳng ngẫu nhiên đi qua gốc tọa độ.

2. Siêu phẳng này chia không gian $\mathbb{R}^n$ thành hai nửa bán không gian:

   • Nếu $v_i \cdot r \ge 0$: Gán $x_i = +1$ (cho $i \in S$).

   • Nếu $v_i \cdot r < 0$: Gán $x_i = -1$ (cho $i \in \bar{S}$).

                             Siêu phẳng ngẫu nhiên (pháp tuyến r)
                                      /
                                v1   /   v2
                                 o  /   o
                                   /
                      ------------+------------
                                 /  o v3
                                /
                               /
                 Nửa không gian (-) | Nửa không gian (+)

Chứng minh Hệ số Xấp xỉ $\alpha_{\text{GW}} \approx 0.878567$

Xét một cặp đỉnh bất kỳ $(i, j)$ với góc tạo bởi hai vector $v_i$ và $v_j$ là $\theta_{i, j} = \arccos(v_i \cdot v_j) \in [0, \pi]$.

1. Siêu phẳng ngẫu nhiên sẽ phân tách hai vector $v_i$ và $v_j$ vào hai phía khác nhau khi và chỉ khi nó cắt qua góc phẳng tạo bởi $v_i$ và $v_j$.

2. Do tính chất đối xứng cầu của phân phối chuẩn, xác suất để siêu phẳng cắt cặp cạnh này tỷ lệ thuận với góc giữa chúng:

   $$\mathbb{P}[x_i \neq x_j] = \frac{\theta_{i, j}}{\pi} = \frac{\arccos(v_i \cdot v_j)}{\pi}$$

3. Kỳ vọng tổng trọng số của lát cắt thu được là:

   $$\mathbb{E}[w(S, \bar{S})] = \sum_{i < j} w_{i, j} \cdot \frac{\arccos(v_i \cdot v_j)}{\pi}$$

4. Đóng góp của cạnh $(i, j)$ trong hàm mục tiêu SDP là $\frac{1}{2}(1 - v_i \cdot v_j) = \frac{1 - \cos \theta_{i, j}}{2}$.

5. Tỷ số giữa đóng góp thực tế và đóng góp SDP trên từng cạnh bị chặn dưới bởi hàm số:

   $$\alpha(\theta) = \frac{\theta / \pi}{(1 - \cos \theta) / 2} = \frac{2}{\pi} \frac{\theta}{1 - \cos \theta}$$

Khảo sát hàm số trên khoảng $\theta \in (0, \pi]$:

• Đạo hàm triệt tiêu tại nghiệm của phương trình $\tan \theta = \theta$, cho giá trị cực tiểu tại $\theta^* \approx 2.331122 \text{ rad } (\approx 133.56^\circ)$.

• Giá trị cực tiểu đạt được chính xác là:

  $$\alpha_{\text{GW}} = \frac{2}{\pi} \frac{\theta^*}{1 - \cos \theta^*} \approx 0.878567$$

Do bất đẳng thức đúng trên từng số hạng, ta suy ra:

$$\mathbb{E}[w(S, \bar{S})] \ge \alpha_{\text{GW}} \cdot \text{SDP}_{\text{opt}} \ge \alpha_{\text{GW}} \cdot \text{OPT} \approx 0.878567 \cdot \text{OPT}$$

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

Để chương trình chạy độc lập mà không cần thư viện bộ giải SDP bên ngoài cồng kềnh, chúng ta sử dụng phương pháp Burer-Monteiro (Riemannian Projected Gradient Descent): biểu diễn trực tiếp các vector $v_i \in \mathbb{R}^d$ ($d \approx \sqrt{2n}$ hoặc $d = n$), tối ưu hóa trên đa tạp hình cầu (Spherical Manifold) để tìm cấu hình vector tối ưu, sau đó áp dụng phép chiếu siêu phẳng ngẫu nhiên:

Go

package main

import (
	"fmt"
	"math"
	"math/rand"
	"time"
)

type Edge struct {
	u, v int
	w    float64
}

type MaxCutGW struct {
	n        int
	edges    []Edge
	adj      [][]float64
	dim      int         // Chiều không gian vector (d <= n)
	vectors  [][]float64 // v[i] là vector đơn vị d-chiều của đỉnh i
	rng      *rand.Rand
}

func NewMaxCutGW(n int, edges []Edge) *MaxCutGW {
	adj := make([][]float64, n)
	for i := 0; i < n; i++ {
		adj[i] = make([]float64, n)
	}
	for _, e := range edges {
		adj[e.u][e.v] += e.w
		adj[e.v][e.u] += e.w
	}

	dim := n
	if dim < 10 {
		dim = 10
	}

	rng := rand.New(rand.NewSource(time.Now().UnixNano()))

	// Khởi tạo ngẫu nhiên các vector đơn vị v_i trên mặt cầu S^{dim-1}
	vectors := make([][]float64, n)
	for i := 0; i < n; i++ {
		vectors[i] = make([]float64, dim)
		var norm float64
		for k := 0; k < dim; k++ {
			vectors[i][k] = rng.NormFloat64()
			norm += vectors[i][k] * vectors[i][k]
		}
		norm = math.Sqrt(norm)
		for k := 0; k < dim; k++ {
			vectors[i][k] /= norm
		}
	}

	return &MaxCutGW{
		n:       n,
		edges:   edges,
		adj:     adj,
		dim:     dim,
		vectors: vectors,
		rng:     rng,
	}
}

// Tích vô hướng giữa hai vector
func dotProduct(a, b []float64) float64 {
	var sum float64
	for i := range a {
		sum += a[i] * b[i]
	}
	return sum
}

// Tối ưu hóa cấu hình vector bằng Riemannian Projected Gradient Descent (Burer-Monteiro)
// Cực đại hóa: sum w_{ij} * (1 - v_i . v_j) <=> Cực tiểu hóa: sum w_{ij} * (v_i . v_j)
func (gw *MaxCutGW) SolveVectorEmbedding(iterations int, lr float64) {
	for iter := 0; iter < iterations; iter++ {
		for i := 0; i < gw.n; i++ {
			// Gradient của hàm chi phí theo v_i: grad = sum_{j} w_{ij} * v_j
			grad := make([]float64, gw.dim)
			for j := 0; j < gw.n; j++ {
				if gw.adj[i][j] > 0 {
					w := gw.adj[i][j]
					for k := 0; k < gw.dim; k++ {
						grad[k] += w * gw.vectors[j][k]
					}
				}
			}

			// Chiếu gradient lên không gian tiếp tuyến của mặt cầu (Riemannian Gradient):
			// proj_grad = grad - (grad . v_i) * v_i
			projCoeff := dotProduct(grad, gw.vectors[i])
			for k := 0; k < gw.dim; k++ {
				grad[k] -= projCoeff * gw.vectors[i][k]
			}

			// Cập nhật vị trí vector (Gradient Descent để giảm sum w_ij * v_i.v_j)
			var newNorm float64
			for k := 0; k < gw.dim; k++ {
				gw.vectors[i][k] -= lr * grad[k]
				newNorm += gw.vectors[i][k] * gw.vectors[i][k]
			}

			// Chiếu ngược lại mặt cầu đơn vị: ||v_i|| = 1
			newNorm = math.Sqrt(newNorm)
			for k := 0; k < gw.dim; k++ {
				gw.vectors[i][k] /= newNorm
			}
		}
	}
}

// Giá trị chặn trên SDP
func (gw *MaxCutGW) ComputeSDPValue() float64 {
	var val float64
	for _, e := range gw.edges {
		dot := dotProduct(gw.vectors[e.u], gw.vectors[e.v])
		val += 0.5 * e.w * (1.0 - dot)
	}
	return val
}

// Thực hiện một lượt làm tròn bằng Siêu phẳng ngẫu nhiên (Random Hyperplane Rounding)
func (gw *MaxCutGW) RandomHyperplaneRound() ([]bool, float64) {
	// 1. Sinh vector pháp tuyến ngẫu nhiên r từ phân phối chuẩn N(0, I)
	r := make([]float64, gw.dim)
	for k := 0; k < gw.dim; k++ {
		r[k] = gw.rng.NormFloat64()
	}

	// 2. Phân hoạch đỉnh dựa trên dấu của tích vô hướng v_i . r
	cutAssignment := make([]bool, gw.n) // true nếu thuộc tập S, false nếu thuộc V \ S
	for i := 0; i < gw.n; i++ {
		cutAssignment[i] = dotProduct(gw.vectors[i], r) >= 0
	}

	// 3. Tính tổng trọng số các cạnh bị cắt
	var cutWeight float64
	for _, e := range gw.edges {
		if cutAssignment[e.u] != cutAssignment[e.v] {
			cutWeight += e.w
		}
	}

	return cutAssignment, cutWeight
}

// Thực hiện lặp lại phép làm tròn để chọn lát cắt tốt nhất
func (gw *MaxCutGW) ApproximateMaxCut(numRounds int) ([]bool, float64, float64) {
	sdpBound := gw.ComputeSDPValue()

	var bestCut []bool
	bestWeight := -1.0

	for round := 0; round < numRounds; round++ {
		assignment, weight := gw.RandomHyperplaneRound()
		if weight > bestWeight {
			bestWeight = weight
			bestCut = assignment
		}
	}

	return bestCut, bestWeight, sdpBound
}

func main() {
	// Khởi tạo đồ thị chu trình ngũ giác C5 (5 đỉnh: 0-1-2-3-4-0)
	// Đây là ví dụ kinh điển:
	// - Chu trình C5 không phải là đồ thị hai phía.
	// - Lát cắt cực đại nguyên tối ưu là 4 cạnh (bỏ 1 cạnh).
	// - SDP nới lỏng cho giá trị: 5 * (1 - cos(4pi/5)) / 2 = 5 * (1 - (-0.809)) / 2 = 4.522 cạnh.
	n := 5
	edges := []Edge{
		{0, 1, 1.0},
		{1, 2, 1.0},
		{2, 3, 1.0},
		{3, 4, 1.0},
		{4, 0, 1.0},
	}

	gw := NewMaxCutGW(n, edges)

	// 1. Giải cấu hình vector SDP bằng Riemannian Gradient Descent
	gw.SolveVectorEmbedding(300, 0.05)

	// 2. Thực hiện làm tròn siêu phẳng ngẫu nhiên (100 lượt lấy mẫu)
	cutAssignment, cutWeight, sdpBound := gw.ApproximateMaxCut(100)

	fmt.Printf("=== KẾT QUẢ THUẬT TOÁN GOEMANS-WILLIAMSON (MAX-CUT) ===\n")
	fmt.Printf("Số lượng đỉnh: %d | Số lượng cạnh: %d\n", n, len(edges))
	fmt.Printf("Chặn trên từ SDP Relaxation: %.4f\n", sdpBound)
	fmt.Printf("Trọng số lát cắt tìm được: %.1f\n", cutWeight)
	fmt.Printf("Tỷ lệ đạt được thực tế so với SDP: %.4f (Lý thuyết: >= 0.8785)\n\n", cutWeight/sdpBound)

	fmt.Println("Phân hoạch các đỉnh:")
	var setS, setNotS []int
	for i, inS := range cutAssignment {
		if inS {
			setS = append(setS, i)
		} else {
			setNotS = append(setNotS, i)
		}
	}
	fmt.Printf("- Tập S:     %v\n", setS)
	fmt.Printf("- Tập V \\ S: %v\n", setNotS)

	fmt.Println("\nChi tiết các cạnh bị cắt:")
	for _, e := range edges {
		if cutAssignment[e.u] != cutAssignment[e.v] {
			fmt.Printf("- Cạnh (%d <---> %d) [Bị cắt]\n", e.u, e.v)
		}
	}
}

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

6. Ý nghĩa Lý thuyết và Giới hạn Không thể vượt qua (Unique Games Conjecture)

Thuật toán Goemans-Williamson không chỉ là một thuật toán xấp xỉ đơn thuần, mà là một cột mốc xác lập ranh giới của lý thuyết độ phức tạp tính toán:

1. Định lý Håstad (1997 - PCP Theorem): Johan Håstad chứng minh rằng bài toán Max-Cut không thể được xấp xỉ với tỷ lệ tốt hơn:

   $$\alpha \ge \frac{16}{17} \approx 0.941176$$

   trừ khi $P = NP$.

2. Đột phá Khot - Kindler - Mossel - O'Donnell (2007): Giả sử Giả thuyết Trò chơi Độc nhất (Unique Games Conjecture - UGC) của Subhash Khot là đúng, thì việc xấp xỉ Max-Cut với tỷ lệ:

   $$\alpha > \alpha_{\text{GW}} = \min_{\theta} \frac{2}{\pi} \frac{\theta}{1 - \cos \theta} \approx 0.878567$$

   là NP-hard!

Hệ quả: Tỷ số $0.878567$ của Goemans-Williamson là tỷ số xấp xỉ tối ưu toàn cục. Thuật toán GW chính là cấu trúc giải thuật tốt nhất có thể tồn tại trong thời gian đa thức cho bài toán Max-Cut.

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

• Vật lý Thống kê & Mô hình Kính Spin (Ising Spin Glass Models): Trong vật lý vật chất ngưng tụ, việc tìm trạng thái năng lượng đáy (Ground State) của mạng lưới spin tương tác phản sắt từ (Antiferromagnetic Ising Model) tương đương với việc tìm lát cắt cực đại trên đồ thị tương tác spin.

• Tối ưu hóa Bố trí Vi mạch VLSI (Two-Layer Via Minimization): Trong thiết kế bo mạch in và chip bán dẫn, các đường dây mạch giao nhau được phân chia vào 2 lớp kim loại khác nhau để tránh ngắn mạch. Bài toán giảm thiểu số lượng lỗ xuyên lớp (vias) tương đương với bài toán Max-Cut trên đồ thị xung đột.

• Phân cụm Tương phản trong Mạng xã hội (Correlation Clustering on Signed Graphs): Trong mạng lưới có các mối quan hệ tích cực (bạn bè, trọng số dương) và tiêu cực (đối thủ, trọng số âm), Max-Cut được sử dụng để phân cực mạng thành hai liên minh đối kháng nhau với mức độ mâu thuẫn nội bộ thấp nhất.