核心理论 - 机器学习理论的深层解析¶

import numpy as np
from collections.abc import Callable
import matplotlib.pyplot as plt

class PACLearningFramework:
    """
    PAC学习框架的实现与验证
    """

    def __init__(self, hypothesis_space_size: int, epsilon: float = 0.1, delta: float = 0.05):  # __init__构造方法，创建对象时自动调用
        """
        初始化PAC学习框架

        Args:
            hypothesis_space_size: 假设空间大小 |H|
            epsilon: 精度参数
            delta: 置信度参数
        """
        self.hypothesis_space_size = hypothesis_space_size
        self.epsilon = epsilon
        self.delta = delta

    def sample_complexity_finite(self) -> int:
        """
        计算有限假设空间的样本复杂度

        Returns:
            所需样本数 m
        """
        # m >= (1/epsilon) * (ln|H| + ln(1/delta))
        m = int(np.ceil((1.0 / self.epsilon) *
                       (np.log(self.hypothesis_space_size) + np.log(1.0 / self.delta))))
        return m

    def sample_complexity_vc(self, vc_dim: int) -> int:
        """
        基于VC维的样本复杂度

        Args:
            vc_dim: VC维

        Returns:
            所需样本数 m
        """
        # m >= (1/epsilon) * (8*VC(H)*log(13/epsilon) + 4*log(2/delta))
        m = int(np.ceil((1.0 / self.epsilon) *
                       (8 * vc_dim * np.log(13.0 / self.epsilon) +
                        4 * np.log(2.0 / self.delta))))
        return m

    def generalization_bound_vc(self, m: int, vc_dim: int) -> float:
        """
        基于VC维的泛化误差上界

        Args:
            m: 样本数
            vc_dim: VC维

        Returns:
            泛化误差上界
        """
        # R(h) <= R_emp(h) + sqrt((2d*ln(em/d) + 2*ln(4/delta))/m)
        bound = np.sqrt((2 * vc_dim * np.log(np.e * m / vc_dim) +
                        2 * np.log(4.0 / self.delta)) / m)
        return bound

    def probably_approximately_correct(self,
                                       empirical_risk: float,
                                       m: int,
                                       vc_dim: int) -> tuple[bool, float]:
        """
        验证PAC条件

        Args:
            empirical_risk: 经验风险
            m: 样本数
            vc_dim: VC维

        Returns:
            (是否满足PAC, 泛化误差上界)
        """
        generalization_gap = self.generalization_bound_vc(m, vc_dim)
        true_risk_upper = empirical_risk + generalization_gap

        # 以至少 1-delta 的概率，真实误差不超过上界
        satisfies_pac = generalization_gap <= self.epsilon

        return satisfies_pac, true_risk_upper

    def visualize_sample_complexity(self, vc_dims: list[int] = None):
        """可视化样本复杂度随VC维的变化"""
        if vc_dims is None:
            vc_dims = list(range(1, 101))

        sample_complexities = [self.sample_complexity_vc(d) for d in vc_dims]  # 列表推导式，简洁创建列表

        plt.figure(figsize=(10, 6))
        plt.plot(vc_dims, sample_complexities, 'b-', linewidth=2)
        plt.xlabel('VC Dimension', fontsize=12)
        plt.ylabel('Sample Complexity m', fontsize=12)
        plt.title(f'PAC Sample Complexity (ε={self.epsilon}, δ={self.delta})', fontsize=14)
        plt.grid(True, alpha=0.3)
        plt.tight_layout()
        plt.savefig('pac_sample_complexity.png', dpi=150)
        plt.show()

# 示例：PAC学习框架验证
def pac_learning_example():
    """PAC学习框架示例"""
    # 创建PAC框架实例
    pac = PACLearningFramework(hypothesis_space_size=1000, epsilon=0.1, delta=0.05)

    # 计算不同VC维下的样本复杂度
    vc_dims = [5, 10, 50, 100]
    print("=" * 60)
    print("PAC学习框架分析")
    print("=" * 60)

    print(f"\n假设空间大小 |H| = {pac.hypothesis_space_size}")
    print(f"精度参数 ε = {pac.epsilon}")
    print(f"置信度参数 δ = {pac.delta}")

    print(f"\n有限假设空间样本复杂度: m >= {pac.sample_complexity_finite()}")

    print("\n不同VC维下的样本复杂度:")
    for d in vc_dims:
        m = pac.sample_complexity_vc(d)
        print(f"  VC维 = {d:3d}: m >= {m:5d}")

    # 验证泛化界
    print("\n泛化误差界分析 (m=1000):")
    for d in [10, 50, 100]:
        bound = pac.generalization_bound_vc(m=1000, vc_dim=d)
        satisfies, upper = pac.probably_approximately_correct(
            empirical_risk=0.05, m=1000, vc_dim=d)
        print(f"  VC维 = {d:3d}: 泛化界 = {bound:.4f}, 满足PAC: {satisfies}")

    return pac

# 运行示例
if __name__ == "__main__":
    pac = pac_learning_example()

