魯棒性認證:隨機平滑的ℓ₁/ℓ₂認證半徑精確分析

對抗樣本的存在猶如深度學習領域的"暗物質",揭示了模型決策邊界中隱藏的脆弱性。面對這一挑戰,隨機平滑技術脱穎而出,成為首個能夠為大規模深度學習模型提供可證明魯棒性保證的實用方法。然而,傳統的隨機平滑分析主要集中於ℓ₂範數威脅模型,在現實世界中更為常見的ℓ₁威脅(如稀疏對抗擾動)面前顯得力不從心。

本文將深入探討隨機平滑在ℓ₁和ℓ₂威脅模型下的認證半徑精確分析,通過嚴格的數學推導和詳細的代碼實現,揭示不同噪聲分佈對認證半徑的影響,並給出最優噪聲選擇的指導原則。我們不僅會展示如何計算緊緻的認證邊界,還將提供實用的認證框架和優化策略。

隨機平滑理論基礎與認證框架

隨機平滑的基本原理與認證機制

隨機平滑通過向輸入添加隨機噪聲,將任意基礎分類器轉換為平滑分類器,從而獲得可證明的魯棒性保證:

import torch
import torch.nn as nn
import numpy as np
from scipy import special
import matplotlib.pyplot as plt
from torch.distributions import Normal, Laplace, Gamma
import math

class RandomSmoothingCertification:
    """隨機平滑認證基礎框架"""
    
    def __init__(self, base_classifier, num_classes=10):
        self.base_classifier = base_classifier
        self.num_classes = num_classes
        
    def smooth_classifier(self, x, noise_std, num_samples=1000):
        """構建平滑分類器"""
        batch_size = x.shape[0]
        
        # 重複輸入以進行蒙特卡洛採樣
        x_expanded = x.repeat(num_samples, 1, 1, 1)
        noise = torch.randn_like(x_expanded) * noise_std
        
        # 添加噪聲並獲取預測
        with torch.no_grad():
            noisy_outputs = self.base_classifier(x_expanded + noise)
            predictions = torch.argmax(noisy_outputs, dim=1)
            
        # 統計每個類的票數
        votes = torch.zeros(batch_size, self.num_classes, device=x.device)
        for i in range(num_samples):
            votes[torch.arange(batch_size), predictions[i*batch_size:(i+1)*batch_size]] += 1
            
        # 計算概率和預測
        probabilities = votes / num_samples
        predicted_classes = torch.argmax(votes, dim=1)
        
        return probabilities, predicted_classes
    
    def _binary_search_radius(self, pA, alpha, norm_type='l2'):
        """二分搜索認證半徑"""
        if norm_type == 'l2':
            return noise_std * special.erfinv(2 * pA - 1)
        elif norm_type == 'l1':
            # 使用更精確的ℓ₁認證
            return self._l1_certification_radius(pA, alpha)
        else:
            raise ValueError(f"不支持的範數類型: {norm_type}")
    
    def _l1_certification_radius(self, pA, alpha, precision=1e-6):
        """計算ℓ₁認證半徑的精確解"""
        # 通過數值方法求解認證半徑
        def equation(r):
            return self._l1_certification_probability(pA, r, alpha) - 0.5
        
        # 二分搜索求解
        left, right = 0, 10.0  # 初始搜索範圍
        while right - left > precision:
            mid = (left + right) / 2
            if equation(mid) > 0:
                left = mid
            else:
                right = mid
                
        return (left + right) / 2
    
    def _l1_certification_probability(self, pA, radius, alpha):
        """計算ℓ₁認證概率"""
        # 基於Laplace噪聲的認證概率計算
        return pA - 0.5 * (1 - torch.exp(-radius / alpha))

ℓ₁與ℓ₂威脅模型的數學基礎

不同威脅模型下的認證需要不同的數學工具和分析方法:

