00 - 数学基础（进阶篇）¶

import numpy as np

def forward_backward_matrix_form(X, Y, weights, biases, activation='relu'):
    """
    使用矩阵形式进行前向和反向传播

    X: 输入 (n_features, m_samples)
    Y: 标签 (n_output, m_samples)
    weights: 权重列表
    biases: 偏置列表
    """
    m = X.shape[1]
    L = len(weights)

    # 存储中间结果
    Z = []
    A = [X]

    # 前向传播
    for l in range(L):
        Z_l = weights[l] @ A[l] + biases[l]
        Z.append(Z_l)

        if l < L - 1:  # 隐藏层
            A_l = np.maximum(0, Z_l) if activation == 'relu' else np.tanh(Z_l)
        else:  # 输出层
            # Softmax
            exp_Z = np.exp(Z_l - np.max(Z_l, axis=0, keepdims=True))
            A_l = exp_Z / np.sum(exp_Z, axis=0, keepdims=True)
        A.append(A_l)

    # 计算损失 (交叉熵)
    loss = -np.sum(Y * np.log(A[-1] + 1e-8)) / m

    # 反向传播
    dW = []
    db = []

    # 输出层梯度
    dZ = A[-1] - Y  # Softmax + CrossEntropy 的简化

    for l in range(L-1, -1, -1):
        # 权重梯度
        dW_l = (dZ @ A[l].T) / m
        dW.insert(0, dW_l)

        # 偏置梯度
        db_l = np.sum(dZ, axis=1, keepdims=True) / m
        db.insert(0, db_l)

        if l > 0:  # 非输入层
            # 传播到前一层
            dA = weights[l].T @ dZ
            # 激活函数导数
            if activation == 'relu':
                dZ = dA * (Z[l-1] > 0)
            else:
                dZ = dA * (1 - np.tanh(Z[l-1])**2)

    return loss, dW, db

# 示例
np.random.seed(42)
n_features, n_hidden, n_output, m = 784, 256, 10, 100

X = np.random.randn(n_features, m)
Y = np.zeros((n_output, m))
for i in range(m):
    Y[np.random.randint(n_output), i] = 1

weights = [
    np.random.randn(n_hidden, n_features) * 0.01,
    np.random.randn(n_output, n_hidden) * 0.01
]
biases = [
    np.zeros((n_hidden, 1)),
    np.zeros((n_output, 1))
]

loss, dW, db = forward_backward_matrix_form(X, Y, weights, biases)
print(f"损失: {loss:.4f}")
print(f"dW[0] shape: {dW[0].shape}")
print(f"dW[1] shape: {dW[1].shape}")

def conv2d_backward(dZ, X, W, stride=1, padding=0):
    """
    卷积层的反向传播

    dZ: 输出梯度 (batch, out_channels, out_h, out_w)
    X: 输入 (batch, in_channels, h, w)
    W: 权重 (out_channels, in_channels, k, k)
    """
    batch, out_channels, out_h, out_w = dZ.shape
    _, in_channels, h, w = X.shape
    k = W.shape[2]

    # 权重梯度
    dW = np.zeros_like(W)
    for i in range(out_channels):
        for j in range(in_channels):
            for b in range(batch):
                dW[i, j] += conv2d_single_channel(
                    X[b, j:j+1],
                    dZ[b, i:i+1],
                    padding=0
                )[0, 0]
    dW /= batch

    # 输入梯度
    dX = np.zeros_like(X)
    W_rot = np.rot90(W, 2, axes=(2, 3))  # 旋转180度

    for b in range(batch):
        for i in range(in_channels):
            for j in range(out_channels):
                dX[b, i] += conv2d_single_channel(
                    dZ[b, j:j+1],
                    W_rot[j, i:i+1],
                    padding=k-1
                )[0, 0]

    return dW, dX