class ComplexityMeasures:
    """
    复杂度度量：VC维和Rademacher复杂度
    """

    @staticmethod  # @staticmethod静态方法，无需实例即可调用
    def vc_dim_linear(d: int) -> int:
        """d维线性分类器的VC维"""
        return d + 1

    @staticmethod
    def vc_dim_neural_network(weights: int, depth: int = None) -> int:
        """
        神经网络的VC维上界

        Args:
            weights: 权重参数数量
            depth: 网络深度（可选）

        Returns:
            VC维上界
        """
        # 简化上界: O(W log W)
        if depth is None:
            return int(weights * np.log2(weights))
        else:
            # 更紧的上界考虑深度
            return int(depth * weights * np.log2(weights))

    @staticmethod
    def rademacher_complexity_bound(m: int, vc_dim: int) -> float:
        """
        基于VC维的Rademacher复杂度上界

        R_m(H) <= sqrt(2d*ln(em/d)/m)
        """
        if m <= vc_dim:
            return 1.0
        return np.sqrt(2 * vc_dim * np.log(np.e * m / vc_dim) / m)

    @staticmethod
    def empirical_rademacher(X: np.ndarray,
                             hypothesis_class: Callable,
                             n_samples: int = 100) -> float:
        """
        估计经验Rademacher复杂度

        Args:
            X: 数据矩阵 (n_samples, n_features)
            hypothesis_class: 假设类函数
            n_samples: Rademacher采样次数

        Returns:
            经验Rademacher复杂度估计
        """
        n = len(X)
        rademacher_values = []

        for _ in range(n_samples):
            # 生成Rademacher随机变量
            sigma = np.random.choice([-1, 1], size=n)

            # 在假设类上最大化
            # 简化版本：随机采样假设
            max_correlation = 0
            for _ in range(50):  # 采样假设
                h = hypothesis_class(X)
                correlation = np.mean(sigma * h)
                max_correlation = max(max_correlation, correlation)

            rademacher_values.append(max_correlation)

        return np.mean(rademacher_values)

    @staticmethod
    def margin_bound(rademacher: float,
                     margin: float,
                     m: int,
                     delta: float = 0.05) -> float:
        """
        基于间隔的泛化界

        对于间隔为margin的分类器，泛化误差上界
        """
        # R(h) <= R_emp^margin(h) + (2/margin)*R_m(H) + sqrt(ln(1/delta)/(2m))
        bound = (2.0 / margin) * rademacher + np.sqrt(np.log(1.0 / delta) / (2 * m))
        return bound

def complexity_analysis_example():
    """复杂度分析示例"""
    print("\n" + "=" * 60)
    print("复杂度度量分析")
    print("=" * 60)

    cm = ComplexityMeasures()

    # VC维分析
    print("\n1. VC维分析")
    print("-" * 40)

    dimensions = [10, 100, 1000, 10000]
    for d in dimensions:
        vc_linear = cm.vc_dim_linear(d)
        print(f"  输入维度 d={d:5d}: 线性分类器VC维 = {vc_linear:6d}")

    # 神经网络VC维
    print("\n  神经网络VC维上界:")
    weight_counts = [1000, 10000, 100000, 1000000]
    for w in weight_counts:
        vc_nn = cm.vc_dim_neural_network(w)
        print(f"    权重数 W={w:8d}: VC维上界 ≈ {vc_nn:10d}")

    # Rademacher复杂度
    print("\n2. Rademacher复杂度分析")
    print("-" * 40)

    sample_sizes = [100, 1000, 10000, 100000]
    vc_dim = 100

    for m in sample_sizes:
        rad_bound = cm.rademacher_complexity_bound(m, vc_dim)
        print(f"  样本数 m={m:6d}, VC维 d={vc_dim:3d}: "
              f"Rademacher上界 = {rad_bound:.4f}")

    # 间隔界
    print("\n3. 基于间隔的泛化界")
    print("-" * 40)

    margins = [0.1, 0.5, 1.0, 2.0]
    m, rad = 1000, 0.1

    for margin in margins:
        bound = cm.margin_bound(rad, margin, m)
        print(f"  间隔 γ={margin:.1f}: 泛化界 = {bound:.4f}")

# 运行示例
if __name__ == "__main__":
    complexity_analysis_example()

class UniformConvergence:
    """
    均匀收敛理论与泛化界
    """

    @staticmethod
    def mcdiarmid_inequality(bounded_differences: list[float],
                            t: float) -> float:
        """
        McDiarmid不等式

        P(|f - E[f]| >= t) <= 2*exp(-2t^2 / sum(c_i^2))

        Args:
            bounded_differences: 有界差分 c_i 列表
            t: 偏离阈值

        Returns:
            概率上界
        """
        sum_c_squared = sum(c**2 for c in bounded_differences)
        bound = 2 * np.exp(-2 * t**2 / sum_c_squared)
        return bound

    @staticmethod
    def covering_number_bound(m: int,
                               covering_number: int,
                               epsilon: float) -> float:
        """
        基于覆盖数的Rademacher复杂度上界

        Args:
            m: 样本数
            covering_number: 覆盖数 N(H, epsilon)
            epsilon: 覆盖半径

        Returns:
            Rademacher复杂度上界
        """
        bound = epsilon + np.sqrt(2 * np.log(covering_number) / m)
        return bound

    @staticmethod
    def uniform_convergence_bound(rademacher: float,
                                   m: int,
                                   delta: float = 0.05) -> float:
        """
        均匀收敛界

        sup|R(h) - R_hat(h)| <= 2*R_m(H) + sqrt(ln(1/delta)/(2m))
        """
        bound = 2 * rademacher + np.sqrt(np.log(1.0 / delta) / (2 * m))
        return bound

    @staticmethod
    def fat_shattering_dimension(margin: float,
                                  data: np.ndarray) -> int:
        """
        胖粉碎维（Fat-shattering dimension）

        实值函数类的推广VC维

        Args:
            margin: 间隔参数
            data: 数据集

        Returns:
            胖粉碎维估计
        """
        # 简化实现：基于间隔的VC维估计
        n_samples, n_features = data.shape

        # 胖粉碎维与间隔成反比
        # 对于线性分类器: fat_gamma = O(1/gamma^2)
        fat_dim = int(n_features / (margin ** 2))

        return min(fat_dim, n_samples)