class ThreatModelAnalysis:
    """ℓ₁和ℓ₂威脅模型的數學分析"""
    
    def __init__(self):
        self.norm_properties = {
            'l1': {'name': 'ℓ₁範數', 'ball_type': '菱形', 'applications': ['稀疏攻擊', '特徵選擇']},
            'l2': {'name': 'ℓ₂範數', 'ball_type': '球形', 'applications': ['一般擾動', '物理攻擊']},
            'linf': {'name': 'ℓ∞範數', 'ball_type': '立方體', 'applications': ['像素級攻擊']}
        }
    
    def norm_ball_geometry(self, radius, norm_type='l2', dim=2):
        """可視化不同範數球的幾何特性"""
        if dim != 2:
            raise ValueError("只支持2維可視化")
            
        theta = np.linspace(0, 2*np.pi, 100)
        
        if norm_type == 'l2':
            # ℓ₂球:圓形
            x = radius * np.cos(theta)
            y = radius * np.sin(theta)
        elif norm_type == 'l1':
            # ℓ₁球:菱形
            x = []
            y = []
            for angle in theta:
                if angle <= np.pi/2:
                    x.append(radius * np.cos(angle))
                    y.append(radius * np.sin(angle))
                elif angle <= np.pi:
                    x.append(-radius * np.cos(np.pi - angle))
                    y.append(radius * np.sin(np.pi - angle))
                elif angle <= 3*np.pi/2:
                    x.append(-radius * np.cos(angle - np.pi))
                    y.append(-radius * np.sin(angle - np.pi))
                else:
                    x.append(radius * np.cos(2*np.pi - angle))
                    y.append(-radius * np.sin(2*np.pi - angle))
        elif norm_type == 'linf':
            # ℓ∞球:正方形
            t_param = np.linspace(0, 1, 25)
            x = np.concatenate([t_param*2*radius - radius, 
                              np.ones(25)*radius,
                              (1-t_param)*2*radius - radius,
                              -np.ones(25)*radius])
            y = np.concatenate([-np.ones(25)*radius,
                              t_param*2*radius - radius,
                              np.ones(25)*radius,
                              (1-t_param)*2*radius - radius])
        
        return x, y
    
    def plot_norm_comparison(self):
        """比較不同範數球的幾何形狀"""
        fig, axes = plt.subplots(1, 3, figsize=(15, 5))
        norms = ['l1', 'l2', 'linf']
        radius = 1.0
        
        for i, norm in enumerate(norms):
            x, y = self.norm_ball_geometry(radius, norm_type=norm)
            axes[i].plot(x, y, 'b-', linewidth=2)
            axes[i].fill(x, y, alpha=0.3)
            axes[i].set_title(f'{self.norm_properties[norm]["name"]}球')
            axes[i].set_xlim(-1.5, 1.5)
            axes[i].set_ylim(-1.5, 1.5)
            axes[i].grid(True, alpha=0.3)
            axes[i].set_aspect('equal')
            
        plt.tight_layout()
        plt.show()
    
    def noise_distribution_analysis(self, noise_type='gaussian', scale=1.0):
        """分析不同噪聲分佈的統計特性"""
        x = np.linspace(-5, 5, 1000)
        
        if noise_type == 'gaussian':
            pdf = np.exp(-0.5 * (x/scale)**2) / (scale * np.sqrt(2*np.pi))
            cdf = 0.5 * (1 + special.erf(x / (scale * np.sqrt(2))))
            name = '高斯分佈'
        elif noise_type == 'laplace':
            pdf = np.exp(-np.abs(x)/scale) / (2 * scale)
            cdf = 0.5 * (1 + np.sign(x) * (1 - np.exp(-np.abs(x)/scale)))
            name = '拉普拉斯分佈'
        else:
            raise ValueError("不支持的噪聲類型")
            
        return x, pdf, cdf, name

# 威脅模型分析
threat_analysis = ThreatModelAnalysis()
threat_analysis.plot_norm_comparison()

# 噪聲分佈分析
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# PDF比較
x_gauss, pdf_gauss, _, name_gauss = threat_analysis.noise_distribution_analysis('gaussian', 1.0)
x_laplace, pdf_laplace, _, name_laplace = threat_analysis.noise_distribution_analysis('laplace', 1.0)

ax1.plot(x_gauss, pdf_gauss, 'b-', label=name_gauss, linewidth=2)
ax1.plot(x_laplace, pdf_laplace, 'r-', label=name_laplace, linewidth=2)
ax1.set_xlabel('x')
ax1.set_ylabel('概率密度')
ax1.set_title('噪聲分佈PDF比較')
ax1.legend()
ax1.grid(True, alpha=0.3)

# CDF比較
_, _, cdf_gauss, _ = threat_analysis.noise_distribution_analysis('gaussian', 1.0)
_, _, cdf_laplace, _ = threat_analysis.noise_distribution_analysis('laplace', 1.0)

ax2.plot(x_gauss, cdf_gauss, 'b-', label=name_gauss, linewidth=2)
ax2.plot(x_laplace, cdf_laplace, 'r-', label=name_laplace, linewidth=2)
ax2.set_xlabel('x')
ax2.set_ylabel('累積分佈')
ax2.set_title('噪聲分佈CDF比較')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

ℓ₂認證半徑的精確分析

高斯隨機平滑的認證理論

對於ℓ₂威脅模型,高斯噪聲提供了最優的認證保證:

class L2Certification:
    """ℓ₂威脅模型的精確認證分析"""
    
    def __init__(self, noise_std=0.5):
        self.noise_std = noise_std
        
    def gaussian_smoothing_certificate(self, pA, alpha=0.001):
        """
        高斯隨機平滑的ℓ₂認證證書
        基於Cohen et al. (2019)的理論
        """
        # 計算認證半徑
        radius = self.noise_std * special.erfinv(2 * pA - 1)
        
        # 計算認證概率下界
        certified_prob = 0.5 * (1 + special.erf(
            (radius - alpha) / (self.noise_std * np.sqrt(2))
        ))
        
        return {
            'certified_radius': radius,
            'certified_probability': certified_prob,
            'noise_std': self.noise_std,
            'confidence_level': pA
        }
    
    def tight_certification_bound(self, pA, delta=1e-5):
        """
        計算緊緻的ℓ₂認證邊界
        基於更精確的概率分析
        """
        def objective_function(r):
            # 基於NIPS 2019的緊緻邊界
            term1 = special.ndtr(special.ndtri(pA) - r / self.noise_std)
            term2 = special.ndtr(-special.ndtri(pA) - r / self.noise_std)
            return term1 + term2 - 1 + delta
        
        # 二分搜索找到最大認證半徑
        left, right = 0, 5 * self.noise_std
        tolerance = 1e-6
        
        while right - left > tolerance:
            mid = (left + right) / 2
            if objective_function(mid) > 0:
                left = mid
            else:
                right = mid
                
        certified_radius = (left + right) / 2
        
        return certified_radius
    
    def optimal_noise_analysis(self, pA_range=np.linspace(0.6, 0.99, 50)):
        """分析不同噪聲水平對認證半徑的影響"""
        noise_stds = [0.1, 0.25, 0.5, 1.0]
        
        plt.figure(figsize=(10, 6))
        
        for noise_std in noise_stds:
            self.noise_std = noise_std
            radii = []
            
            for pA in pA_range:
                radius = self.gaussian_smoothing_certificate(pA)['certified_radius']
                radii.append(radius)
                
            plt.plot(pA_range, radii, 
                    label=f'σ={noise_std}', linewidth=2)
        
        plt.xlabel('基礎分類器置信度 pA')
        plt.ylabel('認證半徑')
        plt.title('ℓ₂認證半徑 vs 置信度 (不同噪聲水平)')
        plt.legend()
        plt.grid(True, alpha=0.3)
        plt.show()
    
    def multidimensional_certification(self, pA, dimensions):
        """高維空間中的認證半徑分析"""
        # 高維空間中的認證半徑縮放
        base_radius = self.gaussian_smoothing_certificate(pA)['certified_radius']
        
        # 考慮維度影響的縮放因子
        scaling_factors = []
        for d in dimensions:
            # 高維空間中的概率集中現象
            scaling_factor = 1.0 / np.sqrt(d)
            scaling_factors.append(scaling_factor)
            
        scaled_radii = [base_radius * s for s in scaling_factors]
        
        return {
            'dimensions': dimensions,
            'base_radius': base_radius,
            'scaling_factors': scaling_factors,
            'scaled_radii': scaled_radii
        }