def conv2d_single_channel(X, W, padding=0):
    """单通道卷积"""
    if padding > 0:
        X = np.pad(X, ((0, 0), (0, 0), (padding, padding), (padding, padding)))

    batch, _, h, w = X.shape
    out_channels, _, k, _ = W.shape
    out_h, out_w = h - k + 1, w - k + 1

    output = np.zeros((batch, out_channels, out_h, out_w))
    for i in range(out_h):
        for j in range(out_w):
            patch = X[:, :, i:i+k, j:j+k]
            for oc in range(out_channels):
                output[:, oc, i, j] = np.sum(patch * W[oc:oc+1], axis=(1, 2, 3))

    return output

import torch

def hessian_vector_product(loss, params, v):
    """
    计算Hessian-vector积 Hv

    利用: Hv = ∇(∇L · v)
    """
    # 计算梯度
    grad = torch.autograd.grad(loss, params, create_graph=True, retain_graph=True)

    # 计算梯度和v的内积
    grad_v = sum((g * v_i).sum() for g, v_i in zip(grad, v))

    # 再求梯度
    Hv = torch.autograd.grad(grad_v, params, retain_graph=True)

    return Hv

# 示例
def test_hessian():
    x = torch.tensor([1.0, 2.0], requires_grad=True)

    # 函数: f(x) = x₁² + 2x₂² + 3x₁x₂
    f = x[0]**2 + 2*x[1]**2 + 3*x[0]*x[1]

    # Hessian应该是:
    # [[2, 3],
    #  [3, 4]]

    v = [torch.tensor([1.0, 0.0])]  # 测试向量
    Hv = hessian_vector_product(f, [x], v)

    print(f"Hv = {Hv[0]}")  # 应该是 [2, 3]

test_hessian()

class AutoDiff:
    """简单的自动微分实现"""

    def __init__(self, value, grad_fn=None):  # __init__构造方法，创建对象时自动调用
        self.value = value
        self.grad = 0.0
        self.grad_fn = grad_fn  # 梯度函数
        self.prev = []  # 前驱节点

    def __add__(self, other):
        result = AutoDiff(self.value + other.value)
        result.prev = [self, other]
        result.grad_fn = lambda g: (g, g)  # d(a+b)/da = 1, d(a+b)/db = 1  # lambda匿名函数，简洁定义单行函数
        return result

    def __mul__(self, other):
        result = AutoDiff(self.value * other.value)
        result.prev = [self, other]
        result.grad_fn = lambda g: (g * other.value, g * self.value)
        return result

    def backward(self, grad=1.0):
        """反向传播"""
        self.grad += grad

        if self.grad_fn is not None:
            grads = self.grad_fn(grad)
            for prev, g in zip(self.prev, grads):
                prev.backward(g)  # 反向传播计算梯度

# 示例: z = (x + y) * x
x = AutoDiff(2.0)
y = AutoDiff(3.0)
z = (x + y) * x

z.backward()

print(f"z = {z.value}")  # (2+3)*2 = 10
print(f"dz/dx = {x.grad}")  # 2x + y = 7
print(f"dz/dy = {y.grad}")  # x = 2

from scipy.optimize import minimize

def constrained_optimization_example():
    """
    示例: 在球面上最小化二次函数
    min x^T Q x
    s.t. ||x||^2 = 1
    """
    Q = np.diag([1, 2, 3])  # 目标函数

    def objective(x):
        return x.T @ Q @ x

    def gradient(x):
        return 2 * Q @ x

    # 约束: ||x||^2 = 1
    def constraint(x):
        return np.sum(x**2) - 1

    # 使用SLSQP方法
    result = minimize(
        objective,
        x0=np.array([1, 0, 0]),  # np.array创建NumPy数组
        method='SLSQP',
        jac=gradient,
        constraints={'type': 'eq', 'fun': constraint}
    )

    print(f"最优解: {result.x}")
    print(f"最优值: {result.fun}")
    print(f"约束满足: {constraint(result.x):.10f}")

    # 解析解应该是Q的最小特征值对应的特征向量
    eigenvalues, eigenvectors = np.linalg.eig(Q)
    min_idx = np.argmin(eigenvalues)
    print(f"\n解析解: {eigenvectors[:, min_idx]}")
    print(f"解析最优值: {eigenvalues[min_idx]}")