def uniform_convergence_example():
    """均匀收敛示例"""
    print("\n" + "=" * 60)
    print("均匀收敛理论")
    print("=" * 60)

    uc = UniformConvergence()

    # McDiarmid不等式
    print("\n1. McDiarmid不等式")
    print("-" * 40)

    # 经验风险的有界差分
    n_samples = 1000
    bounded_diffs = [1.0 / n_samples] * n_samples  # 每个样本影响最大1/n

    thresholds = [0.01, 0.05, 0.1, 0.2]
    for t in thresholds:
        prob_bound = uc.mcdiarmid_inequality(bounded_diffs, t)
        print(f"  偏离阈值 t={t:.2f}: P(|f-E[f]|>={t}) <= {prob_bound:.6f}")

    # 均匀收敛界
    print("\n2. 均匀收敛界")
    print("-" * 40)

    m_values = [100, 500, 1000, 5000]
    rademacher = 0.1

    for m in m_values:
        bound = uc.uniform_convergence_bound(rademacher, m)
        print(f"  样本数 m={m:4d}: 均匀收敛界 = {bound:.4f}")

    # 胖粉碎维
    print("\n3. 胖粉碎维（实值函数类）")
    print("-" * 40)

    margins = [0.1, 0.5, 1.0, 2.0]
    data = np.random.randn(1000, 100)

    for margin in margins:
        fat_dim = uc.fat_shattering_dimension(margin, data)
        print(f"  间隔 γ={margin:.1f}: 胖粉碎维 ≈ {fat_dim}")

# 运行示例
if __name__ == "__main__":
    uniform_convergence_example()

class DeepGeneralizationMystery:
    """
    深度学习泛化之谜的分析
    """

    @staticmethod
    def effective_complexity(parameters: int,
                            training_epochs: int,
                            batch_size: int,
                            learning_rate: float) -> float:
        """
        估计有效复杂度（考虑优化过程）

        启发式：有效复杂度 < 参数数量

        Args:
            parameters: 参数数量
            training_epochs: 训练轮数
            batch_size: 批量大小
            learning_rate: 学习率

        Returns:
            有效复杂度估计
        """
        # 启发式：有效复杂度与训练步数相关
        steps = training_epochs * (parameters / batch_size)

        # 学习率影响收敛到的解的复杂度
        lr_factor = np.tanh(learning_rate * 10)  # 归一化

        # 有效复杂度上界
        effective = min(parameters, steps * lr_factor)

        return effective

    @staticmethod
    def flat_minima_generalization(hessian_eigenvalues: np.ndarray,
                                   empirical_risk: float) -> float:
        """
        基于Hessian特征值的平坦极小值泛化分析

        平坦极小值通常泛化更好

        Args:
            hessian_eigenvalues: Hessian矩阵特征值
            empirical_risk: 经验风险

        Returns:
            泛化能力评分（越高越好）
        """
        # 计算Hessian的迹（特征值之和）
        trace = np.sum(hessian_eigenvalues)

        # 计算特征值分布的集中度
        eigenvalue_variance = np.var(hessian_eigenvalues)

        # 平坦度指标：迹越小、方差越小越平坦
        flatness_score = 1.0 / (1.0 + trace + eigenvalue_variance)

        # 综合泛化评分
        generalization_score = (1 - empirical_risk) * flatness_score

        return generalization_score

    @staticmethod
    def double_descent_curve(n_samples: int,
                             max_params: int,
                             noise_level: float = 0.1) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
        """
        生成双下降曲线

        Args:
            n_samples: 样本数
            max_params: 最大参数数
            noise_level: 噪声水平

        Returns:
            (参数数数组, 训练误差, 测试误差)
        """
        param_counts = np.arange(1, max_params + 1)
        train_errors = []
        test_errors = []

        for n_params in param_counts:
            if n_params < n_samples:
                # 欠参数化区域：传统U型曲线
                underfitting = np.exp(-n_params / (n_samples * 0.3))
                overfitting = 0.1 * (n_params / n_samples) ** 2
                test_err = underfitting + overfitting + noise_level
            elif n_params == n_samples:
                # 插值阈值：峰值
                test_err = 0.5 + noise_level
            else:
                # 过参数化区域：良性过拟合
                interpolation = noise_level * (n_samples / n_params)
                test_err = interpolation + 0.05

            # 训练误差单调递减
            train_err = max(0.01, noise_level * np.exp(-n_params / n_samples))

            train_errors.append(train_err)
            test_errors.append(test_err)

        return param_counts, np.array(train_errors), np.array(test_errors)  # np.array创建NumPy数组

    @staticmethod
    def pac_bayes_bound(empirical_risk: float,
                       kl_divergence: float,
                       m: int,
                       delta: float = 0.05) -> float:
        """
        PAC-Bayes泛化界

        Args:
            empirical_risk: 经验风险
            kl_divergence: KL散度 D(Q||P)
            m: 样本数
            delta: 置信度

        Returns:
            PAC-Bayes泛化界
        """
        # KL项
        kl_term = (kl_divergence + np.log(2 * np.sqrt(m) / delta)) / m

        # PAC-Bayes界
        bound = empirical_risk + np.sqrt(kl_term / 2)

        return bound