# ℓ₂認證分析
l2_cert = L2Certification(noise_std=0.5)

# 基本認證計算
pA = 0.8
cert_result = l2_cert.gaussian_smoothing_certificate(pA)
print(f"ℓ₂認證結果 (pA={pA}):")
print(f"  認證半徑: {cert_result['certified_radius']:.4f}")
print(f"  認證概率: {cert_result['certified_probability']:.4f}")

# 緊緻邊界計算
tight_radius = l2_cert.tight_certification_bound(pA)
print(f"緊緻認證半徑: {tight_radius:.4f}")

# 可視化分析
l2_cert.optimal_noise_analysis()

# 高維分析
dimensions = [10, 100, 1000, 10000]
dim_analysis = l2_cert.multidimensional_certification(pA, dimensions)

plt.figure(figsize=(10, 6))
plt.plot(dim_analysis['dimensions'], dim_analysis['scaled_radii'], 'bo-', linewidth=2)
plt.xscale('log')
plt.xlabel('維度')
plt.ylabel('縮放後的認證半徑')
plt.title('ℓ₂認證半徑的維度縮放')
plt.grid(True, alpha=0.3)
plt.show()

基於Neyman-Pearson引理的緊緻分析

通過統計假設檢驗的理論可以獲得更緊緻的認證邊界:

class NeymanPearsonCertification:
    """基於Neyman-Pearson引理的緊緻認證分析"""
    
    def __init__(self, noise_std=0.5):
        self.noise_std = noise_std
        
    def likelihood_ratio_test(self, x, x_perturbed, norm_type='l2'):
        """計算似然比檢驗統計量"""
        if norm_type == 'l2':
            # 對於高斯噪聲,ℓ₂擾動
            noise_diff = x_perturbed - x
            likelihood_ratio = torch.exp(
                -0.5 * torch.norm(noise_diff, p=2)**2 / self.noise_std**2
            )
        elif norm_type == 'l1':
            # 對於拉普拉斯噪聲,ℓ₁擾動
            noise_diff = x_perturbed - x
            likelihood_ratio = torch.exp(
                -torch.norm(noise_diff, p=1) / self.noise_std
            )
        else:
            raise ValueError("不支持的範數類型")
            
        return likelihood_ratio
    
    def neyman_pearson_bound(self, pA, epsilon, norm_type='l2'):
        """基於Neyman-Pearson引理的認證邊界"""
        
        if norm_type == 'l2':
            # 高斯噪聲的NP引理應用
            # 認證半徑 R = σ * Φ^{-1}(pA)
            radius = self.noise_std * special.ndtri(pA)
            
            # 在擾動ε下的最小概率
            min_prob = special.ndtr(special.ndtri(pA) - epsilon / self.noise_std)
            
        elif norm_type == 'l1':
            # 拉普拉斯噪聲的NP引理應用
            # 需要數值求解
            def probability_bound(r):
                return pA - 0.5 * (1 - np.exp(-r / self.noise_std))
            
            # 找到滿足概率下界的最大半徑
            left, right = 0, 10 * self.noise_std
            tolerance = 1e-6
            
            while right - left > tolerance:
                mid = (left + right) / 2
                if probability_bound(mid) > 0.5:
                    left = mid
                else:
                    right = mid
                    
            radius = (left + right) / 2
            min_prob = probability_bound(epsilon)
            
        return {
            'certified_radius': radius,
            'min_probability_at_epsilon': min_prob,
            'norm_type': norm_type
        }
    
    def plot_neyman_pearson_curves(self, pA_values=[0.6, 0.7, 0.8, 0.9]):
        """繪製Neyman-Pearson認證曲線"""
        
        epsilon_range = np.linspace(0, 2.0, 100)
        
        plt.figure(figsize=(12, 8))
        
        for pA in pA_values:
            min_probs = []
            for epsilon in epsilon_range:
                result = self.neyman_pearson_bound(pA, epsilon, 'l2')
                min_probs.append(result['min_probability_at_epsilon'])
                
            plt.plot(epsilon_range, min_probs, 
                    label=f'pA = {pA}', linewidth=2)
            
            # 標記認證半徑
            cert_radius = self.neyman_pearson_bound(pA, 0, 'l2')['certified_radius']
            plt.axvline(x=cert_radius, linestyle='--', alpha=0.7,
                       color=plt.gca().lines[-1].get_color())
        
        plt.axhline(y=0.5, color='red', linestyle='-', 
                   label='決策邊界 (0.5)', alpha=0.7)
        plt.xlabel('擾動大小 ε')
        plt.ylabel('最小概率下界')
        plt.title('Neyman-Pearson認證邊界')
        plt.legend()
        plt.grid(True, alpha=0.3)
        plt.show()