constrained_optimization_example()

def project_to_sphere(x):
    """投影到单位球面"""
    return x / np.linalg.norm(x)

def riemannian_gradient_descent(f, grad_f, x0, n_iter=100, lr=0.01):
    """
    球面上的黎曼梯度下降

    约束流形: M = {x : ||x|| = 1}
    """
    x = x0.copy()
    x = project_to_sphere(x)  # 确保在流形上

    history = []
    for i in range(n_iter):
        # 计算欧氏梯度
        euclidean_grad = grad_f(x)

        # 投影到切空间: v = (I - xx^T)g
        tangent_grad = euclidean_grad - np.dot(x, euclidean_grad) * x

        # 梯度下降
        x_new = x - lr * tangent_grad

        # 投影回流形（重正交化）
        x = project_to_sphere(x_new)

        history.append(f(x))

    return x, history

# 示例: 在球面上最小化二次型
Q = np.diag([1, 2, 3, 4, 5])
x0 = np.random.randn(5)

x_opt, history = riemannian_gradient_descent(
    lambda x: x.T @ Q @ x,
    lambda x: 2 * Q @ x,
    x0,
    n_iter=200,
    lr=0.1
)

print(f"最优解: {x_opt}")
print(f"最优值: {x_opt.T @ Q @ x_opt:.6f}")
print(f"范数: {np.linalg.norm(x_opt):.6f}")

class BayesianNetwork:
    """简单的贝叶斯网络实现"""

    def __init__(self):
        self.graph = {}  # 邻接表
        self.cpts = {}   # 条件概率表

    def add_node(self, node, parents=None, cpt=None):
        """添加节点"""
        if parents is None:
            parents = []
        self.graph[node] = parents
        self.cpts[node] = cpt

    def get_joint_probability(self, assignment):
        """计算联合概率"""
        prob = 1.0
        for node, value in assignment.items():
            parents = self.graph[node]
            parent_values = tuple(assignment[p] for p in parents)  # 列表推导式，简洁创建列表

            if parents:
                prob *= self.cpts[node][parent_values][value]
            else:
                prob *= self.cpts[node][value]

        return prob

    def get_marginal_probability(self, node, evidence=None):
        """计算条件概率 P(node=1 | evidence)（枚举推理）"""
        if evidence is None:
            evidence = {}

        # 需要边缘化掉的变量（排除查询节点和证据节点）
        all_nodes = list(self.graph.keys())
        marginal_nodes = [n for n in all_nodes if n != node and n not in evidence]

        def sum_over(nodes, fixed_assignment):
            """枚举nodes的所有赋值，累加联合概率"""
            if not nodes:
                assignment = {**evidence, **fixed_assignment}
                return self.get_joint_probability(assignment)

            total = 0.0
            n = nodes[0]
            for value in [0, 1]:  # 假设二值变量
                fixed_assignment[n] = value
                total += sum_over(nodes[1:], fixed_assignment)  # 切片操作取子序列
                del fixed_assignment[n]
            return total

        # 对查询节点的每个取值求边缘概率
        prob_node_1 = sum_over(marginal_nodes, {node: 1})
        prob_node_0 = sum_over(marginal_nodes, {node: 0})

        # 归一化得到条件概率
        total_evidence = prob_node_0 + prob_node_1
        if total_evidence == 0:
            return 0.0
        return prob_node_1 / total_evidence

# 示例: 简单的疾病诊断网络
# 结构: Smoking -> LungCancer <- Pollution
#              -> Bronchitis

bn = BayesianNetwork()

# 先验概率
bn.add_node('Smoking', cpt={0: 0.7, 1: 0.3})
bn.add_node('Pollution', cpt={0: 0.9, 1: 0.1})

