目录
2.1 线性回归与最小二乘
2.1.1 解析解与几何解释
线性回归模型的参数估计可归结为最小化残差平方和的目标函数。解析解的构造依赖于设计矩阵的QR分解或Cholesky分解，其中Cholesky分解利用对称正定矩阵的三角化性质，将正规方程的求解转化为前向与后向替换过程，数值稳定性显著优于直接矩阵求逆。从几何视角审视，最小二乘估计等价于响应向量在设计矩阵列空间上的正交投影，投影矩阵（Hat矩阵）的构造直接导出杠杆值与Cook距离等诊断统计量，为异常点检测提供解析基础。计算复杂度的权衡体现在直接求解与迭代方法的比较：基于稠密矩阵分解的方法具有O(ND2) 的复杂度，其中N 为样本量，D 为特征维度；当N≫D 时，采用Gram矩阵的Cholesky分解可显著降低计算开销，而条件数分析揭示了数据共线性对数值稳定性的影响。
2.1.2 概率视角与高斯噪声假设
在概率框架下，响应变量的高斯噪声假设导出似然函数的对数形式与最小二乘目标函数的等价性。极大似然估计在此设定下退化为残差平方和最小化。引入参数的先验分布后，最大后验估计自然地正则化为Ridge回归，其中高斯先验对应L2惩罚项，拉普拉斯先验对应L1惩罚项。后验预测分布的积分计算涉及高斯分布的卷积性质，解析解的获得依赖于对权重空间的边缘化，这为预测不确定性的量化提供了贝叶斯解释。
脚本1：线性回归数值实现与概率解释
Python
"""
脚本1：线性回归的数值方法与概率视角实现
涉及内容：
- 2.1.1.1 正规方程的Cholesky分解求解
- 2.1.1.2 Hat矩阵计算与杠杆值可视化
- 2.1.2.1 MLE与最小二乘等价性数值验证
- 2.1.2.2 MAP估计与Ridge回归的贝叶斯实现
使用方式：
python linear_regression_numerical.py
依赖：numpy, scipy, matplotlib, seaborn
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy import linalg
from scipy.stats import multivariate_normal
import seaborn as sns
# 设置样式
sns.set_palette("husl")
plt.style.use('seaborn-v0_8-whitegrid')
class LinearRegressionNumerical:
"""
线性回归的数值稳定实现
包含直接求逆、Cholesky分解、SVD三种求解策略
"""
def __init__(self, method='cholesky'):
self.method = method
self.coef_ = None
self.intercept_ = None
self.H = None # Hat矩阵
self.leverage_ = None
self.condition_number_ = None
def fit(self, X, y):
"""
拟合线性回归模型
实现复杂度O(ND^2)的Cholesky分解与直接O(D^3)方法对比
"""
n_samples, n_features = X.shape
# 增广设计矩阵 [1, X]
X_aug = np.column_stack([np.ones(n_samples), X])
# 计算Gram矩阵 A = X^T X
A = X_aug.T @ X_aug
b = X_aug.T @ y
# 记录条件数用于数值稳定性分析
self.condition_number_ = np.linalg.cond(A)
if self.method == 'direct':
# 直接求逆（数值不稳定，仅用于教学对比）
self.coef_ = np.linalg.inv(A) @ b
elif self.method == 'cholesky':
# Cholesky分解: A = L L^T
# 复杂度O(D^3/3)，但数值稳定性优于直接求逆
try:
L = linalg.cholesky(A, lower=True)
# 前向替换解 L y = b
y_temp = linalg.solve_triangular(L, b, lower=True)
# 后向替换解 L^T x = y
self.coef_ = linalg.solve_triangular(L.T, y_temp, lower=False)
except np.linalg.LinAlgError:
# 矩阵非正定，添加微小正则化
A_reg = A + 1e-8 * np.eye(n_features + 1)
L = linalg.cholesky(A_reg, lower=True)
y_temp = linalg.solve_triangular(L, b, lower=True)
self.coef_ = linalg.solve_triangular(L.T, y_temp, lower=False)
elif self.method == 'svd':
# SVD分解用于共线性处理
U, s, Vh = linalg.svd(X_aug, full_matrices=False)
# 伪逆求解
self.coef_ = Vh.T @ np.diag(1/s) @ U.T @ y
self.intercept_ = self.coef_[0]
self.coef_ = self.coef_[1:]
# 计算Hat矩阵 H = X (X^T X)^{-1} X^T
# 用于杠杆值与影响分析
X_aug_pinv = np.linalg.pinv(X_aug)
self.H = X_aug @ X_aug_pinv
self.leverage_ = np.diag(self.H)
return self
def predict(self, X):
"""预测响应值"""
return X @ self.coef_ + self.intercept_
def bayesian_posterior(self, X, y, alpha=1.0, beta=1.0):
"""
2.1.2.2 贝叶斯线性回归的后验分布计算
alpha: 噪声精度先验 (1/sigma^2)
beta: 权重先验精度 (1/lambda^2)
"""
n_samples, n_features = X.shape
X_aug = np.column_stack([np.ones(n_samples), X])
# 后验精度矩阵 S^{-1} = beta I + alpha X^T X
S_inv = beta * np.eye(n_features + 1) + alpha * X_aug.T @ X_aug
# Cholesky分解用于稳定计算
L = linalg.cholesky(S_inv, lower=True)
# 后验均值 mu = alpha S X^T y
temp = linalg.solve_triangular(L, X_aug.T @ y, lower=True)
mu = alpha * linalg.solve_triangular(L.T, temp, lower=False)
# 计算后验协方差 S = (S^{-1})^{-1}
S = linalg.cho_solve((L, True), np.eye(n_features + 1))
return mu, S
def posterior_predictive(self, X_new, X_train, y_train, alpha=1.0, beta=1.0):
"""
2.1.2.3 后验预测分布计算
返回预测分布的均值与方差
"""
mu_post, S_post = self.bayesian_posterior(X_train, y_train, alpha, beta)
X_new_aug = np.column_stack([np.ones(X_new.shape[0]), X_new])
# 预测均值
y_pred_mean = X_new_aug @ mu_post
# 预测方差: sigma^2 + X_new S X_new^T
pred_var = 1.0/alpha + np.sum(X_new_aug @ S_post * X_new_aug, axis=1)
return y_pred_mean, pred_var
def demo_numerical_stability():
"""
演示Cholesky分解与直接求逆的数值稳定性差异
特别是当设计矩阵具有较高条件数时
"""
np.random.seed(42)
# 构造病态矩阵（高度共线性）
n_samples = 100
x_base = np.random.randn(n_samples)
X = np.column_stack([x_base, x_base + 1e-10 * np.random.randn(n_samples)])
y = 2 * x_base + 1 + 0.1 * np.random.randn(n_samples)
methods = ['direct', 'cholesky', 'svd']
results = {}
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for idx, method in enumerate(methods):
model = LinearRegressionNumerical(method=method)
model.fit(X, y)
results[method] = {
'coef': model.coef_,
'intercept': model.intercept_,
'condition': model.condition_number_,
'residual': np.linalg.norm(y - model.predict(X))
}
# 可视化杠杆值
axes[idx].scatter(range(len(model.leverage_)), model.leverage_,
c=model.leverage_, cmap='viridis', alpha=0.6)
axes[idx].axhline(y=2*(X.shape[1]+1)/len(X), color='r',
linestyle='--', label='Threshold: 2(p+1)/n')
axes[idx].set_title(f'{method.upper()} Method\n'
f'Cond: {model.condition_number_:.2e}\n'
f'Residual: {results[method]["residual"]:.4f}')
axes[idx].set_xlabel('Sample Index')
axes[idx].set_ylabel('Leverage (Hat Diagonal)')
axes[idx].legend()
axes[idx].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('numerical_stability_comparison.png', dpi=150, bbox_inches='tight')
plt.show()
print("数值稳定性对比:")
for method, res in results.items():
print(f"{method:10s}: Coef=[{res['coef'][0]:.6f}, {res['coef'][1]:.6f}], "
f"Intercept={res['intercept']:.6f}, Cond={res['condition']:.2e}")
def demo_bayesian_approach():
"""
演示贝叶斯线性回归与后验预测分布
对比MLE与MAP估计的差异
"""
np.random.seed(42)
# 生成数据
n_samples = 50
X = np.linspace(-3, 3, n_samples).reshape(-1, 1)
y_true = np.sin(X).ravel()
y = y_true + 0.3 * np.random.randn(n_samples)
# 多项式特征扩展以展示正则化效果
from sklearn.preprocessing import PolynomialFeatures
poly = PolynomialFeatures(degree=7, include_bias=False)
X_poly = poly.fit_transform(X)
# MLE估计（无正则化，对应alpha->inf, beta->0）
model_mle = LinearRegressionNumerical(method='svd')
model_mle.fit(X_poly, y)
y_mle = model_mle.predict(X_poly)
# MAP估计（Ridge回归，对应有限alpha, beta）
alpha = 10.0 # 噪声精度
beta = 0.1 # 权重先验精度（正则化强度）
model_map = LinearRegressionNumerical(method='cholesky')
mu_post, S_post = model_map.bayesian_posterior(X_poly, y, alpha=alpha, beta=beta)
# 后验预测
y_pred_mean, y_pred_var = model_map.posterior_predictive(
X_poly, X_poly, y, alpha=alpha, beta=beta
)
y_pred_std = np.sqrt(y_pred_var)
# 采样后验权重以展示不确定性
n_posterior_samples = 100
post_samples = multivariate_normal.rvs(mean=mu_post, cov=S_post,
size=n_posterior_samples)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左图：MLE vs MAP对比
axes[0].scatter(X, y, c='black', s=30, alpha=0.6, label='Observations')
axes[0].plot(X, y_true, 'g-', linewidth=2, label='True Function')
axes[0].plot(X, y_mle, 'r--', linewidth=2, label='MLE (No Regularization)')
axes[0].plot(X, y_pred_mean, 'b-', linewidth=2, label='MAP (Ridge)')
axes[0].fill_between(X.ravel(),
y_pred_mean - 2*y_pred_std,
y_pred_mean + 2*y_pred_std,
alpha=0.3, color='blue', label='95% Predictive Interval')
axes[0].set_title('MLE vs MAP Estimation\n(High-degree Polynomial)')
axes[0].legend()
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].grid(True, alpha=0.3)
# 右图：后验权重分布
axes[1].bar(range(len(mu_post)), mu_post, yerr=2*np.sqrt(np.diag(S_post)),
capsize=5, alpha=0.7, color='steelblue', edgecolor='black')
axes[1].set_title('Posterior Mean of Weights\n(with 2σ Uncertainty)')
axes[1].set_xlabel('Weight Index')
axes[1].set_ylabel('Posterior Value')
axes[1].grid(True, alpha=0.3, axis='y')
plt.tight_layout()
plt.savefig('bayesian_linear_regression.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"\n贝叶斯分析结果:")
print(f"噪声精度参数 alpha: {alpha}")
print(f"权重先验精度 beta: {beta} (等价于Ridge的lambda={alpha/beta})")
print(f"后验权重均值: {mu_post[:4]} ...") # 只显示前4个
if __name__ == "__main__":
demo_numerical_stability()
demo_bayesian_approach() 
2.2 正则化与稀疏性
2.2.1 L2正则化与Ridge回归
Ridge回归通过引入L2惩罚项处理设计矩阵的奇异或近奇异情形。奇异值分解揭示了Ridge估计量对最小二乘解的收缩效应：滤波系数以σj2+λσj2 的形式几何衰减，其中σj 为设计矩阵的奇异值。这种收缩机制有效地改善了条件数，使病态问题的数值求解趋于稳定。核化Ridge回归的对偶形式将计算复杂度从O(D3) 转换为O(N3) ，在高维特征空间（D≫N ）中实现计算效率的逆转，为后续核方法的引入奠定算法基础。
2.2.2 L1正则化与Lasso回归
Lasso回归采用L1惩罚诱导参数稀疏性，优化问题不再具有解析解，需借助数值方法求解。软阈值算子作为近端映射的具体形式，构成了迭代收缩阈值算法（ISTA）的核心步骤。坐标下降法通过逐维度更新策略规避Hessian矩阵的显式构造，活跃集策略进一步加速收敛，仅对非零系数的子集执行迭代更新。最小角回归（LARS）则提供了整条正则化路径的高效计算，其几何解释与最小二乘的前向选择密切相关，在每一步选择与当前残差具有最大相关性的预测变量，并沿等角方向更新系数，直至所有变量进入模型或相关系数降至零。
2.2.3 弹性网络与结构化稀疏
弹性网络结合L1与L2惩罚，在保持稀疏性的同时处理组效应，避免高度相关变量选择时的随机排除现象。稀疏组Lasso扩展至多层级结构，通过分组惩罚项实现变量的层次化选择。快速迭代收缩阈值算法（FISTA）引入Nesterov动量加速策略，将收敛速率从O(1/k) 提升至O(1/k2) ，在大规模稀疏优化问题中显著减少迭代次数。
脚本2：正则化方法与稀疏优化实现
Python
"""
脚本2：正则化方法与稀疏优化算法实现
涉及内容：
- 2.2.1.1 Ridge回归的SVD视角与滤波系数分析
- 2.2.2.1 软阈值算子与坐标下降法实现Lasso
- 2.2.2.3 LARS算法完整路径求解
- 2.2.3.3 FISTA加速与Nesterov动量实现
使用方式：
python regularization_sparse_optimization.py
依赖：numpy, matplotlib, scikit-learn (仅用于数据生成对比)
"""
import numpy as np
import matplotlib.pyplot as plt
from itertools import cycle
import warnings
class RidgeRegressionSVD:
"""
2.2.1.1 Ridge回归的SVD实现与滤波系数分析
展示收缩效应的几何衰减特性
"""
def __init__(self, lambdas=[0.01, 0.1, 1.0, 10.0]):
self.lambdas = lambdas
self.results_ = {}
def fit(self, X, y):
"""
使用SVD分解实现Ridge回归
X = U S V^T, w_ridge = V (S^2 + lambda I)^{-1} S U^T y
"""
# 中心化处理
self.mean_X = X.mean(axis=0)
self.mean_y = y.mean()
X_centered = X - self.mean_X
y_centered = y - self.mean_y
# SVD分解
U, s, Vh = np.linalg.svd(X_centered, full_matrices=False)
V = Vh.T
for lam in self.lambdas:
# 滤波系数: s_j^2 / (s_j^2 + lambda)
filter_coef = s**2 / (s**2 + lam)
# Ridge系数
w = V @ np.diag(filter_coef / s) @ U.T @ y_centered
self.results_[lam] = {
'coef': w,
'filter': filter_coef,
'singular_values': s,
'effective_df': np.sum(filter_coef) # 有效自由度
}
return self
def plot_shrinkage(self):
"""可视化奇异值收缩效应"""
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左图：滤波系数随lambda变化
s = self.results_[self.lambdas[0]]['singular_values']
for lam in self.lambdas:
axes[0].plot(range(1, len(s)+1),
self.results_[lam]['filter'],
'o-', label=f'λ={lam}', markersize=6)
axes[0].set_xlabel('Singular Value Index')
axes[0].set_ylabel('Filter Coefficient σ²/(σ²+λ)')
axes[0].set_title('Geometric Decay of Filter Coefficients\n(Ridge Shrinkage Effect)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[0].set_ylim(0, 1.05)
# 右图：有效自由度与lambda关系
eff_dfs = [self.results_[lam]['effective_df'] for lam in self.lambdas]
axes[1].semilogx(self.lambdas, eff_dfs, 'bo-', markersize=8, linewidth=2)
axes[1].set_xlabel('Regularization Parameter λ')
axes[1].set_ylabel('Effective Degrees of Freedom')
axes[1].set_title('Effective DoF vs Regularization Strength')
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('ridge_shrinkage_analysis.png', dpi=150, bbox_inches='tight')
plt.show()
class LassoCoordinateDescent:
"""
2.2.2.1 坐标下降法实现Lasso
结合软阈值算子与活跃集策略
"""
def __init__(self, lambda_reg=0.1, max_iter=1000, tol=1e-4,
active_set=True):
self.lambda_reg = lambda_reg
self.max_iter = max_iter
self.tol = tol
self.active_set = active_set
self.coef_path_ = []
def _soft_threshold(self, z, gamma):
"""软阈值算子: sign(z) * max(|z| - gamma, 0)"""
return np.sign(z) * np.maximum(np.abs(z) - gamma, 0)
def fit(self, X, y):
"""
坐标下降法求解Lasso
活跃集策略：仅更新非零系数以加速收敛
"""
n_samples, n_features = X.shape
# 标准化
self.X_mean = X.mean(axis=0)
self.X_std = X.std(axis=0) + 1e-8
self.y_mean = y.mean()
X_std = (X - self.X_mean) / self.X_std
y_std = y - self.y_mean
# 初始化
w = np.zeros(n_features)
residuals = y_std.copy()
# 预计算X的列范数（用于标准化lambda）
x_norms = np.sum(X_std**2, axis=0)
for iteration in range(self.max_iter):
w_old = w.copy()
if self.active_set and iteration > 0:
# 活跃集策略：仅更新当前非零的系数
active_indices = np.where(w != 0)[0]
if len(active_indices) == 0:
active_indices = range(n_features)
else:
active_indices = range(n_features)
# 坐标下降循环
for j in active_indices:
# 残差扣除第j个特征的影响
residuals += X_std[:, j] * w[j]
# 偏相关
rho_j = np.dot(X_std[:, j], residuals)
# 软阈值更新
w[j] = self._soft_threshold(rho_j / x_norms[j],
self.lambda_reg / x_norms[j])
# 残差更新
residuals -= X_std[:, j] * w[j]
# 记录路径
if iteration % 10 == 0:
self.coef_path_.append(w.copy())
# 收敛判断
if np.max(np.abs(w - w_old)) < self.tol:
break
self.coef_ = w / self.X_std # 反标准化
self.intercept_ = self.y_mean - np.dot(self.X_mean, self.coef_)
self.n_iter_ = iteration + 1
return self
def predict(self, X):
return X @ self.coef_ + self.intercept_
class LARS:
"""
2.2.2.3 最小角回归（LARS）实现
计算完整正则化路径，等角方向更新
"""
def __init__(self):
self.coef_path_ = None
self.lambda_path_ = None
def fit(self, X, y):
"""
LARS算法实现
基于Efron et al. (2004) Annals of Statistics
"""
n_samples, n_features = X.shape
# 标准化
X = X / np.sqrt(np.sum(X**2, axis=0))
# 初始化
beta = np.zeros(n_features)
mu = np.zeros(n_samples) # 当前预测值
c = X.T @ y # 当前相关性
# 路径记录
self.coef_path_ = [beta.copy()]
self.lambda_path_ = [np.max(np.abs(c))]
active_set = []
signs = np.zeros(n_features)
n_iter = min(n_features, n_samples - 1)
for _ in range(n_iter):
# 寻找最大相关性变量
c_abs = np.abs(c)
C_max = np.max(c_abs)
if C_max < 1e-12:
break
# 加入活跃集
j = np.argmax(c_abs)
if j not in active_set:
active_set.append(j)
signs[j] = np.sign(c[j])
# 计算等角方向
X_A = X[:, active_set] * signs[active_set]
G_A = X_A.T @ X_A # Gram矩阵
try:
# 求解等角方向 A_A = (G_A)^{-1} 1_A
ones_A = np.ones(len(active_set))
w_A = np.linalg.solve(G_A, ones_A)
# 归一化
A_A = 1 / np.sqrt(np.sum(w_A))
u_A = X_A @ w_A * A_A
# 计算步长
a = X.T @ u_A
# 计算加入新变量的gamma
gamma_candidates = []
for j_inactive in range(n_features):
if j_inactive not in active_set:
cc = c[j_inactive]
aa = a[j_inactive]
# 方程 (C_max - cc) / (A_A - aa) 和 (C_max + cc) / (A_A + aa)
if abs(A_A - aa) > 1e-12:
gamma_candidates.append((C_max - cc) / (A_A - aa))
if abs(A_A + aa) > 1e-12:
gamma_candidates.append((C_max + cc) / (A_A + aa))
if len(gamma_candidates) > 0:
gamma = np.min([g for g in gamma_candidates if g > 1e-12])
else:
gamma = C_max / A_A
except np.linalg.LinAlgError:
# 数值问题，使用伪逆
break
# 更新系数和预测值
for idx, var_idx in enumerate(active_set):
beta[var_idx] += gamma * signs[var_idx] * w_A[idx] * A_A
mu += gamma * u_A
c = X.T @ (y - mu)
self.coef_path_.append(beta.copy())
self.lambda_path_.append(np.max(np.abs(c)) if len(c) > 0 else 0)
self.coef_path_ = np.array(self.coef_path_)
self.lambda_path_ = np.array(self.lambda_path_)
self.coef_ = beta
return self
class FISTA:
"""
2.2.3.3 FISTA加速近端梯度下降
用于Elastic Net优化，收敛率O(1/k^2)
"""
def __init__(self, lambda_reg=0.1, alpha=0.5, max_iter=500, tol=1e-6):
"""
alpha: Elastic Net混合参数 (0=纯Ridge, 1=纯Lasso)
"""
self.lambda_reg = lambda_reg
self.alpha = alpha # L1比例
self.max_iter = max_iter
self.tol = tol
def _soft_threshold(self, x, gamma):
return np.sign(x) * np.maximum(np.abs(x) - gamma, 0)
def fit(self, X, y):
"""
FISTA算法实现
基于Nesterov加速: t_{k+1} = (1 + sqrt(1 + 4*t_k^2)) / 2
"""
n_samples, n_features = X.shape
# 初始化
w = np.zeros(n_features)
w_old = w.copy()
t = 1.0
# Lipschitz常数估计（谱范数）
L = np.linalg.norm(X, ord=2)**2 + 1e-8
for k in range(self.max_iter):
# 外推点
y_k = w + ((t - 1) / (t + 1)) * (w - w_old) # Nesterov动量
# 梯度计算（平滑部分：最小二乘 + Ridge）
grad = X.T @ (X @ y_k - y) / n_samples + \
(1 - self.alpha) * self.lambda_reg * y_k
# 近端映射（L1部分）
w_new = self._soft_threshold(
y_k - grad / L,
self.alpha * self.lambda_reg / L
)
# 收敛判断
if np.linalg.norm(w_new - w) < self.tol:
break
# 更新动量参数
t_new = (1 + np.sqrt(1 + 4 * t**2)) / 2
w_old = w.copy()
w = w_new
t = t_new
self.coef_ = w
self.n_iter_ = k + 1
return self
def demo_regularization_paths():
"""
对比Ridge、Lasso、LARS的正则化路径
"""
np.random.seed(42)
# 构造稀疏数据
n_samples, n_features = 50, 10
true_coef = np.array([5, -4, 3, 0, 0, 0, 0, 2, 0, 0])
X = np.random.randn(n_samples, n_features)
y = X @ true_coef + 0.5 * np.random.randn(n_samples)
# Ridge路径（使用SVD）
lambdas = np.logspace(-3, 2, 50)
ridge_coefs = []
U, s, Vh = np.linalg.svd(X, full_matrices=False)
for lam in lambdas:
filter_coef = s / (s**2 + lam)
w = Vh.T @ np.diag(filter_coef) @ U.T @ y
ridge_coefs.append(w)
ridge_coefs = np.array(ridge_coefs)
# LARS路径
lars = LARS()
lars.fit(X, y)
# 绘制对比
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# Ridge路径
for i in range(n_features):
axes[0].plot(lambdas, ridge_coefs[:, i], label=f'Feature {i+1}' if i < 3 else '')
axes[0].set_xscale('log')
axes[0].set_xlabel('Regularization Parameter λ')
axes[0].set_ylabel('Coefficient Value')
axes[0].set_title('Ridge Regression Path\n(Smooth Shrinkage)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# LARS路径
for i in range(n_features):
axes[1].plot(lars.lambda_path_, lars.coef_path_[:, i],
marker='o', markersize=3, label=f'Feature {i+1}' if i < 3 else '')
axes[1].set_xlabel('Max Correlation (λ equivalent)')
axes[1].set_ylabel('Coefficient Value')
axes[1].set_title('LARS Path\n(Piecewise Linear with Knots)')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('regularization_paths.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"LARS路径长度: {len(lars.lambda_path_)}")
print(f"LARS最终非零系数数量: {np.sum(lars.coef_ != 0)}")
if __name__ == "__main__":
# 生成测试数据
np.random.seed(42)
X = np.random.randn(100, 8)
true_w = np.array([3, -2, 1.5, 0, 0, 0, 0.5, 0])
y = X @ true_w + 0.3 * np.random.randn(100)
# Ridge SVD分析
print("=== Ridge SVD分析 ===")
ridge = RidgeRegressionSVD(lambdas=[0.001, 0.01, 0.1, 1.0, 10.0])
ridge.fit(X, y)
ridge.plot_shrinkage()
# Lasso坐标下降
print("\n=== Lasso坐标下降 ===")
lasso = LassoCoordinateDescent(lambda_reg=0.2, active_set=True)
lasso.fit(X, y)
print(f"迭代次数: {lasso.n_iter_}")
print(f"非零系数: {np.sum(lasso.coef_ != 0)} / {len(lasso.coef_)}")
print(f"估计系数: {np.round(lasso.coef_, 3)}")
print(f"真实系数: {true_w}")
# LARS
print("\n=== LARS算法 ===")
lars = LARS()
lars.fit(X, y)
print(f"路径步数: {len(lars.lambda_path_)}")
# FISTA
print("\n=== FISTA加速 ===")
fista = FISTA(lambda_reg=0.2, alpha=0.8, max_iter=1000)
fista.fit(X, y)
print(f"FISTA迭代次数: {fista.n_iter_}")
print(f"非零系数: {np.sum(fista.coef_ != 0)}")
# 路径可视化
demo_regularization_paths()运行结果
=== Ridge SVD分析 ===
=== Lasso坐标下降 ===
迭代次数: 6
非零系数: 7 / 8
估计系数: [ 3.008e+00 -1.952e+00 1.498e+00 -1.500e-02 -3.000e-03 0.000e+00
4.620e-01 1.800e-02]
真实系数: [ 3. -2. 1.5 0. 0. 0. 0.5 0. ]
=== LARS算法 ===
路径步数: 9
=== FISTA加速 ===
FISTA迭代次数: 1000
非零系数: 4
LARS路径长度: 11
LARS最终非零系数数量: 4
2.3 分类模型：从线性到概率
2.3.1 感知机与逻辑回归
感知机算法通过迭代更新权重向量寻找分离超平面，Novikoff收敛定理保证了在数据线性可分时算法在有限步内终止，间隔（Margin）的大小直接影响收敛速度的上界。对数几率回归将分类问题纳入概率框架，通过极大似然估计优化对数似然函数，其目标函数为凸函数保证全局最优解的存在。迭代重加权最小二乘（IRLS）算法将非线性优化问题转化为一系列加权最小二乘子问题，每步迭代基于当前参数估计计算工作响应与权重矩阵，利用Newton-Raphson方法实现二次收敛速率，Hessian矩阵的解析形式为XTWX ，其中W 为对角的方差权重矩阵。
2.3.2 多类分类与Softmax回归
Softmax回归将二分类逻辑回归扩展至多类情形，输出层的概率分布通过指数归一化获得。数值稳定性实践中需采用Log-Sum-Exp技巧避免指数运算的上溢或下溢。多项分布的广义线性模型形式与交叉熵损失的等价性通过最大似然估计建立，多类逻辑回归的Hessian矩阵具有分块结构，利用Kronecker积或稀疏结构可优化高维参数优化中的存储与计算。
2.3.3 概率生成模型
高斯判别分析假设类条件分布服从多元高态，通过共享或独立协方差矩阵分别导出线性判别分析（LDA）与二次判别分析（QDA）。协方差共享机制减少了待估参数数量，在样本量有限时降低过拟合风险。朴素贝叶斯基于特征条件独立性假设，通过拉普拉斯平滑处理零概率问题。生成式与判别式模型的渐近误差分析表明，生成式模型在模型设定正确时收敛更快，而判别式模型对模型误设具有更强的鲁棒性。
脚本3：分类模型与概率实现
"""
脚本3：分类模型与概率生成/判别式实现
涉及内容：
- 2.3.1.3 IRLS算法Newton-Raphson实现与Hessian计算
- 2.3.2.1 Softmax数值稳定性技巧（Log-Sum-Exp）
- 2.3.2.3 多类Hessian分块结构优化
- 2.3.3.1 GDA/QDA协方差共享推导实现
使用方式：
python classification_models.py
依赖：numpy, matplotlib, scipy
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import block_diag
from scipy.special import logsumexp
import warnings
class LogisticRegressionIRLS:
"""
2.3.1.3 迭代重加权最小二乘（IRLS）实现
Newton-Raphson方法，显式Hessian计算
"""
def __init__(self, max_iter=100, tol=1e-6, fit_intercept=True):
self.max_iter = max_iter
self.tol = tol
self.fit_intercept = fit_intercept
self.coef_ = None
self.intercept_ = None
self.convergence_history_ = []
def _sigmoid(self, z):
"""数值稳定的sigmoid实现"""
# 防止溢出
return np.where(z >= 0,
1 / (1 + np.exp(-z)),
np.exp(z) / (1 + np.exp(z)))
def fit(self, X, y):
"""
IRLS算法实现
基于Fisher评分：w_new = w + (X^T W X)^{-1} X^T (y - p)
其中W是权重矩阵，对角线元素为 p(1-p)
"""
n_samples, n_features = X.shape
# 添加截距项
if self.fit_intercept:
X_aug = np.column_stack([np.ones(n_samples), X])
else:
X_aug = X
# 初始化
w = np.zeros(X_aug.shape[1])
for i in range(self.max_iter):
# 线性预测
z = X_aug @ w
p = self._sigmoid(z)
# 工作响应（Working Response）
# z_working = z + (y - p) / (p * (1 - p))
# 但直接计算权重矩阵W更稳定
W_diag = p * (1 - p)
W_diag = np.maximum(W_diag, 1e-8) # 防止除零
# 梯度（Score Function）: X^T (y - p)
gradient = X_aug.T @ (y - p)
# Hessian: -X^T W X
# 使用Cholesky分解求解线性系统避免直接求逆
try:
# 构造加权设计矩阵 sqrt(W) * X
X_weighted = X_aug * np.sqrt(W_diag)[:, np.newaxis]
# 正规方程: (X^T W X) delta = X^T (y - p)
A = X_weighted.T @ X_weighted
# 添加微小正则化确保正定性
A += 1e-8 * np.eye(A.shape[0])
# 求解更新方向
L = np.linalg.cholesky(A)
temp = np.linalg.solve(L, gradient)
delta = np.linalg.solve(L.T, temp)
except np.linalg.LinAlgError:
# 回退到伪逆
H = -X_aug.T @ np.diag(W_diag) @ X_aug
delta = -np.linalg.pinv(H) @ gradient
# 更新参数
w_new = w + delta
# 记录收敛历史
loss = -np.sum(y * np.log(p + 1e-10) + (1-y) * np.log(1-p + 1e-10))
self.convergence_history_.append(loss)
# 收敛判断
if np.linalg.norm(delta) < self.tol:
break
w = w_new
if self.fit_intercept:
self.intercept_ = w[0]
self.coef_ = w[1:]
else:
self.coef_ = w
self.n_iter_ = i + 1
return self
def predict_proba(self, X):
z = X @ self.coef_ + (self.intercept_ if self.fit_intercept else 0)
return self._sigmoid(z)
def predict(self, X):
return (self.predict_proba(X) >= 0.5).astype(int)
class SoftmaxRegression:
"""
2.3.2.1 Softmax回归实现
包含Log-Sum-Exp数值稳定性技巧
"""
def __init__(self, max_iter=100, tol=1e-6, reg_lambda=0.01):
self.max_iter = max_iter
self.tol = tol
self.reg_lambda = reg_lambda
def _softmax_stable(self, Z):
"""
数值稳定的Softmax实现
减去每行的最大值防止指数上溢
"""
# Log-Sum-Exp技巧
Z_max = np.max(Z, axis=1, keepdims=True)
exp_Z = np.exp(Z - Z_max)
return exp_Z / np.sum(exp_Z, axis=1, keepdims=True)
def _compute_hessian_block(self, X, P, n_classes):
"""
2.3.2.3 计算多类Hessian的分块结构
Hessian是块对角矩阵，每块为 -X^T W_k X
总维度 (n_classes * n_features) x (n_classes * n_features)
"""
n_samples, n_features = X.shape
# 对于多类逻辑回归，Hessian是Kronecker积结构
# 这里实现简化的块对角近似用于教学演示
hessian_blocks = []
for k in range(n_classes):
# 第k类的权重矩阵
W_k = P[:, k] * (1 - P[:, k])
H_k = -X.T @ np.diag(W_k) @ X - self.reg_lambda * np.eye(n_features)
hessian_blocks.append(H_k)
# 构造块对角Hessian（近似）
return block_diag(*hessian_blocks)
def fit(self, X, y):
"""
多类逻辑回归的Newton-Raphson实现
参数W为 n_features x n_classes 矩阵
"""
n_samples, n_features = X.shape
n_classes = len(np.unique(y))
# One-hot编码
Y_onehot = np.eye(n_classes)[y]
# 初始化参数 (增广X以包含截距)
X_aug = np.column_stack([np.ones(n_samples), X])
n_aug_features = X_aug.shape[1]
# 向量化参数: (n_aug_features * n_classes) x 1
W = np.zeros((n_aug_features, n_classes))
for iteration in range(self.max_iter):
# 线性预测 n_samples x n_classes
Z = X_aug @ W
# 概率预测（数值稳定）
P = self._softmax_stable(Z)
# 梯度: X_aug^T (P - Y_onehot) + reg
G = X_aug.T @ (P - Y_onehot) + self.reg_lambda * W
# 向量化梯度用于收敛判断
g_vec = G.ravel()
# Hessian计算与求解（分块Newton-Raphson）
# 使用IRLS风格的加权最小二乘近似
W_old = W.copy()
for k in range(n_classes):
# 第k类的IRLS更新
W_k_diag = P[:, k] * (1 - P[:, k])
# 加权正规方程
X_weighted = X_aug * np.sqrt(W_k_diag)[:, np.newaxis]
try:
# Cholesky求解
A = X_weighted.T @ X_weighted + self.reg_lambda * np.eye(n_aug_features)
L = np.linalg.cholesky(A)
temp = np.linalg.solve(L, G[:, k])
update_k = np.linalg.solve(L.T, temp)
W[:, k] -= 0.5 * update_k # 步长0.5确保收敛
except:
# 梯度下降回退
W[:, k] -= 0.01 * G[:, k]
# 收敛判断
if np.linalg.norm(W - W_old) < self.tol:
break
self.coef_ = W[1:, :]
self.intercept_ = W[0, :]
self.classes_ = np.unique(y)
return self
def predict_proba(self, X):
n_samples = X.shape[0]
X_aug = np.column_stack([np.ones(n_samples), X])
W = np.row_stack([self.intercept_, self.coef_])
Z = X_aug @ W
return self._softmax_stable(Z)
def predict(self, X):
proba = self.predict_proba(X)
return self.classes_[np.argmax(proba, axis=1)]
class GaussianDiscriminantAnalysis:
"""
2.3.3.1 高斯判别分析（GDA）
实现LDA（共享协方差）与QDA（独立协方差）
"""
def __init__(self, shared_covariance=True):
self.shared_covariance = shared_covariance
self.classes_ = None
self.phi_ = None # 类先验
self.mu_ = None # 类均值
self.sigma_ = None # 共享协方差（LDA）
self.sigma_k_ = None # 独立协方差（QDA）
def fit(self, X, y):
"""
拟合GDA模型
共享协方差时为LDA，独立时为QDA
"""
self.classes_ = np.unique(y)
n_samples, n_features = X.shape
n_classes = len(self.classes_)
# 估计类先验
self.phi_ = np.array([np.mean(y == k) for k in self.classes_])
# 估计类均值
self.mu_ = np.array([X[y == k].mean(axis=0) for k in self.classes_])
if self.shared_covariance:
# LDA: 共享协方差矩阵（加权平均）
sigma = np.zeros((n_features, n_features))
for idx, k in enumerate(self.classes_):
X_k = X[y == k]
diff = X_k - self.mu_[idx]
sigma += self.phi_[idx] * (diff.T @ diff / len(X_k))
self.sigma_ = sigma
# 预计算逆矩阵与行列式
self.sigma_inv_ = np.linalg.pinv(self.sigma_)
self.sigma_det_ = np.linalg.det(self.sigma_)
else:
# QDA: 独立协方差矩阵
self.sigma_k_ = []
self.sigma_inv_k_ = []
self.sigma_det_k_ = []
for idx, k in enumerate(self.classes_):
X_k = X[y == k]
diff = X_k - self.mu_[idx]
sigma_k = (diff.T @ diff) / len(X_k)
self.sigma_k_.append(sigma_k)
self.sigma_inv_k_.append(np.linalg.pinv(sigma_k))
self.sigma_det_k_.append(np.linalg.det(sigma_k))
return self
def _log_gaussian_prob(self, X, mu, sigma_inv, sigma_det):
"""计算多元高斯对数概率"""
diff = X - mu
# log N(x; mu, sigma) = -0.5 * (log|sigma| + (x-mu)^T sigma^{-1} (x-mu) + d*log(2*pi))
mahalanobis = np.sum(diff @ sigma_inv * diff, axis=1)
log_prob = -0.5 * (np.log(sigma_det + 1e-8) + mahalanobis + X.shape[1] * np.log(2*np.pi))
return log_prob
def predict_log_proba(self, X):
"""预测对数概率"""
n_samples = X.shape[0]
log_probs = np.zeros((n_samples, len(self.classes_)))
for idx, k in enumerate(self.classes_):
if self.shared_covariance:
# LDA: 共享协方差，线性决策边界
log_likelihood = self._log_gaussian_prob(
X, self.mu_[idx], self.sigma_inv_, self.sigma_det_
)
else:
# QDA: 独立协方差，二次决策边界
log_likelihood = self._log_gaussian_prob(
X, self.mu_[idx],
self.sigma_inv_k_[idx],
self.sigma_det_k_[idx]
)
log_probs[:, idx] = np.log(self.phi_[idx] + 1e-10) + log_likelihood
# 使用log-sum-exp技巧归一化
log_probs -= logsumexp(log_probs, axis=1, keepdims=True)
return log_probs
def predict(self, X):
log_probs = self.predict_log_proba(X)
return self.classes_[np.argmax(log_probs, axis=1)]
def demo_classification_models():
"""
演示二分类与多分类模型的对比
"""
np.random.seed(42)
# 构造二分类数据
n_samples = 200
X_class0 = np.random.multivariate_normal([0, 0], [[1, 0.5], [0.5, 1]], n_samples//2)
X_class1 = np.random.multivariate_normal([3, 3], [[1, -0.5], [-0.5, 1]], n_samples//2)
X_binary = np.vstack([X_class0, X_class1])
y_binary = np.array([0]*(n_samples//2) + [1]*(n_samples//2))
# 训练Logistic回归
log_reg = LogisticRegressionIRLS(max_iter=50)
log_reg.fit(X_binary, y_binary)
# 训练LDA与QDA
lda = GaussianDiscriminantAnalysis(shared_covariance=True)
lda.fit(X_binary, y_binary)
qda = GaussianDiscriminantAnalysis(shared_covariance=False)
qda.fit(X_binary, y_binary)
# 可视化决策边界
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
models = [
(log_reg, 'Logistic Regression (IRLS)', 'coef_'),
(lda, 'LDA (Shared Covariance)', 'mu_'),
(qda, 'QDA (Independent Covariance)', 'mu_')
]
for idx, (model, title, attr) in enumerate(models):
# 绘制数据点
axes[idx].scatter(X_binary[y_binary==0, 0], X_binary[y_binary==0, 1],
c='blue', alpha=0.5, label='Class 0')
axes[idx].scatter(X_binary[y_binary==1, 0], X_binary[y_binary==1, 1],
c='red', alpha=0.5, label='Class 1')
# 绘制决策边界
x_min, x_max = X_binary[:, 0].min() - 1, X_binary[:, 0].max() + 1
y_min, y_max = X_binary[:, 1].min() - 1, X_binary[:, 1].max() + 1
xx, yy = np.meshgrid(np.linspace(x_min, x_max, 100),
np.linspace(y_min, y_max, 100))
grid = np.c_[xx.ravel(), yy.ravel()]
if isinstance(model, LogisticRegressionIRLS):
Z = model.predict_proba(grid).reshape(xx.shape)
else:
Z = np.exp(model.predict_log_proba(grid))[:, 1].reshape(xx.shape)
axes[idx].contourf(xx, yy, Z, levels=20, alpha=0.3, cmap='RdBu')
axes[idx].contour(xx, yy, Z, levels=[0.5], colors='k', linewidths=2)
axes[idx].set_title(title)
axes[idx].legend()
axes[idx].set_xlabel('Feature 1')
axes[idx].set_ylabel('Feature 2')
axes[idx].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('classification_boundaries.png', dpi=150, bbox_inches='tight')
plt.show()
# 多分类演示
print("\n=== 多分类Softmax回归 ===")
X_multi = np.random.randn(300, 4)
true_w = np.array([[2, -1, 0.5], [1, 2, -1], [-0.5, 1, 1.5], [0, 0, 1]])
logits = X_multi @ true_w
y_multi = np.argmax(logits + 0.5*np.random.randn(300, 3), axis=1)
softmax = SoftmaxRegression(max_iter=200)
softmax.fit(X_multi, y_multi)
acc = np.mean(softmax.predict(X_multi) == y_multi)
print(f"Softmax训练准确率: {acc:.4f}")
print(f"Softmax迭代次数: {softmax.n_iter_}")
if __name__ == "__main__":
demo_classification_models()
2.4 基函数与非线性扩展
2.4.1 多项式回归与样条
多项式回归通过范德蒙德矩阵构造特征映射，其条件数随阶数增长呈指数恶化，导致数值不稳定。正交多项式基（Legendre、Chebyshev）通过Gram-Schmidt正交化或递推关系构造，在区间上关于特定权函数正交，显著改善设计矩阵的条件数。三次样条插值要求分段三次多项式在节点处保持连续的二阶导数，自然边界条件设定端点二阶导数为零，将插值问题转化为三对角线性系统的求解，利用Thomas算法可在O(N) 时间内完成。B-样条基函数通过de Boor-Cox递推公式构造，具备局部支撑性质，即每个基函数仅在有限区间内非零，这一特性使B-样条在局部数据修改时具有计算稳定性，de Boor算法提供了高效的递归计算框架。
2.4.2 径向基函数与核方法初步
径向基函数网络以输入到中心点的径向距离作为非线性映射，中心点的选择通过K-Means聚类或正交最小二乘（OLS）优化。高斯过程回归建立了函数空间视角与权重空间视角的等价性，通过Mercer定理保证核函数对应无穷维特征空间的合法映射，协方差函数（核函数）直接定义函数空间的相似性度量，规避了显式特征映射的计算。
脚本4：基函数方法与核方法实现
"""
脚本4：基函数方法与核方法实现
涉及内容：
- 2.4.1.1 正交多项式基（Legendre/Chebyshev）与条件数改善
- 2.4.1.3 B-样条的de Boor-Cox递推实现
- 2.4.2.1 RBF网络中心点选择（K-Means/OLS）
- 2.4.2.2 高斯过程的权重空间与函数空间视角
使用方式：
python basis_functions_kernels.py
依赖：numpy, matplotlib, scipy
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import legendre, chebyt
from scipy.linalg import solve_banded
from scipy.spatial.distance import cdist
from sklearn.cluster import KMeans
class OrthogonalPolynomialBasis:
"""
2.4.1.1 正交多项式基实现
对比范德蒙德矩阵与正交多项式的条件数
"""
def __init__(self, degree=10, poly_type='legendre'):
self.degree = degree
self.poly_type = poly_type
def _vandermonde(self, x):
"""标准范德蒙德矩阵（病态）"""
return np.vander(x, self.degree + 1, increasing=True)
def _legendre_basis(self, x):
"""
Legendre多项式基（区间[-1,1]正交）
通过递推关系构造数值稳定
"""
# 映射到[-1, 1]
x_scaled = 2 * (x - x.min()) / (x.max() - x.min()) - 1
basis = np.zeros((len(x), self.degree + 1))
for d in range(self.degree + 1):
# 使用scipy的Legendre函数
coeffs = [0] * d + [1]
P_d = legendre(d)
basis[:, d] = P_d(x_scaled)
return basis
def _chebyshev_basis(self, x):
"""
Chebyshev多项式基（第一类）
在[-1,1]上关于权函数(1-x^2)^{-1/2}正交
"""
x_scaled = 2 * (x - x.min()) / (x.max() - x.min()) - 1
basis = np.zeros((len(x), self.degree + 1))
for d in range(self.degree + 1):
T_d = chebyt(d)
basis[:, d] = T_d(x_scaled)
return basis
def fit_transform(self, x):
"""根据类型返回相应的基"""
if self.poly_type == 'vandermonde':
return self._vandermonde(x)
elif self.poly_type == 'legendre':
return self._legendre_basis(x)
elif self.poly_type == 'chebyshev':
return self._chebyshev_basis(x)
else:
raise ValueError(f"Unknown type: {self.poly_type}")
def compare_condition_numbers(self, x):
"""对比不同多项式基的条件数"""
types = ['vandermonde', 'legendre', 'chebyshev']
cond_numbers = []
for t in types:
self.poly_type = t
X_basis = self.fit_transform(x)
cond = np.linalg.cond(X_basis)
cond_numbers.append(cond)
return dict(zip(types, cond_numbers))
class BSplineBasis:
"""
2.4.1.3 B-样条基函数实现
de Boor-Cox递推公式
"""
def __init__(self, degree=3, n_knots=10):
self.degree = degree
self.n_knots = n_knots
def _cox_deboor_recursion(self, x, i, k, knots):
"""
de Boor-Cox递推公式实现
B_{i,0}(x) = 1 if knots[i] <= x < knots[i+1], else 0
B_{i,k}(x) = w_{i,k}(x) B_{i,k-1}(x) + (1-w_{i+1,k}(x)) B_{i+1,k-1}(x)
其中 w_{i,k}(x) = (x - knots[i]) / (knots[i+k] - knots[i])
"""
if k == 0:
# 0次基函数（指示函数）
if i < len(knots) - 1:
return np.where((x >= knots[i]) & (x < knots[i+1]), 1.0, 0.0)
return np.zeros_like(x)
# 递归计算
B_im1 = self._cox_deboor_recursion(x, i, k-1, knots)
B_ip1 = self._cox_deboor_recursion(x, i+1, k-1, knots)
# 计算权重（处理除零）
left_denom = knots[i+k] - knots[i]
right_denom = knots[i+k+1] - knots[i+1]
w_left = np.zeros_like(x)
w_right = np.zeros_like(x)
if abs(left_denom) > 1e-10:
w_left = (x - knots[i]) / left_denom
if abs(right_denom) > 1e-10:
w_right = (knots[i+k+1] - x) / right_denom
return w_left * B_im1 + w_right * B_ip1
def fit_transform(self, x):
"""
生成B-样条基矩阵
使用重复的边界节点实现自然边界条件
"""
x_min, x_max = x.min(), x.max()
# 创建节点向量（Clamped，重复边界）
interior_knots = np.linspace(x_min, x_max, self.n_knots)
# 边界节点重复degree+1次以实现插值性质
knots = np.concatenate([
[x_min] * (self.degree + 1),
interior_knots,
[x_max] * (self.degree + 1)
])
n_basis = len(knots) - self.degree - 1
basis_matrix = np.zeros((len(x), n_basis))
for i in range(n_basis):
basis_matrix[:, i] = self._cox_deboor_recursion(x, i, self.degree, knots)
# 处理最后一个点的包含性（右闭区间）
basis_matrix[x == x_max, n_basis-1] = 1.0
return basis_matrix, knots
def fit(self, x, y):
"""拟合B-样条回归"""
self.x_range_ = (x.min(), x.max())
B, self.knots_ = self.fit_transform(x)
# 最小二乘拟合（可使用稀疏矩阵优化）
self.coef_ = np.linalg.lstsq(B, y, rcond=None)[0]
return self
def predict(self, x):
B, _ = self.fit_transform(x)
return B @ self.coef_
class RBFnetwork:
"""
2.4.2.1 RBF网络实现
中心点选择：K-Means与正交最小二乘（OLS）
"""
def __init__(self, n_centers=10, method='kmeans', sigma=None):
self.n_centers = n_centers
self.method = method
self.sigma = sigma # 高斯RBF宽度
self.centers_ = None
self.weights_ = None
def _rbf_kernel(self, x, centers):
"""高斯径向基函数"""
dists = cdist(x, centers, 'euclidean')
if self.sigma is None:
# 启发式设置sigma
self.sigma = np.median(dists) / np.sqrt(2)
return np.exp(-dists**2 / (2 * self.sigma**2))
def fit(self, x, y):
"""
拟合RBF网络
method='kmeans': 使用K-Means选择中心点
method='ols': 正交最小二乘逐步选择中心点
"""
if self.method == 'kmeans':
# K-Means选择中心点
kmeans = KMeans(n_clusters=self.n_centers, random_state=42, n_init=10)
kmeans.fit(x.reshape(-1, 1) if x.ndim == 1 else x)
self.centers_ = kmeans.cluster_centers_
elif self.method == 'ols':
# 正交最小二乘（简化版）
# 从候选池中选择使残差减少最大的中心点
candidates = x[np.random.choice(len(x), min(50, len(x)), replace=False)]
if x.ndim == 1:
candidates = candidates.reshape(-1, 1)
selected = []
residuals = y.copy()
phi_full = self._rbf_kernel(x, candidates)
for _ in range(self.n_centers):
# 选择解释最多方差的中心点
correlations = np.abs(phi_full.T @ residuals)
best_idx = np.argmax(correlations)
selected.append(candidates[best_idx])
# 更新残差（Gram-Schmidt正交化思想）
if len(selected) == 1:
phi_selected = phi_full[:, best_idx:best_idx+1]
else:
phi_selected = self._rbf_kernel(x, np.array(selected))
# 投影到已选空间
proj = phi_selected @ np.linalg.pinv(phi_selected) @ residuals
residuals = y - proj
self.centers_ = np.array(selected).reshape(-1, 1) if x.ndim == 1 else np.array(selected)
# 拟合权重（Ridge正则化确保稳定）
Phi = self._rbf_kernel(x, self.centers_)
# 添加微小正则化
self.weights_ = np.linalg.lstsq(Phi, y, rcond=1e-6)[0]
self.phi_matrix_ = Phi
return self
def predict(self, x):
Phi = self._rbf_kernel(x, self.centers_)
return Phi @ self.weights_
class GaussianProcessRegression:
"""
2.4.2.2 高斯过程回归（RBF核）
展示函数空间视角（无需显式权重）
"""
def __init__(self, length_scale=1.0, noise_var=0.1):
self.length_scale = length_scale
self.noise_var = noise_var
def _rbf_kernel(self, x1, x2):
"""平方指数核（高斯核）"""
sqdist = cdist(x1, x2, 'sqeuclidean')
return np.exp(-0.5 * sqdist / self.length_scale**2)
def fit(self, X, y):
"""
GP回归的函数空间视角
计算后验分布 p(f|X,y) ~ N(mu, Sigma)
mu = K(X,X) (K(X,X) + sigma^2 I)^{-1} y
"""
self.X_train_ = X
self.y_train_ = y
# 训练协方差矩阵
K = self._rbf_kernel(X, X)
K_noise = K + self.noise_var * np.eye(len(X))
# 使用Cholesky分解求解（稳定）
try:
L = np.linalg.cholesky(K_noise)
self.alpha_ = np.linalg.solve(L.T, np.linalg.solve(L, y))
except:
# 如果矩阵接近奇异，添加微小抖动
L = np.linalg.cholesky(K_noise + 1e-6 * np.eye(len(X)))
self.alpha_ = np.linalg.solve(L.T, np.linalg.solve(L, y))
self.K_inv_solver_ = (L, L.T) # 保存求解器用于预测
return self
def predict(self, X, return_std=True):
"""
预测：后验高斯过程
mu* = K(X*, X) alpha
var* = K(X*, X*) - K(X*, X) (K(X,X) + sigma^2 I)^{-1} K(X, X*)
"""
K_star = self._rbf_kernel(X, self.X_train_)
mu = K_star @ self.alpha_
if return_std:
# 计算预测方差
K_star_star = self._rbf_kernel(X, X)
# 使用已分解的矩阵求解
L, Lt = self.K_inv_solver_
v = np.linalg.solve(L, K_star.T)
var = np.diag(K_star_star) - np.sum(v**2, axis=0)
std = np.sqrt(np.maximum(var, 0))
return mu, std
return mu
def demo_basis_functions():
"""
演示正交多项式、B-样条、RBF与GP的对比
"""
np.random.seed(42)
# 生成复杂函数数据
x = np.linspace(-3, 3, 100)
y_true = np.sin(x) * np.exp(-x**2/10) + 0.5 * np.cos(2*x)
y = y_true + 0.2 * np.random.randn(len(x))
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
# 1. 多项式条件数对比
degrees = range(5, 21)
vand_cond = []
leg_cond = []
cheb_cond = []
for d in degrees:
ortho = OrthogonalPolynomialBasis(degree=d)
conds = ortho.compare_condition_numbers(x)
vand_cond.append(conds['vandermonde'])
leg_cond.append(conds['legendre'])
cheb_cond.append(conds['chebyshev'])
axes[0, 0].semilogy(degrees, vand_cond, 'o-', label='Vandermonde', markersize=6)
axes[0, 0].semilogy(degrees, leg_cond, 's-', label='Legendre', markersize=6)
axes[0, 0].semilogy(degrees, cheb_cond, '^-', label='Chebyshev', markersize=6)
axes[0, 0].set_xlabel('Polynomial Degree')
axes[0, 0].set_ylabel('Condition Number (log scale)')
axes[0, 0].set_title('Condition Number Comparison\n(Orthogonal vs Monomial Basis)')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)
# 2. B-样条拟合
bspline = BSplineBasis(degree=3, n_knots=8)
bspline.fit(x, y)
x_fine = np.linspace(-3, 3, 300)
y_bspline = bspline.predict(x_fine)
axes[0, 1].scatter(x, y, c='gray', alpha=0.5, s=30, label='Data')
axes[0, 1].plot(x_fine, y_true, 'g-', linewidth=2, alpha=0.7, label='True Function')
axes[0, 1].plot(x_fine, y_bspline, 'r--', linewidth=2, label='B-Spline (cubic)')
axes[0, 1].scatter(bspline.knots_[4:-4],
np.interp(bspline.knots_[4:-4], x_fine, y_bspline),
c='red', s=100, marker='x', label='Interior Knots', zorder=5)
axes[0, 1].set_title('Cubic B-Spline Regression\n(de Boor-Cox Recursion)')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
# 3. B-样条基函数可视化
x_vis = np.linspace(-3, 3, 300)
B_vis, _ = bspline.fit_transform(x_vis)
for i in range(min(6, B_vis.shape[1])):
axes[0, 2].plot(x_vis, B_vis[:, i], alpha=0.7, label=f'$B_{{{i},3}}$')
axes[0, 2].set_title('B-Spline Basis Functions\n(Local Support Property)')
axes[0, 2].set_xlabel('x')
axes[0, 2].legend(ncol=2, fontsize=8)
axes[0, 2].grid(True, alpha=0.3)
# 4. RBF网络对比（K-Means vs OLS中心选择）
X = x.reshape(-1, 1)
rbf_kmeans = RBFnetwork(n_centers=8, method='kmeans')
rbf_kmeans.fit(X, y)
rbf_ols = RBFnetwork(n_centers=8, method='ols')
rbf_ols.fit(X, y)
y_rbf_kmeans = rbf_kmeans.predict(X)
y_rbf_ols = rbf_ols.predict(X)
axes[1, 0].scatter(x, y, c='gray', alpha=0.5, s=30)
axes[1, 0].plot(x, y_true, 'g-', linewidth=2, alpha=0.7, label='True Function')
axes[1, 0].plot(x, y_rbf_kmeans, 'b-', linewidth=2, label='RBF (K-Means Centers)')
axes[1, 0].scatter(rbf_kmeans.centers_,
rbf_kmeans.predict(rbf_kmeans.centers_),
c='blue', s=100, marker='x', zorder=5)
axes[1, 0].set_title('RBF Network: K-Means Center Selection')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
axes[1, 1].scatter(x, y, c='gray', alpha=0.5, s=30)
axes[1, 1].plot(x, y_true, 'g-', linewidth=2, alpha=0.7, label='True Function')
axes[1, 1].plot(x, y_rbf_ols, 'purple', linewidth=2, label='RBF (OLS Centers)')
axes[1, 1].scatter(rbf_ols.centers_,
rbf_ols.predict(rbf_ols.centers_),
c='purple', s=100, marker='x', zorder=5)
axes[1, 1].set_title('RBF Network: OLS Center Selection')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
# 5. 高斯过程回归（不确定性量化）
gp = GaussianProcessRegression(length_scale=0.5, noise_var=0.1)
gp.fit(X, y)
mu, std = gp.predict(x_fine.reshape(-1, 1), return_std=True)
axes[1, 2].scatter(x, y, c='gray', alpha=0.5, s=30, label='Observations')
axes[1, 2].plot(x_fine, y_true, 'g-', linewidth=2, alpha=0.7, label='True Function')
axes[1, 2].plot(x_fine, mu, 'b-', linewidth=2, label='GP Mean')
axes[1, 2].fill_between(x_fine, mu - 2*std, mu + 2*std,
alpha=0.3, color='blue', label='95% Confidence')
axes[1, 2].set_title('Gaussian Process Regression\n(Function Space Perspective)')
axes[1, 2].legend()
axes[1, 2].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('basis_functions_comparison.png', dpi=150, bbox_inches='tight')
plt.show()
# 打印数值分析结果
print("=== 正交多项式条件数分析 ===")
print(f"10阶范德蒙德条件数: {vand_cond[5]:.2e}")
print(f"10阶Legendre条件数: {leg_cond[5]:.2e}")
print(f"10阶Chebyshev条件数: {cheb_cond[5]:.2e}")
print(f"\n条件数改善比: {vand_cond[5]/leg_cond[5]:.2f}x")
if __name__ == "__main__":
demo_basis_functions()以上脚本提供了第2章监督学习核心内容的完整实现，涵盖从基础线性回归到高级非线性方法的数值技术细节。每个实现均包含计算复杂度优化、数值稳定性处理及可视化验证，可直接用于教学演示与算法原型开发。
1065
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