# Neyman-Pearson認證分析
np_cert = NeymanPearsonCertification(noise_std=0.5)
np_cert.plot_neyman_pearson_curves()

# 比較不同方法的認證半徑
pA = 0.8
np_result = np_cert.neyman_pearson_bound(pA, 0, 'l2')
basic_result = l2_cert.gaussian_smoothing_certificate(pA)

print("認證半徑比較:")
print(f"Neyman-Pearson方法: {np_result['certified_radius']:.4f}")
print(f"基礎高斯方法: {basic_result['certified_radius']:.4f}")

ℓ₁認證半徑的精確分析

拉普拉斯隨機平滑的理論基礎

對於ℓ₁威脅模型,拉普拉斯噪聲提供了更合適的認證基礎:

class L1Certification:
    """ℓ₁威脅模型的精確認證分析"""
    
    def __init__(self, noise_scale=0.5):
        self.noise_scale = noise_scale  # 拉普拉斯分佈的尺度參數
        
    def laplace_smoothing_certificate(self, pA, method='exact'):
        """
        拉普拉斯隨機平滑的ℓ₁認證證書
        基於Lecuyer et al. (2019)的理論
        """
        if method == 'exact':
            # 精確認證半徑
            if pA <= 0.5:
                certified_radius = 0.0
            else:
                certified_radius = self.noise_scale * np.log(2 * pA)
                
        elif method == 'improved':
            # 改進的認證邊界
            certified_radius = self._improved_l1_certificate(pA)
        else:
            raise ValueError("不支持的認證方法")
            
        # 計算認證概率
        certified_prob = self._l1_certification_probability(pA, certified_radius)
        
        return {
            'certified_radius': certified_radius,
            'certified_probability': certified_prob,
            'noise_scale': self.noise_scale,
            'confidence_level': pA
        }
    
    def _improved_l1_certificate(self, pA):
        """改進的ℓ₁認證邊界"""
        # 基於更精細的概率分析
        if pA <= 0.5:
            return 0.0
        
        # 數值求解認證半徑
        def probability_equation(r):
            return pA - 0.5 * np.exp(-r / self.noise_scale) - 0.5
        
        # 二分搜索
        left, right = 0, 10 * self.noise_scale
        tolerance = 1e-6
        
        while right - left > tolerance:
            mid = (left + right) / 2
            if probability_equation(mid) > 0:
                left = mid
            else:
                right = mid
                
        return (left + right) / 2
    
    def _l1_certification_probability(self, pA, radius):
        """計算ℓ₁認證概率"""
        if radius == 0:
            return pA
        
        # 基於拉普拉斯分佈的認證概率
        prob = pA - 0.5 * (1 - np.exp(-radius / self.noise_scale))
        return max(prob, 0.5)  # 確保不低於隨機猜測
    
    def sparse_attack_certification(self, pA, sparsity_level):
        """
        針對稀疏攻擊的專用認證
        sparsity_level: 攻擊者可以修改的最大特徵數
        """
        # 基礎認證半徑
        base_radius = self.laplace_smoothing_certificate(pA)['certified_radius']
        
        # 稀疏攻擊下的有效半徑
        # 在稀疏約束下,攻擊者可以集中擾動在少數特徵上
        effective_radius = base_radius / np.sqrt(sparsity_level)
        
        return {
            'base_radius': base_radius,
            'sparsity_level': sparsity_level,
            'effective_radius': effective_radius,
            'certification_strength': '強' if effective_radius > base_radius * 0.5 else '弱'
        }
    
    def compare_noise_distributions(self, pA_range=np.linspace(0.55, 0.95, 50)):
        """比較不同噪聲分佈在ℓ₁認證中的效果"""
        
        noise_scales = [0.1, 0.25, 0.5, 1.0]
        
        plt.figure(figsize=(12, 8))
        
        for scale in noise_scales:
            self.noise_scale = scale
            radii_exact = []
            radii_improved = []
            
            for pA in pA_range:
                cert_exact = self.laplace_smoothing_certificate(pA, 'exact')
                cert_improved = self.laplace_smoothing_certificate(pA, 'improved')
                
                radii_exact.append(cert_exact['certified_radius'])
                radii_improved.append(cert_improved['certified_radius'])
            
            plt.plot(pA_range, radii_exact, '--', 
                    label=f'精確 (b={scale})', linewidth=2)
            plt.plot(pA_range, radii_improved, '-',
                    label=f'改進 (b={scale})', linewidth=2)
        
        plt.xlabel('基礎分類器置信度 pA')
        plt.ylabel('ℓ₁認證半徑')
        plt.title('拉普拉斯隨機平滑的ℓ₁認證半徑比較')
        plt.legend()
        plt.grid(True, alpha=0.3)
        plt.show()

# ℓ₁認證分析
l1_cert = L1Certification(noise_scale=0.5)