# 条件概率表
bn.add_node('LungCancer',
            parents=['Smoking', 'Pollution'],
            cpt={
                (0, 0): {0: 0.99, 1: 0.01},
                (0, 1): {0: 0.95, 1: 0.05},
                (1, 0): {0: 0.90, 1: 0.10},
                (1, 1): {0: 0.70, 1: 0.30}
            })

bn.add_node('Bronchitis',
            parents=['Smoking'],
            cpt={
                0: {0: 0.95, 1: 0.05},
                1: {0: 0.60, 1: 0.40}
            })

# 查询: 给定吸烟，患肺癌的概率
prob = bn.get_marginal_probability('LungCancer', evidence={'Smoking': 1})
print(f"P(LungCancer | Smoking) = {prob:.4f}")

class MarkovRandomField:
    """马尔可夫随机场"""

    def __init__(self, variables, cliques):
        """
        variables: 变量列表
        cliques: 团列表，每个团是变量的子集
        """
        self.variables = variables
        self.cliques = cliques
        self.potentials = {}

    def set_potential(self, clique_idx, potential_func):
        """设置势函数"""
        self.potentials[clique_idx] = potential_func

    def compute_partition_function(self):
        """计算配分函数（枚举所有赋值）"""
        n = len(self.variables)
        Z = 0.0

        # 枚举所有2^n种赋值
        for assignment_bits in range(2**n):
            assignment = {}
            for i, var in enumerate(self.variables):  # enumerate同时获取索引和元素
                assignment[var] = (assignment_bits >> i) & 1

            # 计算非归一化概率
            energy = 0.0
            for idx, clique in enumerate(self.cliques):
                clique_assignment = tuple(assignment[v] for v in clique)
                if idx in self.potentials:
                    energy += self.potentials[idx](clique_assignment)

            Z += np.exp(energy)

        return Z

    def gibbs_sampling(self, n_samples=1000, burn_in=100):
        """Gibbs采样"""
        # 随机初始化
        state = {var: np.random.randint(2) for var in self.variables}
        samples = []

        for i in range(n_samples + burn_in):
            for var in self.variables:
                # 计算条件概率 P(var | 其他变量)
                probs = []
                for val in [0, 1]:
                    state[var] = val

                    # 计算能量
                    energy = 0.0
                    for idx, clique in enumerate(self.cliques):
                        if var in clique:
                            clique_assignment = tuple(state[v] for v in clique)
                            if idx in self.potentials:
                                energy += self.potentials[idx](clique_assignment)

                    probs.append(np.exp(energy))

                # 归一化并采样
                probs = np.array(probs) / sum(probs)
                state[var] = np.random.choice([0, 1], p=probs)

            if i >= burn_in:
                samples.append(state.copy())

        return samples

# 示例: Ising模型（格子上的MRF）
def create_ising_model(n, J=1.0, h=0.0):
    """
    创建Ising模型
    能量: E = -J * Σ_{i~j} x_i x_j - h * Σ_i x_i
    """
    variables = [(i, j) for i in range(n) for j in range(n)]

    # 团：相邻的变量对
    cliques = []
    for i in range(n):
        for j in range(n):
            if i < n - 1:
                cliques.append([(i, j), (i+1, j)])
            if j < n - 1:
                cliques.append([(i, j), (i, j+1)])

    mrf = MarkovRandomField(variables, cliques)

    # 设置势函数
    for idx in range(len(cliques)):
        # 工厂函数make_potential解决Python闭包延迟绑定问题：若直接用lambda，循环结束后所有势函数的j都会指向最后一个值
        def make_potential(j):
            return lambda assignment: J * (2*assignment[0] - 1) * (2*assignment[1] - 1)
        mrf.set_potential(idx, make_potential(idx))

    return mrf

# 创建3x3 Ising模型
ising = create_ising_model(3, J=0.5)
samples = ising.gibbs_sampling(n_samples=1000)

# 计算平均磁化
magnetization = np.mean([np.mean(list(s.values())) for s in samples])
print(f"平均磁化: {magnetization:.4f}")