def deep_generalization_example():
    """深度学习泛化之谜示例"""
    print("\n" + "=" * 60)
    print("深度学习泛化之谜")
    print("=" * 60)

    dgm = DeepGeneralizationMystery()

    # 有效复杂度
    print("\n1. 有效复杂度分析")
    print("-" * 40)

    configs = [
        {"params": 10000, "epochs": 100, "batch": 32, "lr": 0.001},
        {"params": 1000000, "epochs": 50, "batch": 128, "lr": 0.01},
        {"params": 100000000, "epochs": 10, "batch": 256, "lr": 0.1},
    ]

    for config in configs:
        effective = dgm.effective_complexity(
            config["params"], config["epochs"],
            config["batch"], config["lr"]
        )
        ratio = effective / config["params"]
        print(f"  参数: {config['params']:10d}, 有效复杂度: {effective:12.1f}, "
              f"比例: {ratio:.4f}")

    # 双下降曲线
    print("\n2. 双下降现象")
    print("-" * 40)

    n_samples = 100
    max_params = 300

    params, train_err, test_err = dgm.double_descent_curve(
        n_samples, max_params, noise_level=0.1)

    print(f"  样本数: {n_samples}")
    print(f"  插值阈值（参数数=样本数）: {n_samples}")
    print(f"  阈值处测试误差: {test_err[n_samples-1]:.4f}")
    print(f"  最大参数处测试误差: {test_err[-1]:.4f}")

    # 平坦极小值
    print("\n3. 平坦极小值分析")
    print("-" * 40)

    # 模拟Hessian特征值
    sharp_eigenvalues = np.array([10, 8, 5, 3, 2, 1, 0.5, 0.1])  # 尖锐
    flat_eigenvalues = np.array([2, 1.5, 1, 0.8, 0.5, 0.3, 0.2, 0.1])  # 平坦

    sharp_score = dgm.flat_minima_generalization(sharp_eigenvalues, 0.05)
    flat_score = dgm.flat_minima_generalization(flat_eigenvalues, 0.05)

    print(f"  尖锐极小值泛化评分: {sharp_score:.4f}")
    print(f"  平坦极小值泛化评分: {flat_score:.4f}")
    print(f"  平坦/尖锐比值: {flat_score/sharp_score:.2f}x")

    # PAC-Bayes
    print("\n4. PAC-Bayes界")
    print("-" * 40)

    kl_divs = [1, 10, 100, 1000]
    for kl in kl_divs:
        bound = dgm.pac_bayes_bound(0.05, kl, m=1000)
        print(f"  KL散度={kl:4d}: PAC-Bayes界 = {bound:.4f}")

# 运行示例
if __name__ == "__main__":
    deep_generalization_example()

class ManifoldLearning:
    """
    流形学习理论与实现
    """

    @staticmethod
    def lle(X: np.ndarray, n_components: int = 2, n_neighbors: int = 10) -> np.ndarray:
        """
        局部线性嵌入 (LLE)

        Args:
            X: 高维数据 (n_samples, n_features)
            n_components: 目标维度
            n_neighbors: 近邻数

        Returns:
            低维嵌入 (n_samples, n_components)
        """
        n_samples = X.shape[0]

        # 1. 计算k近邻
        from scipy.spatial.distance import cdist
        distances = cdist(X, X)
        neighbors = np.argsort(distances, axis=1)[:, 1:n_neighbors+1]

        # 2. 计算重构权重
        W = np.zeros((n_samples, n_samples))

        for i in range(n_samples):
            # 中心化的近邻点
            Xi = X[neighbors[i]] - X[i]

            # 局部协方差矩阵
            C = Xi @ Xi.T

            # 添加正则化
            C += np.eye(n_neighbors) * 1e-3 * np.trace(C)

            # 求解权重
            w = np.linalg.solve(C, np.ones(n_neighbors))
            w /= np.sum(w)  # 归一化

            W[i, neighbors[i]] = w

        # 3. 计算低维嵌入
        M = (np.eye(n_samples) - W).T @ (np.eye(n_samples) - W)

        # 特征值分解
        eigenvalues, eigenvectors = np.linalg.eigh(M)

        # 选择最小的非零特征值对应的特征向量
        Y = eigenvectors[:, 1:n_components+1]

        return Y

    @staticmethod
    def intrinsic_dimension(X: np.ndarray,
                           k: int = 10,
                           method: str = 'correlation') -> float:
        """
        估计数据的内在维度

        Args:
            X: 数据矩阵
            k: 近邻数
            method: 估计方法

        Returns:
            内在维度估计
        """
        from scipy.spatial.distance import cdist

        n_samples = X.shape[0]
        distances = cdist(X, X)

        # 对每个点计算到k近邻的距离
        neighbor_distances = np.sort(distances, axis=1)[:, 1:k+1]

        if method == 'correlation':
            # 基于距离相关性的方法
            # 在d维流形上，体积按r^d增长

            # 计算平均距离
            r = np.mean(neighbor_distances, axis=0)
            k_values = np.arange(1, k+1)

            # 拟合 log(k) ~ d * log(r)
            log_r = np.log(r + 1e-10)
            log_k = np.log(k_values)

            # 线性回归估计斜率
            slope = np.polyfit(log_r, log_k, 1)[0]

            return slope

        elif method == 'mle':
            # 最大似然估计
            # d_hat = (1/(N-1)) * sum(log(r_k / r_i))^{-1}

            dims = []
            for i in range(n_samples):
                r_i = neighbor_distances[i, :-1]
                r_k = neighbor_distances[i, -1]

                if r_k > 0:
                    d_hat = 1.0 / (np.mean(np.log(r_k / (r_i + 1e-10))) + 1e-10)
                    dims.append(d_hat)

            return np.median(dims)

        return 0.0

    @staticmethod
    def manifold_hypothesis_validation(X: np.ndarray,
                                       embedding_methods: list[str] = None) -> dict:
        """
        验证流形假设

        Args:
            X: 数据矩阵
            embedding_methods: 嵌入方法列表

        Returns:
            验证结果字典
        """
        if embedding_methods is None:
            embedding_methods = ['pca', 'lle', 'tsne']

        results = {}

        # 估计内在维度
        intrinsic_dim = ManifoldLearning.intrinsic_dimension(X)
        results['intrinsic_dimension'] = intrinsic_dim
        results['ambient_dimension'] = X.shape[1]
        results['dimension_ratio'] = intrinsic_dim / X.shape[1]

        print(f"环境维度: {X.shape[1]}")
        print(f"估计内在维度: {intrinsic_dim:.2f}")
        print(f"维度比: {intrinsic_dim/X.shape[1]:.4f}")

        return results