# 基本認證計算
pA = 0.8
l1_result = l1_cert.laplace_smoothing_certificate(pA, 'exact')
print(f"ℓ₁認證結果 (pA={pA}):")
print(f"  認證半徑: {l1_result['certified_radius']:.4f}")
print(f"  認證概率: {l1_result['certified_probability']:.4f}")

# 改進認證
l1_improved = l1_cert.laplace_smoothing_certificate(pA, 'improved')
print(f"改進ℓ₁認證半徑: {l1_improved['certified_radius']:.4f}")

# 稀疏攻擊認證
sparse_result = l1_cert.sparse_attack_certification(pA, sparsity_level=10)
print(f"稀疏攻擊認證:")
print(f"  基礎半徑: {sparse_result['base_radius']:.4f}")
print(f"  有效半徑: {sparse_result['effective_radius']:.4f}")
print(f"  認證強度: {sparse_result['certification_strength']}")

# 比較不同噪聲分佈
l1_cert.compare_noise_distributions()

ℓ₁認證的緊緻下界與優化

通過更精細的概率分析,我們可以獲得ℓ₁認證的緊緻下界:

class TightL1Certification:
    """ℓ₁認證的緊緻下界分析"""
    
    def __init__(self, noise_scale=0.5):
        self.noise_scale = noise_scale
        
    def compute_tight_bound(self, pA, delta=1e-5):
        """計算ℓ₁認證的緊緻下界"""
        
        def certification_condition(r):
            """
            認證條件: pA - 0.5 * (1 - exp(-r/b)) >= 0.5 + delta
            這確保了在擾動r下,正確類的概率仍然超過0.5
            """
            return pA - 0.5 * (1 - np.exp(-r / self.noise_scale)) - 0.5 - delta
        
        if certification_condition(0) < 0:
            return 0.0
        
        # 二分搜索找到最大認證半徑
        left, right = 0, 10 * self.noise_scale
        tolerance = 1e-6
        
        while right - left > tolerance:
            mid = (left + right) / 2
            if certification_condition(mid) > 0:
                left = mid
            else:
                right = mid
                
        return (left + right) / 2
    
    def multidimensional_l1_certification(self, pA, dimensions, attack_budget):
        """
        高維空間中的ℓ₁認證分析
        attack_budget: 攻擊者的總ℓ₁預算
        """
        # 基礎認證半徑
        base_radius = self.compute_tight_bound(pA)
        
        # 高維空間中的認證半徑縮放
        # 在ℓ₁約束下,維度的影響與ℓ₂不同
        dimensional_scaling = 1.0 / np.sqrt(dimensions)
        
        scaled_radii = base_radius * dimensional_scaling
        
        # 考慮攻擊預算的有效認證
        effective_radii = np.minimum(scaled_radii, attack_budget)
        
        return {
            'dimensions': dimensions,
            'base_radius': base_radius,
            'scaled_radii': scaled_radii,
            'effective_radii': effective_radii,
            'attack_budget': attack_budget
        }
    
    def plot_tight_bounds_comparison(self):
        """繪製緊緻邊界與基礎邊界的比較"""
        pA_range = np.linspace(0.55, 0.95, 100)
        
        basic_radii = []
        tight_radii = []
        
        for pA in pA_range:
            # 基礎邊界
            basic_radius = self.noise_scale * np.log(2 * pA) if pA > 0.5 else 0.0
            basic_radii.append(basic_radius)
            
            # 緊緻邊界
            tight_radius = self.compute_tight_bound(pA)
            tight_radii.append(tight_radius)
        
        plt.figure(figsize=(10, 6))
        plt.plot(pA_range, basic_radii, 'r--', label='基礎邊界', linewidth=2)
        plt.plot(pA_range, tight_radii, 'b-', label='緊緻邊界', linewidth=2)
        
        plt.xlabel('基礎分類器置信度 pA')
        plt.ylabel('ℓ₁認證半徑')
        plt.title('ℓ₁認證的緊緻邊界 vs 基礎邊界')
        plt.legend()
        plt.grid(True, alpha=0.3)
        plt.show()
        
        # 計算改進比例
        improvement_ratio = np.array(tight_radii) / np.array(basic_radii)
        avg_improvement = np.mean(improvement_ratio[np.isfinite(improvement_ratio)])
        
        print(f"緊緻邊界平均改進比例: {avg_improvement:.4f}")
        
        return avg_improvement

# 緊緻ℓ₁認證分析
tight_l1 = TightL1Certification(noise_scale=0.5)

# 計算緊緻邊界
pA = 0.8
tight_radius = tight_l1.compute_tight_bound(pA)
print(f"緊緻ℓ₁認證半徑: {tight_radius:.4f}")

# 高維分析
dimensions = [10, 100, 1000, 10000]
dim_analysis = tight_l1.multidimensional_l1_certification(pA, dimensions, attack_budget=2.0)

plt.figure(figsize=(10, 6))
plt.plot(dim_analysis['dimensions'], dim_analysis['scaled_radii'], 'bo-', 
         label='縮放半徑', linewidth=2)
plt.plot(dim_analysis['dimensions'], dim_analysis['effective_radii'], 'ro-',
         label='有效半徑', linewidth=2)
plt.axhline(y=dim_analysis['attack_budget'], color='green', linestyle='--',
           label='攻擊預算', alpha=0.7)
plt.xscale('log')
plt.xlabel('維度')
plt.ylabel('認證半徑')
plt.title('高維空間中的ℓ₁認證半徑')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# 邊界比較
improvement = tight_l1.plot_tight_bounds_comparison()