class FactorGraph:
    """因子图与和积算法"""

    def __init__(self):
        self.variables = {}  # 变量 -> 邻居因子
        self.factors = {}    # 因子 -> (变量列表, 因子函数)
        self.messages = {}   # 消息缓存

    def add_variable(self, var, domain):
        """添加变量"""
        self.variables[var] = {'domain': domain, 'neighbors': []}

    def add_factor(self, factor_id, variables, factor_func):
        """添加因子"""
        self.factors[factor_id] = (variables, factor_func)
        for var in variables:
            self.variables[var]['neighbors'].append(factor_id)

    def sum_product(self, max_iter=10):
        """和积算法（信念传播）"""
        # 初始化消息
        for var in self.variables:
            for factor in self.variables[var]['neighbors']:
                self.messages[('v2f', var, factor)] = np.ones(
                    len(self.variables[var]['domain'])
                )
                self.messages[('f2v', factor, var)] = np.ones(
                    len(self.variables[var]['domain'])
                )

        # 迭代传递消息
        for _ in range(max_iter):
            # 变量到因子
            for var in self.variables:
                for factor in self.variables[var]['neighbors']:
                    msg = np.ones(len(self.variables[var]['domain']))
                    for other_factor in self.variables[var]['neighbors']:
                        if other_factor != factor:
                            msg *= self.messages[('f2v', other_factor, var)]
                    self.messages[('v2f', var, factor)] = msg / msg.sum()

            # 因子到变量
            for factor_id, (vars_in_factor, factor_func) in self.factors.items():
                for target_var in vars_in_factor:
                    # 收集其他变量的消息
                    other_vars = [v for v in vars_in_factor if v != target_var]

                    # 计算消息
                    msg = np.zeros(len(self.variables[target_var]['domain']))

                    # 枚举其他变量的所有赋值
                    def enumerate_assignments(idx, current_assignment):
                        if idx == len(other_vars):
                            # 计算因子值
                            for val_idx, val in enumerate(self.variables[target_var]['domain']):
                                full_assignment = {**current_assignment, target_var: val}
                                factor_val = factor_func(full_assignment)

                                # 乘以其他变量的消息
                                msg_factor = 1.0
                                for other_var in other_vars:
                                    other_val = full_assignment[other_var]
                                    other_val_idx = self.variables[other_var]['domain'].index(other_val)
                                    msg_factor *= self.messages[('v2f', other_var, factor_id)][other_val_idx]

                                msg[val_idx] += factor_val * msg_factor
                            return

                        var = other_vars[idx]
                        for val in self.variables[var]['domain']:
                            current_assignment[var] = val
                            enumerate_assignments(idx + 1, current_assignment)
                            del current_assignment[var]

                    enumerate_assignments(0, {})
                    self.messages[('f2v', factor_id, target_var)] = msg / msg.sum()

        # 计算边缘分布
        marginals = {}
        for var in self.variables:
            belief = np.ones(len(self.variables[var]['domain']))
            for factor in self.variables[var]['neighbors']:
                belief *= self.messages[('f2v', factor, var)]
            marginals[var] = belief / belief.sum()

        return marginals

# 示例: 简单的纠错码
def create_ldpc_code():
    """创建LDPC码的因子图"""
    fg = FactorGraph()

    # 变量（比特）
    for i in range(4):
        fg.add_variable(f'x{i}', [0, 1])

    # 校验因子（奇偶校验）
    def parity_check(assignment):
        """奇偶校验: 和为偶数时返回1，否则返回0"""
        return 1.0 if sum(assignment.values()) % 2 == 0 else 0.0

    fg.add_factor('f0', ['x0', 'x1', 'x2'], parity_check)
    fg.add_factor('f1', ['x1', 'x2', 'x3'], parity_check)

    # 观测因子（信道输出）
    def likelihood(x_idx, p_flip):
        """似然因子"""
        def factor(assignment):
            x = assignment[f'x{x_idx}']
            # 假设观测到0
            return 1 - p_flip if x == 0 else p_flip
        return factor

    fg.add_factor('obs0', ['x0'], likelihood(0, 0.1))
    fg.add_factor('obs1', ['x1'], likelihood(1, 0.1))
    fg.add_factor('obs2', ['x2'], likelihood(2, 0.1))
    fg.add_factor('obs3', ['x3'], likelihood(3, 0.1))

    return fg