def manifold_learning_example():
    """流形学习示例"""
    print("\n" + "=" * 60)
    print("流形学习理论")
    print("=" * 60)

    ml = ManifoldLearning()

    # 生成瑞士卷数据
    from sklearn.datasets import make_swiss_roll
    X, t = make_swiss_roll(n_samples=1000, noise=0.1, random_state=42)

    print("\n1. 瑞士卷数据流形分析")
    print("-" * 40)
    print(f"  数据形状: {X.shape}")

    # 估计内在维度
    intrinsic_dim = ml.intrinsic_dimension(X, k=10)
    print(f"  估计内在维度: {intrinsic_dim:.2f} (真实: 2)")

    # LLE降维
    print("\n2. LLE降维")
    print("-" * 40)

    X_lle = ml.lle(X, n_components=2, n_neighbors=10)
    print(f"  降维后形状: {X_lle.shape}")

    # 验证流形假设
    print("\n3. 流形假设验证")
    print("-" * 40)

    results = ml.manifold_hypothesis_validation(X)
    print(f"  维度压缩比: {results['dimension_ratio']:.4f}")

    return X, X_lle

# 运行示例
if __name__ == "__main__":
    X_swiss, X_lle = manifold_learning_example()

class DistributedRepresentation:
    """
    分布式表示与层次表示分析
    """

    @staticmethod
    def representation_capacity(n_neurons: int,
                                n_layers: int,
                                binary: bool = True) -> int:
        """
        计算网络的表示能力

        Args:
            n_neurons: 每层神经元数
            n_layers: 层数
            binary: 是否二元激活

        Returns:
            可表示的模式数
        """
        if binary:
            # 二元激活：每层有 2^n 种状态
            patterns_per_layer = 2 ** n_neurons
            # 总表示能力（考虑层次组合）
            total_capacity = patterns_per_layer ** n_layers
        else:
            # 连续激活：无限可能，考虑精度
            precision = 16  # 假设16位精度
            states_per_neuron = 2 ** precision
            total_capacity = states_per_neuron ** (n_neurons * n_layers)

        return total_capacity

    @staticmethod
    def compositionality_analysis(features: np.ndarray,
                                   compositions: list[tuple[int, ...]]) -> dict:
        """
        分析表示的组合性

        Args:
            features: 特征矩阵 (n_samples, n_features)
            compositions: 组合关系列表

        Returns:
            组合性分析结果
        """
        results = {}

        # 计算特征之间的相关性
        feature_corr = np.corrcoef(features.T)

        # 分析组合关系
        composition_scores = []
        for comp in compositions:
            if len(comp) >= 2:
                # 检查组合特征是否线性可分
                comp_features = features[:, comp]

                # 简单线性可分性检验
                # 如果组合特征的方差大于单个特征，说明有组合效应
                joint_variance = np.var(np.sum(comp_features, axis=1))
                individual_variances = [np.var(features[:, i]) for i in comp]

                composition_score = joint_variance / (sum(individual_variances) + 1e-10)
                composition_scores.append(composition_score)

        results['mean_composition_score'] = np.mean(composition_scores)
        results['feature_correlation_mean'] = np.mean(np.abs(feature_corr[np.triu_indices_from(feature_corr, k=1)]))

        return results

    @staticmethod
    def hierarchical_decomposition(X: np.ndarray,
                                   layer_sizes: list[int]) -> list[np.ndarray]:
        """
        层次分解数据表示

        Args:
            X: 输入数据
            layer_sizes: 每层的大小

        Returns:
            各层表示列表
        """
        representations = [X]
        current = X

        for size in layer_sizes:
            # 使用PCA作为简单的层次分解
            from sklearn.decomposition import PCA

            pca = PCA(n_components=min(size, current.shape[1]))
            current = pca.fit_transform(current)
            representations.append(current)

        return representations

    @staticmethod
    def disentanglement_score(representations: np.ndarray,
                              ground_truth_factors: np.ndarray) -> float:
        """
        计算解耦分数

        衡量表示中每个维度是否对应独立的生成因子

        Args:
            representations: 学习到的表示 (n_samples, n_features)
            ground_truth_factors: 真实因子 (n_samples, n_factors)

        Returns:
            解耦分数 [0, 1]
        """
        from sklearn.ensemble import RandomForestRegressor
        from sklearn.metrics import r2_score

        n_factors = ground_truth_factors.shape[1]
        n_features = representations.shape[1]

        # 对每个真实因子，找到最相关的表示维度
        importance_matrix = np.zeros((n_factors, n_features))

        for i in range(n_factors):
            # 使用随机森林评估每个特征的重要性
            rf = RandomForestRegressor(n_estimators=10, random_state=42)
            rf.fit(representations, ground_truth_factors[:, i])
            importance_matrix[i] = rf.feature_importances_

        # 计算解耦程度：每个因子应该主要由一个特征解释
        max_importance_per_factor = np.max(importance_matrix, axis=1)
        disentanglement = np.mean(max_importance_per_factor)

        return disentanglement