混合範數認證與最優噪聲選擇

ℓ₁/ℓ₂混合認證框架

在實際應用中,我們經常需要同時考慮多種威脅模型:

class HybridCertification:
    """ℓ₁/ℓ₂混合認證框架"""
    
    def __init__(self, l1_scale=0.5, l2_std=0.5):
        self.l1_scale = l1_scale
        self.l2_std = l2_std
        
    def hybrid_noise_certification(self, pA, threat_model='both'):
        """
        混合噪聲下的認證分析
        threat_model: 'l1', 'l2', 或 'both'
        """
        if threat_model == 'l1':
            # 純ℓ₁認證
            l1_cert = L1Certification(self.l1_scale)
            return l1_cert.laplace_smoothing_certificate(pA, 'improved')
            
        elif threat_model == 'l2':
            # 純ℓ₂認證
            l2_cert = L2Certification(self.l2_std)
            return l2_cert.gaussian_smoothing_certificate(pA)
            
        elif threat_model == 'both':
            # 混合認證 - 取較弱的保證
            l1_cert = L1Certification(self.l1_scale)
            l2_cert = L2Certification(self.l2_std)
            
            l1_radius = l1_cert.laplace_smoothing_certificate(pA, 'improved')['certified_radius']
            l2_radius = l2_cert.gaussian_smoothing_certificate(pA)['certified_radius']
            
            # 對於混合威脅,我們取最小認證半徑作為保守估計
            certified_radius = min(l1_radius, l2_radius)
            
            return {
                'certified_radius': certified_radius,
                'l1_radius': l1_radius,
                'l2_radius': l2_radius,
                'threat_model': 'mixed_l1_l2'
            }
    
    def optimal_noise_selection(self, pA, dimension, threat_weights=None):
        """
        最優噪聲參數選擇
        threat_weights: [w_l1, w_l2] 威脅模型權重
        """
        if threat_weights is None:
            threat_weights = [0.5, 0.5]  # 默認均衡權重
            
        # 網格搜索最優噪聲參數
        scale_range = np.linspace(0.1, 2.0, 20)
        std_range = np.linspace(0.1, 2.0, 20)
        
        best_score = -np.inf
        best_params = None
        results = []
        
        for scale in scale_range:
            for std in std_range:
                self.l1_scale = scale
                self.l2_std = std
                
                cert_result = self.hybrid_noise_certification(pA, 'both')
                
                if 'l1_radius' in cert_result and 'l2_radius' in cert_result:
                    # 加權認證分數
                    score = (threat_weights[0] * cert_result['l1_radius'] + 
                            threat_weights[1] * cert_result['l2_radius'])
                    
                    results.append((scale, std, score, 
                                  cert_result['l1_radius'], cert_result['l2_radius']))
                    
                    if score > best_score:
                        best_score = score
                        best_params = (scale, std, cert_result['l1_radius'], cert_result['l2_radius'])
        
        return {
            'optimal_scale': best_params[0],
            'optimal_std': best_params[1],
            'optimal_l1_radius': best_params[2],
            'optimal_l2_radius': best_params[3],
            'best_score': best_score,
            'all_results': results
        }
    
    def plot_optimal_noise_landscape(self, pA=0.8, threat_weights=[0.5, 0.5]):
        """繪製最優噪聲參數的可視化"""
        
        result = self.optimal_noise_selection(pA, 100, threat_weights)
        all_results = result['all_results']
        
        # 提取數據
        scales = [r[0] for r in all_results]
        stds = [r[1] for r in all_results]
        scores = [r[2] for r in all_results]
        
        # 創建網格數據
        unique_scales = sorted(set(scales))
        unique_stds = sorted(set(stds))
        
        score_grid = np.zeros((len(unique_scales), len(unique_stds)))
        
        for i, scale in enumerate(unique_scales):
            for j, std in enumerate(unique_stds):
                for r in all_results:
                    if abs(r[0] - scale) < 1e-6 and abs(r[1] - std) < 1e-6:
                        score_grid[i, j] = r[2]
                        break
        
        # 繪製熱力圖
        plt.figure(figsize=(10, 8))
        im = plt.imshow(score_grid, extent=[min(unique_stds), max(unique_stds), 
                                          min(unique_scales), max(unique_scales)], 
                       origin='lower', aspect='auto', cmap='viridis')
        
        # 標記最優點
        plt.plot(result['optimal_std'], result['optimal_scale'], 'ro', 
                markersize=10, label='最優點')
        
        plt.colorbar(im, label='加權認證分數')
        plt.xlabel('高斯噪聲標準差 σ')
        plt.ylabel('拉普拉斯噪聲尺度 b')
        plt.title('最優噪聲參數選擇')
        plt.legend()
        plt.show()
        
        return result

# 混合認證分析
hybrid_cert = HybridCertification()

# 混合認證計算
pA = 0.8
mixed_result = hybrid_cert.hybrid_noise_certification(pA, 'both')
print(f"混合認證結果 (pA={pA}):")
print(f"  ℓ₁認證半徑: {mixed_result['l1_radius']:.4f}")
print(f"  ℓ₂認證半徑: {mixed_result['l2_radius']:.4f}")
print(f"  最終認證半徑: {mixed_result['certified_radius']:.4f}")