# 运行和积算法
fg = create_ldpc_code()
marginals = fg.sum_product(max_iter=5)

print("边缘分布:")
for var, dist in marginals.items():
    print(f"{var}: P(0)={dist[0]:.4f}, P(1)={dist[1]:.4f}")

def analyze_convergence():
    """分析不同条件下的收敛率"""

    # 问题设置: 二次函数 f(x) = 0.5 * x^T A x - b^T x
    # 其中 A 的特征值在 [mu, L] 之间

    mu, L = 1, 100  # 强凸参数和平滑参数
    kappa = L / mu  # 条件数

    # 创建病态矩阵
    n = 10
    eigenvalues = np.linspace(mu, L, n)
    A = np.diag(eigenvalues)
    b = np.ones(n)

    # 最优解
    x_star = np.linalg.solve(A, b)
    f_star = 0.5 * x_star.T @ A @ x_star - b.T @ x_star

    def gradient_descent(x0, eta, n_iter):
        """梯度下降"""
        x = x0.copy()
        history = []

        for _ in range(n_iter):
            grad = A @ x - b
            x = x - eta * grad

            f = 0.5 * x.T @ A @ x - b.T @ x
            history.append(f - f_star)

        return x, history

    # 不同步长的比较
    x0 = np.random.randn(n)

    etas = [1/L, 2/(mu + L), 1/mu]
    labels = ['η=1/L (保守)', 'η=2/(μ+L) (最优)', 'η=1/μ (激进)']

    import matplotlib.pyplot as plt
    plt.figure(figsize=(10, 6))

    for eta, label in zip(etas, labels):
        _, history = gradient_descent(x0, eta, 200)
        plt.semilogy(history, label=label)

    # 理论收敛率
    T = np.arange(200)
    theory = 0.5 * L * np.linalg.norm(x0 - x_star)**2 / (T + 1)
    plt.semilogy(theory, '--', label='Theory: O(1/T)', alpha=0.5)

    plt.xlabel('Iteration')
    plt.ylabel('f(x) - f* (log scale)')
    plt.title('Gradient Descent Convergence')
    plt.legend()
    plt.grid(True)
    plt.show()

    print(f"条件数 κ = {kappa:.2f}")

analyze_convergence()