def distributed_representation_example():
    """分布式表示示例"""
    print("\n" + "=" * 60)
    print("分布式表示与层次表示")
    print("=" * 60)

    dr = DistributedRepresentation()

    # 表示能力分析
    print("\n1. 表示能力分析")
    print("-" * 40)

    configs = [
        {"neurons": 100, "layers": 1},
        {"neurons": 100, "layers": 3},
        {"neurons": 1000, "layers": 3},
        {"neurons": 1000, "layers": 10},
    ]

    for config in configs:
        capacity = dr.representation_capacity(config["neurons"], config["layers"])
        log_capacity = np.log10(capacity) if capacity > 0 else 0
        print(f"  神经元={config['neurons']:4d}, 层数={config['layers']:2d}: "
              f"log10(容量) = {log_capacity:.2f}")

    # 层次分解
    print("\n2. 层次分解")
    print("-" * 40)

    # 生成层次数据
    np.random.seed(42)
    n_samples = 1000

    # 底层特征（如边缘）
    low_level = np.random.randn(n_samples, 100)

    # 中层特征（如形状）- 底层特征的组合
    mid_level = np.hstack([
        np.mean(low_level[:, :50], axis=1, keepdims=True),
        np.std(low_level[:, 50:], axis=1, keepdims=True),
    ])

    # 高层特征（如物体）- 中层特征的组合
    high_level = np.hstack([
        (mid_level[:, 0] > 0).astype(float).reshape(-1, 1),  # 链式调用，连续执行多个方法  # reshape重塑张量形状
        (mid_level[:, 1] < 0.5).astype(float).reshape(-1, 1),
    ])

    X = np.hstack([low_level, mid_level, high_level])

    representations = dr.hierarchical_decomposition(X, [50, 20, 5])

    print(f"  原始数据维度: {X.shape}")
    for i, rep in enumerate(representations[1:], 1):  # enumerate同时获取索引和元素  # 切片操作取子序列
        print(f"  第{i}层表示维度: {rep.shape}")

    # 组合性分析
    print("\n3. 组合性分析")
    print("-" * 40)

    # 模拟组合关系
    compositions = [(0, 1), (2, 3), (4, 5, 6)]
    comp_results = dr.compositionality_analysis(X[:, :10], compositions)

    print(f"  平均组合分数: {comp_results['mean_composition_score']:.4f}")
    print(f"  特征相关性: {comp_results['feature_correlation_mean']:.4f}")

# 运行示例
if __name__ == "__main__":
    distributed_representation_example()