# 最優噪聲選擇
optimal_result = hybrid_cert.optimal_noise_selection(pA, 100, [0.5, 0.5])
print(f"最優噪聲參數:")
print(f"  拉普拉斯尺度: {optimal_result['optimal_scale']:.4f}")
print(f"  高斯標準差: {optimal_result['optimal_std']:.4f}")
print(f"  最優ℓ₁半徑: {optimal_result['optimal_l1_radius']:.4f}")
print(f"  最優ℓ₂半徑: {optimal_result['optimal_l2_radius']:.4f}")

# 可視化噪聲參數景觀
landscape_result = hybrid_cert.plot_optimal_noise_landscape()

實際應用與性能評估

完整認證流程實現

下面實現一個完整的隨機平滑認證流程:

class CompleteCertificationPipeline:
    """完整的隨機平滑認證流程"""
    
    def __init__(self, base_classifier, num_classes=10):
        self.base_classifier = base_classifier
        self.num_classes = num_classes
        self.certification_history = []
        
    def monte_carlo_certification(self, x, y_true, noise_params, 
                                num_samples=10000, alpha=0.001):
        """
        蒙特卡洛認證流程
        基於Cohen et al. (2019)的算法
        """
        batch_size = x.shape[0]
        
        # 步驟1: 估計top-1和top-2類的概率
        pA, pB = self._estimate_class_probabilities(
            x, noise_params, num_samples, alpha)
        
        # 步驟2: 計算認證半徑
        certification_results = []
        for i in range(batch_size):
            if pA[i] > 0.5:  # 只有置信度超過0.5才能認證
                if noise_params['type'] == 'gaussian':
                    certifier = L2Certification(noise_params['std'])
                    radius = certifier.gaussian_smoothing_certificate(
                        pA[i], alpha)['certified_radius']
                elif noise_params['type'] == 'laplace':
                    certifier = L1Certification(noise_params['scale'])
                    radius = certifier.laplace_smoothing_certificate(
                        pA[i], 'improved')['certified_radius']
                else:
                    raise ValueError("不支持的噪聲類型")
                
                is_certified = radius > 0
                correct_prediction = torch.argmax(pA[i]) == y_true[i]
            else:
                radius = 0.0
                is_certified = False
                correct_prediction = False
            
            certification_results.append({
                'radius': radius,
                'is_certified': is_certified,
                'pA': pA[i],
                'pB': pB[i],
                'correct_prediction': correct_prediction,
                'noise_type': noise_params['type']
            })
        
        self.certification_history.extend(certification_results)
        return certification_results
    
    def _estimate_class_probabilities(self, x, noise_params, num_samples, alpha):
        """估計類別概率"""
        batch_size = x.shape[0]
        
        # 生成噪聲樣本
        x_expanded = x.repeat(num_samples, 1, 1, 1)
        
        if noise_params['type'] == 'gaussian':
            noise = torch.randn_like(x_expanded) * noise_params['std']
        elif noise_params['type'] == 'laplace':
            noise = torch.distributions.Laplace(
                0, noise_params['scale']).sample(x_expanded.shape)
        else:
            raise ValueError("不支持的噪聲類型")
        
        # 獲取預測
        with torch.no_grad():
            noisy_outputs = self.base_classifier(x_expanded + noise)
            predictions = torch.argmax(noisy_outputs, dim=1)
        
        # 統計概率
        pA = torch.zeros(batch_size, device=x.device)  # top-1概率
        pB = torch.zeros(batch_size, device=x.device)  # top-2概率
        
        for i in range(batch_size):
            batch_predictions = predictions[i::batch_size]
            counts = torch.bincount(batch_predictions, minlength=self.num_classes)
            sorted_counts, _ = torch.sort(counts, descending=True)
            
            pA[i] = sorted_counts[0].float() / num_samples
            pB[i] = sorted_counts[1].float() / num_samples
        
        # 應用置信區間修正
        pA_lower = self._clopper_pearson_bound(pA, num_samples, alpha/2)
        pB_upper = self._clopper_pearson_bound(pB, num_samples, 1-alpha/2)
        
        return pA_lower, pB_upper
    
    def _clopper_pearson_bound(self, p, n, alpha):
        """Clopper-Pearson置信區間"""
        # 簡化實現 - 實際應使用精確的Beta分佈
        z = torch.tensor(special.ndtri(1 - alpha))
        se = torch.sqrt(p * (1 - p) / n)
        bound = p - z * se
        return torch.clamp(bound, 0.0, 1.0)
    
    def certification_statistics(self):
        """計算認證統計信息"""
        if not self.certification_history:
            return {}
        
        radii = [r['radius'] for r in self.certification_history]
        certified = [r['is_certified'] for r in self.certification_history]
        correct = [r['correct_prediction'] for r in self.certification_history]
        
        certification_rate = np.mean(certified)
        accuracy = np.mean(correct)
        avg_radius = np.mean(radii)
        median_radius = np.median(radii)
        
        # 認證樣本的平均半徑
        certified_radii = [r for r, c in zip(radii, certified) if c]
        avg_certified_radius = np.mean(certified_radii) if certified_radii else 0
        
        return {
            'certification_rate': certification_rate,
            'accuracy': accuracy,
            'average_radius': avg_radius,
            'median_radius': median_radius,
            'average_certified_radius': avg_certified_radius,
            'total_samples': len(self.certification_history)
        }

# 模擬認證流程
class MockClassifier(nn.Module):
    """模擬分類器用於演示"""
    def __init__(self, num_classes=10):
        super().__init__()
        self.num_classes = num_classes
        
    def forward(self, x):
        # 返回隨機logits
        return torch.randn(x.shape[0], self.num_classes)