def sgd_convergence_analysis():
    """SGD收敛性分析"""

    np.random.seed(42)

    # 线性回归: min ||Ax - b||^2
    n, d = 1000, 20
    A = np.random.randn(n, d)
    x_true = np.random.randn(d)
    b = A @ x_true + 0.1 * np.random.randn(n)

    # 计算最优解
    x_star = np.linalg.lstsq(A, b, rcond=None)[0]

    # 强凸参数
    mu = 2 * np.min(np.linalg.eigvals(A.T @ A)) / n
    L = 2 * np.max(np.linalg.eigvals(A.T @ A)) / n

    def sgd(x0, lr_schedule, n_epochs, batch_size=1):
        """SGD"""
        x = x0.copy()
        history = []

        for epoch in range(n_epochs):
            indices = np.random.permutation(n)

            for i in range(0, n, batch_size):
                batch_idx = indices[i:i+batch_size]
                A_batch = A[batch_idx]
                b_batch = b[batch_idx]

                # 随机梯度
                grad = 2 * A_batch.T @ (A_batch @ x - b_batch) / len(batch_idx)

                # 学习率
                t = epoch * (n // batch_size) + i // batch_size
                lr = lr_schedule(t)

                x = x - lr * grad

            # 记录误差
            error = np.linalg.norm(x - x_star)**2
            history.append(error)

        return history

    x0 = np.zeros(d)
    n_epochs = 50

    # 不同学习率调度
    # 策略模式：字典值为lambda函数，通过schedules['decay'](t)调用对应的学习率计算公式
    schedules = {
        'constant': lambda t: 0.01,
        'decay': lambda t: 0.1 / (1 + 0.01 * t),
        'optimal': lambda t: 1 / (mu * (t + 1))
    }

    plt.figure(figsize=(10, 6))

    for name, schedule in schedules.items():
        history = sgd(x0, schedule, n_epochs)
        plt.semilogy(history, label=name)

    # 理论收敛率 O(1/T)
    T = np.arange(1, n_epochs + 1)
    theory = 1 / T
    plt.semilogy(theory / theory[0] * 10, '--', label='Theory: O(1/T)', alpha=0.5)

    plt.xlabel('Epoch')
    plt.ylabel('||x - x*||^2 (log scale)')
    plt.title('SGD Convergence Analysis')
    plt.legend()
    plt.grid(True)
    plt.show()

sgd_convergence_analysis()

def compare_optimization_methods():
    """比较不同优化方法的收敛速度"""

    # 病态二次函数
    n = 100
    kappa = 1000  # 条件数

    eigenvalues = np.linspace(1, kappa, n)
    A = np.diag(eigenvalues)
    b = np.ones(n)

    x_star = np.linalg.solve(A, b)

    def gradient_descent(x0, eta, n_iter):
        """普通梯度下降"""
        x = x0.copy()
        history = []

        for _ in range(n_iter):
            grad = A @ x - b
            x = x - eta * grad
            error = np.linalg.norm(x - x_star)
            history.append(error)

        return history

    def nesterov_acceleration(x0, eta, n_iter):
        """Nesterov加速梯度下降"""
        x = x0.copy()
        y = x0.copy()
        history = []

        for t in range(n_iter):
            x_new = y - eta * (A @ y - b)
            y = x_new + (t / (t + 3)) * (x_new - x)
            x = x_new

            error = np.linalg.norm(x - x_star)
            history.append(error)

        return history

    def conjugate_gradient(x0, n_iter):
        """共轭梯度法"""
        x = x0.copy()
        r = b - A @ x
        p = r.copy()
        history = []

        for _ in range(n_iter):
            Ap = A @ p
            alpha = np.dot(r, r) / np.dot(p, Ap)
            x = x + alpha * p
            r_new = r - alpha * Ap

            beta = np.dot(r_new, r_new) / np.dot(r, r)
            p = r_new + beta * p
            r = r_new

            error = np.linalg.norm(x - x_star)
            history.append(error)

        return history

    x0 = np.random.randn(n)
    n_iter = 200

    # 运行各种方法
    gd_history = gradient_descent(x0, 2/(1+kappa), n_iter)
    nesterov_history = nesterov_acceleration(x0, 1/kappa, n_iter)
    cg_history = conjugate_gradient(x0, n_iter)

    # 绘制结果
    plt.figure(figsize=(12, 6))

    plt.subplot(1, 2, 1)
    plt.semilogy(gd_history, label='Gradient Descent')
    plt.semilogy(nesterov_history, label='Nesterov Acceleration')
    plt.semilogy(cg_history, label='Conjugate Gradient')
    plt.xlabel('Iteration')
    plt.ylabel('||x - x*||')
    plt.title('Convergence Comparison')
    plt.legend()
    plt.grid(True)

    plt.subplot(1, 2, 2)
    # 前50次迭代
    plt.semilogy(gd_history[:50], label='Gradient Descent')
    plt.semilogy(nesterov_history[:50], label='Nesterov Acceleration')
    plt.semilogy(cg_history[:50], label='Conjugate Gradient')
    plt.xlabel('Iteration')
    plt.ylabel('||x - x*||')
    plt.title('First 50 Iterations')
    plt.legend()
    plt.grid(True)

    plt.tight_layout()
    plt.show()

    print(f"条件数: {kappa}")
    print(f"GD最终误差: {gd_history[-1]:.2e}")
    print(f"Nesterov最终误差: {nesterov_history[-1]:.2e}")
    print(f"CG最终误差: {cg_history[-1]:.2e}")

compare_optimization_methods()

class VarianceReductionMethods:
    """方差缩减方法实现"""

    def __init__(self, A, b):
        self.A = A
        self.b = b
        self.n = A.shape[0]

    def full_gradient(self, x):
        """计算全梯度"""
        return 2 * self.A.T @ (self.A @ x - self.b) / self.n

    def stochastic_gradient(self, x, i):
        """计算随机梯度"""
        a_i = self.A[i]
        return 2 * a_i * (a_i @ x - self.b[i])

    def svrg(self, x0, n_epochs, inner_iter, lr):
        """SVRG算法"""
        x = x0.copy()
        history = []

        for epoch in range(n_epochs):
            # 计算全梯度
            x_tilde = x.copy()
            mu_tilde = self.full_gradient(x_tilde)

            for t in range(inner_iter):
                # 随机采样
                i = np.random.randint(self.n)

                # 方差缩减梯度
                grad = (self.stochastic_gradient(x, i) -
                       self.stochastic_gradient(x_tilde, i) + mu_tilde)

                x = x - lr * grad

            error = np.linalg.norm(self.full_gradient(x))
            history.append(error)

        return x, history

    def saga(self, x0, n_epochs, lr):
        """SAGA算法"""
        x = x0.copy()
        history = []

        # 存储历史梯度
        grad_history = np.zeros((self.n, len(x0)))
        for i in range(self.n):
            grad_history[i] = self.stochastic_gradient(x, i)

        grad_avg = grad_history.mean(axis=0)

        for epoch in range(n_epochs):
            for _ in range(self.n):
                i = np.random.randint(self.n)

                # 新梯度
                grad_new = self.stochastic_gradient(x, i)

                # SAGA更新
                grad = grad_new - grad_history[i] + grad_avg
                x = x - lr * grad

                # 更新历史
                grad_avg += (grad_new - grad_history[i]) / self.n
                grad_history[i] = grad_new

            error = np.linalg.norm(self.full_gradient(x))
            history.append(error)

        return x, history

def test_variance_reduction():
    """测试方差缩减方法"""

    np.random.seed(42)

    # 生成数据
    n, d = 1000, 50
    A = np.random.randn(n, d)
    x_true = np.random.randn(d)
    b = A @ x_true + 0.1 * np.random.randn(n)

    vr = VarianceReductionMethods(A, b)
    x0 = np.zeros(d)

    # 运行SVRG
    x_svrg, history_svrg = vr.svrg(x0, n_epochs=20, inner_iter=n, lr=0.01)

    # 运行SAGA
    x_saga, history_saga = vr.saga(x0, n_epochs=20, lr=0.01)

    # 运行普通SGD
    def sgd(x0, n_epochs, lr):
        x = x0.copy()
        history = []

        for epoch in range(n_epochs):
            for _ in range(n):
                i = np.random.randint(n)
                grad = vr.stochastic_gradient(x, i)
                x = x - lr * grad

            error = np.linalg.norm(vr.full_gradient(x))
            history.append(error)

        return history

    history_sgd = sgd(x0, 20, 0.001)

    # 绘制结果
    plt.figure(figsize=(10, 6))
    plt.semilogy(history_sgd, label='SGD', alpha=0.7)
    plt.semilogy(history_svrg, label='SVRG', linewidth=2)
    plt.semilogy(history_saga, label='SAGA', linewidth=2)

    plt.xlabel('Epoch')
    plt.ylabel('Gradient Norm (log scale)')
    plt.title('Variance Reduction Methods')
    plt.legend()
    plt.grid(True)
    plt.show()

    print(f"SGD最终梯度范数: {history_sgd[-1]:.4f}")
    print(f"SVRG最终梯度范数: {history_svrg[-1]:.4f}")
    print(f"SAGA最终梯度范数: {history_saga[-1]:.4f}")

test_variance_reduction()