class RepresentationLearningTheory:
    """
    表示学习的理论分析
    """

    @staticmethod
    def mutual_information_lower_bound(X: np.ndarray,
                                        Z: np.ndarray,
                                        method: str = 'knn') -> float:
        """
        估计互信息的下界

        Args:
            X: 原始数据
            Z: 学习到的表示
            method: 估计方法

        Returns:
            互信息估计
        """
        if method == 'knn':
            # 基于k近邻的互信息估计（Kraskov-Stögbauer-Grassberger）
            from sklearn.neighbors import NearestNeighbors

            n_samples = len(X)
            k = min(5, n_samples // 10)

            # 联合空间 (X, Z)
            XZ = np.hstack([X, Z])

            # 计算近邻
            nn_x = NearestNeighbors(n_neighbors=k+1).fit(X)
            nn_z = NearestNeighbors(n_neighbors=k+1).fit(Z)
            nn_xz = NearestNeighbors(n_neighbors=k+1).fit(XZ)

            distances_x, _ = nn_x.kneighbors(X)
            distances_z, _ = nn_z.kneighbors(Z)
            distances_xz, _ = nn_xz.kneighbors(XZ)

            # 互信息估计
            eps_x = 2 * distances_x[:, k]
            eps_z = 2 * distances_z[:, k]
            eps_xz = 2 * distances_xz[:, k]

            # 简化估计
            mi_estimate = (np.mean(np.log(eps_x + 1e-10)) +
                          np.mean(np.log(eps_z + 1e-10)) -
                          np.mean(np.log(eps_xz + 1e-10)))

            return max(0, mi_estimate)

        return 0.0

    @staticmethod
    def information_bottleneck_objective(X: np.ndarray,
                                          Y: np.ndarray,
                                          Z: np.ndarray,
                                          beta: float = 1.0) -> float:
        """
        信息瓶颈目标函数

        min I(X; Z) - beta * I(Z; Y)

        Args:
            X: 输入
            Y: 目标
            Z: 表示
            beta: 权衡参数

        Returns:
            IB目标值
        """
        # 估计各项互信息
        I_XZ = RepresentationLearningTheory.mutual_information_lower_bound(X, Z)
        I_ZY = RepresentationLearningTheory.mutual_information_lower_bound(Z, Y.reshape(-1, 1) if Y.ndim == 1 else Y)

        # IB目标
        objective = I_XZ - beta * I_ZY

        return objective

    @staticmethod
    def contrastive_learning_bound(n_negative: int) -> float:
        """
        对比学习的互信息下界

        InfoNCE >= -log(N) - I(Z; Y)

        Args:
            n_negative: 负样本数

        Returns:
            下界常数
        """
        return -np.log(n_negative)

    @staticmethod
    def representation_quality_metrics(Z: np.ndarray,
                                        Y: np.ndarray,
                                        task_type: str = 'classification') -> dict:
        """
        评估表示质量的多个指标

        Args:
            Z: 表示
            Y: 标签
            task_type: 任务类型

        Returns:
            质量指标字典
        """
        metrics = {}

        # 1. 线性可分性
        from sklearn.linear_model import LogisticRegression
        from sklearn.model_selection import cross_val_score

        if task_type == 'classification':
            clf = LogisticRegression(max_iter=1000, random_state=42)
            scores = cross_val_score(clf, Z, Y, cv=5)
            metrics['linear_separability'] = np.mean(scores)

        # 2. 表示的紧凑性
        # 计算特征值的熵
        from sklearn.decomposition import PCA
        pca = PCA().fit(Z)
        eigenvalues = pca.explained_variance_ratio_

        # 有效维度
        entropy = -np.sum(eigenvalues * np.log(eigenvalues + 1e-10))
        metrics['effective_dimensionality'] = np.exp(entropy)
        metrics['explained_variance_top10'] = np.sum(eigenvalues[:10])

        # 3. 类内/类间方差比
        if task_type == 'classification':
            classes = np.unique(Y)
            overall_mean = np.mean(Z, axis=0)

            # 类间散度
            S_b = np.zeros((Z.shape[1], Z.shape[1]))
            # 类内散度
            S_w = np.zeros((Z.shape[1], Z.shape[1]))

            for c in classes:
                Z_c = Z[Y == c]
                mean_c = np.mean(Z_c, axis=0)

                S_b += len(Z_c) * np.outer(mean_c - overall_mean, mean_c - overall_mean)
                S_w += np.cov(Z_c.T) * (len(Z_c) - 1)

            # Fisher ratio
            if np.linalg.det(S_w) > 1e-10:
                fisher_ratio = np.trace(S_b) / np.trace(S_w)
                metrics['fisher_ratio'] = fisher_ratio

        return metrics

def representation_theory_example():
    """表示学习理论示例"""
    print("\n" + "=" * 60)
    print("表示学习理论")
    print("=" * 60)

    rlt = RepresentationLearningTheory()

    # 生成模拟数据
    np.random.seed(42)
    n_samples = 500

    # 输入数据
    X = np.random.randn(n_samples, 50)

    # 目标（基于X的某些维度）
    Y = (X[:, 0] + X[:, 1] > 0).astype(int)

    # 学习到的表示（模拟）
    Z = np.hstack([
        X[:, :2] + 0.1 * np.random.randn(n_samples, 2),  # 相关信息
        0.1 * np.random.randn(n_samples, 8)  # 噪声
    ])

    # 互信息估计
    print("\n1. 互信息分析")
    print("-" * 40)

    I_XZ = rlt.mutual_information_lower_bound(X, Z)
    I_ZY = rlt.mutual_information_lower_bound(Z, Y.reshape(-1, 1))

    print(f"  I(X; Z) ≈ {I_XZ:.4f}")
    print(f"  I(Z; Y) ≈ {I_ZY:.4f}")

    # 信息瓶颈
    print("\n2. 信息瓶颈分析")
    print("-" * 40)

    for beta in [0.1, 1.0, 10.0]:
        ib_obj = rlt.information_bottleneck_objective(X, Y, Z, beta)
        print(f"  β={beta:5.1f}: IB目标 = {ib_obj:.4f}")

    # 对比学习界
    print("\n3. 对比学习理论界")
    print("-" * 40)

    for n_neg in [10, 100, 1000, 10000]:
        bound = rlt.contrastive_learning_bound(n_neg)
        print(f"  负样本数={n_neg:5d}: 下界 = {bound:.4f}")

    # 表示质量评估
    print("\n4. 表示质量指标")
    print("-" * 40)

    metrics = rlt.representation_quality_metrics(Z, Y)

    for metric, value in metrics.items():
        print(f"  {metric}: {value:.4f}")

# 运行示例
if __name__ == "__main__":
    representation_theory_example()

class CriticalPointAnalysis:
    """
    临界点分析与优化景观
    """

    @staticmethod
    def classify_critical_point(hessian_eigenvalues: np.ndarray,
                                 tolerance: float = 1e-6) -> str:
        """
        根据Hessian特征值分类临界点

        Args:
            hessian_eigenvalues: Hessian矩阵特征值
            tolerance: 数值容差

        Returns:
            临界点类型
        """
        n_positive = np.sum(hessian_eigenvalues > tolerance)
        n_negative = np.sum(hessian_eigenvalues < -tolerance)
        n_zero = len(hessian_eigenvalues) - n_positive - n_negative

        if n_zero > len(hessian_eigenvalues) * 0.1:
            return "degenerate"  # 退化临界点
        elif n_positive > 0 and n_negative > 0:
            return "saddle"  # 鞍点
        elif n_positive == len(hessian_eigenvalues):
            return "local_minimum"  # 局部极小值
        elif n_negative == len(hessian_eigenvalues):
            return "local_maximum"  # 局部极大值
        else:
            return "unknown"

    @staticmethod
    def hessian_spectrum_analysis(model: torch.nn.Module,
                                   data_loader: torch.utils.data.DataLoader,
                                   loss_fn: Callable,
                                   n_top_eigenvalues: int = 20) -> dict:
        """
        分析Hessian矩阵的谱

        Args:
            model: PyTorch模型
            data_loader: 数据加载器
            loss_fn: 损失函数
            n_top_eigenvalues: 计算的最大特征值数量

        Returns:
            谱分析结果
        """
        model.eval()  # eval()开启评估模式（关闭Dropout等）

        # 收集梯度
        for batch in data_loader:
            inputs, targets = batch

            model.zero_grad()  # 清零梯度，防止梯度累积
            outputs = model(inputs)
            loss = loss_fn(outputs, targets)
            loss.backward()  # 反向传播计算梯度

            # 提取梯度
            grads = []
            for param in model.parameters():
                if param.grad is not None:
                    grads.append(param.grad.flatten())

            gradient = torch.cat(grads)
            break  # 只用一个batch

        # 使用幂迭代估计最大特征值
        def hessian_vector_product(v):
            """计算Hessian-向量积"""
            model.zero_grad()
            for batch in data_loader:
                inputs, targets = batch
                outputs = model(inputs)
                loss = loss_fn(outputs, targets)

                # 一阶梯度
                grads = torch.autograd.grad(loss, model.parameters(), create_graph=True)
                flat_grads = torch.cat([g.flatten() for g in grads])

                # 梯度与v的内积
                grad_v = torch.sum(flat_grads * v)

                # 二阶梯度
                hvp = torch.autograd.grad(grad_v, model.parameters(), retain_graph=True)
                flat_hvp = torch.cat([g.flatten() for g in hvp])

                return flat_hvp

        # 幂迭代求最大特征值
        v = torch.randn_like(gradient)
        v = v / torch.norm(v)

        for _ in range(50):  # 迭代次数
            Hv = hessian_vector_product(v)
            v = Hv / torch.norm(Hv)

        # 瑞利商估计特征值
        Hv = hessian_vector_product(v)
        max_eigenvalue = torch.dot(v, Hv).item()  # .item()将单元素张量转为Python数值

        results = {
            'max_eigenvalue': max_eigenvalue,
            'gradient_norm': torch.norm(gradient).item(),
        }

        return results

    @staticmethod
    def loss_landscape_1d(model: torch.nn.Module,
                          direction: torch.Tensor,
                          data_loader: torch.utils.data.DataLoader,
                          loss_fn: Callable,
                          alpha_range: tuple[float, float] = (-1, 1),
                          n_points: int = 50) -> tuple[np.ndarray, np.ndarray]:
        """
        沿一维方向分析损失景观

        Args:
            model: 模型
            direction: 方向向量
            data_loader: 数据加载器
            loss_fn: 损失函数
            alpha_range: 参数范围
            n_points: 采样点数

        Returns:
            (alpha值, 损失值)
        """
        # 保存原始参数
        original_params = {}
        for name, param in model.named_parameters():
            original_params[name] = param.data.clone()

        alphas = np.linspace(alpha_range[0], alpha_range[1], n_points)
        losses = []

        for alpha in alphas:
            # 沿方向移动参数
            idx = 0
            for name, param in model.named_parameters():
                param_size = param.numel()
                param_direction = direction[idx:idx+param_size].reshape(param.shape)
                param.data = original_params[name] + alpha * param_direction
                idx += param_size

            # 计算损失
            total_loss = 0
            n_batches = 0
            with torch.no_grad():  # 禁用梯度计算，节省内存
                for batch in data_loader:
                    inputs, targets = batch
                    outputs = model(inputs)
                    loss = loss_fn(outputs, targets)
                    total_loss += loss.item()
                    n_batches += 1

            losses.append(total_loss / n_batches)

        # 恢复原始参数
        for name, param in model.named_parameters():
            param.data = original_params[name]

        return alphas, np.array(losses)

def critical_point_example():
    """临界点分析示例"""
    print("\n" + "=" * 60)
    print("临界点分析与优化景观")
    print("=" * 60)

    cpa = CriticalPointAnalysis()

    # 临界点分类
    print("\n1. 临界点分类")
    print("-" * 40)

    # 局部极小值
    min_eigenvalues = np.array([1.0, 0.5, 0.3, 0.1])
    print(f"  特征值 {min_eigenvalues}: {cpa.classify_critical_point(min_eigenvalues)}")

    # 局部极大值
    max_eigenvalues = np.array([-1.0, -0.5, -0.3, -0.1])
    print(f"  特征值 {max_eigenvalues}: {cpa.classify_critical_point(max_eigenvalues)}")

    # 鞍点
    saddle_eigenvalues = np.array([1.0, 0.5, -0.3, -0.1])
    print(f"  特征值 {saddle_eigenvalues}: {cpa.classify_critical_point(saddle_eigenvalues)}")

    # 退化临界点
    degenerate_eigenvalues = np.array([1.0, 0.0, 0.0, -0.1])
    print(f"  特征值 {degenerate_eigenvalues}: {cpa.classify_critical_point(degenerate_eigenvalues)}")

    # 高维分析
    print("\n2. 高维Hessian分析")
    print("-" * 40)

    dimensions = [10, 100, 1000, 10000]

    for dim in dimensions:
        # 模拟Hessian特征值分布
        # 大多数特征值接近0（鞍点特征）
        eigenvalues = np.random.randn(dim) * 0.1
        eigenvalues[0] = 1.0  # 一个正的大特征值
        eigenvalues[1] = -0.5  # 一个负的特征值

        critical_type = cpa.classify_critical_point(eigenvalues)
        n_positive = np.sum(eigenvalues > 1e-6)
        n_negative = np.sum(eigenvalues < -1e-6)

        print(f"  维度={dim:5d}: 类型={critical_type:15s}, "
              f"正特征值={n_positive:4d}, 负特征值={n_negative:4d}")

# 运行示例
if __name__ == "__main__":
    critical_point_example()