# 演示認證流程
mock_classifier = MockClassifier()
cert_pipeline = CompleteCertificationPipeline(mock_classifier)

# 模擬數據
batch_size = 10
x_test = torch.randn(batch_size, 3, 32, 32)
y_test = torch.randint(0, 10, (batch_size,))

# 高斯噪聲認證
gaussian_params = {'type': 'gaussian', 'std': 0.5}
gaussian_results = cert_pipeline.monte_carlo_certification(
    x_test, y_test, gaussian_params, num_samples=1000)

# 拉普拉斯噪聲認證  
laplace_params = {'type': 'laplace', 'scale': 0.5}
laplace_results = cert_pipeline.monte_carlo_certification(
    x_test, y_test, laplace_params, num_samples=1000)

# 統計結果
stats = cert_pipeline.certification_statistics()
print("認證統計結果:")
for key, value in stats.items():
    print(f"  {key}: {value:.4f}")

# 可視化認證結果
def plot_certification_results(gaussian_results, laplace_results):
    """可視化認證結果比較"""
    
    gaussian_radii = [r['radius'] for r in gaussian_results if r['is_certified']]
    laplace_radii = [r['radius'] for r in laplace_results if r['is_certified']]
    
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
    
    # 認證半徑分佈
    ax1.hist(gaussian_radii, alpha=0.7, label='高斯噪聲', bins=20)
    ax1.hist(laplace_radii, alpha=0.7, label='拉普拉斯噪聲', bins=20)
    ax1.set_xlabel('認證半徑')
    ax1.set_ylabel('頻數')
    ax1.set_title('認證半徑分佈')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 認證率比較
    gaussian_cert_rate = np.mean([r['is_certified'] for r in gaussian_results])
    laplace_cert_rate = np.mean([r['is_certified'] for r in laplace_results])
    
    ax2.bar(['高斯', '拉普拉斯'], [gaussian_cert_rate, laplace_cert_rate],
           color=['blue', 'orange'], alpha=0.7)
    ax2.set_ylabel('認證率')
    ax2.set_title('噪聲類型對認證率的影響')
    ax2.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()

plot_certification_results(gaussian_results, laplace_results)

理論擴展與未來方向

現有方法的侷限性及改進方向

class TheoreticalExtensions:
    """隨機平滑認證的理論擴展"""
    
    def __init__(self):
        self.current_limitations = [
            "保守的概率邊界",
            "高維度的認證半徑衰減", 
            "對分佈假設的敏感性",
            "計算複雜度高"
        ]
        
        self.research_directions = [
            "更緊緻的認證邊界",
            "自適應噪聲分佈",
            "黑盒認證方法",
            "認證與訓練的聯合優化"
        ]
    
    def advanced_certification_techniques(self):
        """高級認證技術概覽"""
        
        techniques = {
            'denoised_smoothing': {
                'name': '去噪平滑',
                'description': '在添加噪聲前先進行去噪',
                'improvement': '提升認證半徑20-30%',
                'complexity': '中等'
            },
            'randomized_smoothing': {
                'name': '隨機化平滑', 
                'description': '使用隨機化認證框架',
                'improvement': '更緊緻的邊界',
                'complexity': '高'
            },
            'distributionally_robust': {
                'name': '分佈魯棒認證',
                'description': '考慮噪聲分佈的不確定性',
                'improvement': '更好的分佈外泛化',
                'complexity': '高'
            }
        }
        
        print("高級認證技術:")
        for tech_id, tech_info in techniques.items():
            print(f"\n{tech_info['name']}:")
            print(f"  描述: {tech_info['description']}")
            print(f"  改進: {tech_info['improvement']}")
            print(f"  複雜度: {tech_info['complexity']}")
            
        return techniques
    
    def future_research_agenda(self):
        """未來研究議程"""
        
        agenda = {
            'short_term': [
                '混合噪聲分佈的認證理論',
                '針對特定架構的認證優化',
                '認證效率的工程改進'
            ],
            'mid_term': [
                '認證感知的模型訓練',
                '多模態威脅模型的統一認證',
                '認證強度的可解釋性'
            ],
            'long_term': [
                '完全可認證的深度學習框架',
                '認證與泛化的理論統一',
                '自適應認證系統'
            ]
        }
        
        print("\n未來研究議程:")
        for timeframe, topics in agenda.items():
            print(f"\n{timeframe.replace('_', ' ').title()}:")
            for topic in topics:
                print(f"  • {topic}")
                
        return agenda

# 理論擴展分析
extensions = TheoreticalExtensions()
techniques = extensions.advanced_certification_techniques()
agenda = extensions.future_research_agenda()

# 繪製研究路線圖
fig, ax = plt.subplots(figsize=(10, 6))
timeframes = ['短期', '中期', '長期']
topic_counts = [len(agenda['short_term']), len(agenda['mid_term']), len(agenda['long_term'])]

bars = ax.bar(timeframes, topic_counts, color=['lightblue', 'lightgreen', 'lightcoral'])
ax.set_ylabel('研究方向數量')
ax.set_title('隨機平滑認證研究路線圖')
ax.grid(True, alpha=0.3)

# 添加標籤
for bar, count in zip(bars, topic_counts):
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.1,
           f'{count}個方向', ha='center', va='bottom')

plt.tight_layout()
plt.show()

結論

隨機平滑為深度學習模型提供了一種實用且可證明的魯棒性認證框架。通過本文的詳細分析,我們可以得出以下關鍵結論: