简介：围绕投影追踪与投影寻踪回归实现的Scikit-Learn兼容估计器源码包，面向机器学习研究者、统计建模人员及想熟悉自定义估计器开发流程的Python开发者。项目基于Friedman与Stuetzle的Projection Pursuit Regression经典算法，提供支持多元估计与降维的回归器，以及将回归思想迁移到单变量投影表示上的分类器，适合处理高维数据的可解释建模与特征探索。包内共23个文件，以12个Python源码文件为核心，覆盖核心估计器、对比基准、训练过程可视化与测试，以及可视化脚本；辅助配置包含Travis CI、Coveralls、Sphinx、PyTest及PyPI/GitHub Pages部署等工程化文件，并遵循PEP 8与Scikit-Learn检查规范，可作为开发自定义估计器的完整精简示例。压缩包仅31KB，结构紧凑，便于快速通读和改造。目前已有149人学习，适合希望把投影追踪算法落地为sklearn组件的进阶读者。
1. 投影追踪：给高维数据一个「看得出结构」的回归与分类方案
搞过特征数量超过 30 的数据集的人都懂，模型调参走到一半，最让人头疼的不是损失函数不收敛，而是你根本说不清模型在拿哪些特征做判断。这个项目实现了一套基于投影追踪（Projection Pursuit）的回归与分类模块：把高维数据沿着若干个最优方向投影到低维空间，再在每个方向上做非线性拟合，从而让「特征太多、结构太隐」的数据集重新变得可解释、可预测。我拿到手第一反应是，它比直接堆一个多层神经网络更快收敛，而且能把每个方向对应的特征权重完整交给你看。适合谁？被高维非线性建模折腾过、想在 scikit-learn 生态里换一条路走的人。
2. 项目脉络与技术栈：一个机器学习模块该有的工程底座
2.1 它在补什么短板：高维非线性建模的另一种路径
高维数据常见的死法有两种。第一种是特征之间根本不是线性关系，线性回归不管怎么加交叉项都追不上真实结构；第二种是维度一高，无核的 KNN、朴素贝叶斯这类简单模型直接失效，而核方法又得为核函数的选择和带宽交一笔「玄学费」。投影追踪回归（Projection Pursuit Regression，PPR）走的是完全不同的路：它假设模型可以写成若干个方向函数之和
f(x) = g_1(a_1^T x) + g_2(a_2^T x) + ... + g_m(a_m^T x)
其中 a_m 是单位投影方向，g_m 是该方向上的岭函数（ridge function）。模型的表达能力来自两个部分：方向 a_m 负责挑选特征组合，岭函数 g_m 负责在这个组合上做非线性变换。比起神经网络的黑匣子，这套结构能直接把方向和函数曲线拿出来检查，这也是我一开始愿意拆这个源码的原因。
这个项目把 PPR 和基于投影方向的单变量分类都实现了，并且整体套用了 scikit-learn 的 Estimator 风格：有 fit、predict、score，能直接塞进 GridSearchCV 和 Pipeline。也就是说它不是一份贴在文档里的公式，而是一个拿到就能训练、能评估、能部署的模块。对要解决实际问题的人来说，这点比算法本身更值钱：你不需要自己处理交叉验证、样本划分、评分函数这些边角料。
2.2 技术栈盘点：一个 ML 模块的工程化底座
从依赖和工程配置可以很清楚地看到作者是按一个「正经开源项目」的标准来写的，不是临时脚本。数值部分以 NumPy 为主，所有矩阵运算、方向迭代、距离计算都直接落在 numpy 数组上；
回归的岭函数拟合里用到了 SciPy 的窗口查询与多项式拟合工具；最外层接口完全对齐 scikit-learn，包括检查输入维度的 check_array、模拟 parameter grid 的 get_params/set_params 约定。
| 依赖 | 用途 |
|---|---|
| numpy | 投影方向迭代、矩阵运算、随机重启 |
| scipy | 局部加权回归的距离窗口与多项式拟合 |
| scikit-learn | 接口规范、评分器、交叉验证与 Pipeline 对接 |
工程侧也有配套：测试目录完整，CI 流水线在每次提交后跑全量测试并把覆盖率上报，文档目录里放的是 API 说明和示例图生成脚本。这套搭配对二次开发非常友好——你改一行算法逻辑，CI 会立刻告诉你有没有把回归精度或分类准确率弄坏，覆盖率数字则会暴露哪些分支你压根没测到。我个人的习惯是，拿到源码先看 tests 目录，测试设计得细不细，直接决定了这个项目敢不敢在 production 里用。
2.3 目录结构：从 clone 到摸熟骨架需要知道的文件
整个仓库的骨架并不大，但每个文件的位置都有讲究：
projection_pursuit/ ├── projection_pursuit/ │ ├── __init__.py │ ├── _ppr.py # 多元投影追踪回归核心 │ ├── _classifier.py # 基于投影方向的单变量分类 │ ├── _smoother.py # 岭函数拟合：局部加权回归 │ └── _optimizer.py # 方向寻优：多起点重启+交替优化 ├── tests/ │ ├── test_ppr.py │ ├── test_classifier.py │ └── test_smoother.py ├── docs/ │ ├── conf.py │ └── index.rst ├── examples/ │ └── plot_ppr_demo.py ├── setup.py └── README.md其中_optimizer.py是最值得花时间读的：PPR 的拟合质量几乎全看方向寻优做得好不好，后面算法部分也会集中讲它。_smoother.py是岭函数的具体实现，改动它可以换掉整个非线性拟合策略。_classifier.py则是在回归核心之上套了一层分类损失，把连续响应改造成类别判定，两者共用同一套方向寻优逻辑，代码复用做得干净。
3. 核心算法拆解：多元投影追踪回归的完整实现链路
3.1 方向寻优：残差平方和是我唯一信得过的目标
PPR 的拟合过程用交替优化：固定岭函数时更新方向，固定方向时更新岭函数。项目里的做法是，先把当前轮残差算出来，然后找一个使残差平方和下降最多的投影方向。方向寻优的难点在于目标函数对方向参数是非凸的，直接梯度下降极容易掉进局部最优点，所以源码里做了多起点随机重启。
def _fit_one_term(self, X, residual): best = None for _ in range(self.n_restarts): direction = self._random_unit_vector(X.shape[1]) direction, ridge = self._alternating_optimize(X, residual, direction) score = self._term_score(residual, ridge) if best is None or score < best[0]: best = (score, direction, ridge) return best[1], best[2]这段代码的核心是_alternating_optimize：它会用当前方向算投影值，拟合岭函数，再反过来用岭函数对方向求梯度更新方向，循环若干圈直到收敛。direction始终被归一化成单位向量，避免优化过程中方向长度漫无目的地增长。多起点重启的作用相当于给非凸问题买了多重后悔药，n_restarts越大，越有机会找到更优的方向，但训练时间线性上涨。
参数说明：n_restarts默认一般在 10 到 30 之间，数据维度越高越应该加大；_term_score用的是该方向上拟合后的残差平方和，分数越低代表这一项吃掉的信息越多。如果发现训练结果随随机种子变化很大，第一反应不是去调学习率，而是把n_restarts调大。
3.2 岭函数：局部加权回归让每个投影方向都有表达力
方向定下来之后，每个样本在这个方向上会得到一个投影值 z，接下来要在 z 和残差之间拟合一个非线性函数 g。项目默认用局部加权回归（lowess 风格）来做这件事：对每个点，取它周围的邻域样本，用加权二次多项式拟合出该点的预测值。
def _lowess_smooth(z, y, bandwidth=0.3): n = len(z) fitted = np.zeros(n) for i, z0 in enumerate(z): dist = np.abs(z - z0) window = dist <= np.quantile(dist, bandwidth) if window.sum() < 5: window[np.argsort(dist)[:5]] = True d_max = dist[window].max() weight = (1 - (dist[window] / (d_max + 1e-12)) ** 2) ** 2 coeffs = np.polyfit(z[window], y[window], 2, w=weight) fitted[i] = np.polyval(coeffs, z0) return fitted这里bandwidth表示分位数窗口：0.3 意味着每个点只看离它最近的 30% 样本。权重是 tricube 核的平方，距离越近贡献越大；np.polyfit(..., 2)拟合的是局部二次曲线，比一次线性拟合更能抓住弯曲结构。window.sum() < 5的兜底是为了处理边界点或样本极稀疏的区域，强制至少取 5 个邻居，防止polyfit因为矩阵病态直接崩溃。
调这个参数对结果影响很直接。带宽太大，岭函数会趋向直线，模型退化成近似线性回归；带宽太小，每个点只靠身边三五个样本撑腰，拟合曲线剧烈震荡，训练集误差很低但验证集一塌糊涂。我一般先从 0.3 起步，再根据验证集残差调整。
3.3 主循环：逐项扣残差与早停策略
有了方向和岭函数，剩下的就是一层层把信息从残差里剥出来。主循环逻辑很简单：拟合一项，预测，扣掉残差，继续下一项。
def fit(self, X, y): X = self._validate_input(X) residual = y.copy() self.terms_ = [] for _ in range(self.max_terms): direction, ridge = self._fit_one_term(X, residual) pred = ridge(np.dot(X, direction)) self.terms_.append((direction, ridge)) residual -= pred if np.linalg.norm(residual) < self.tol * np.linalg.norm(y): break return self每一轮生成的(direction, ridge)会被保存下来，预测阶段把所有项的岭函数输出加总就是最终预测值。tol是相对残差阈值，默认设成 1e-3 左右；如果某轮拟合完后残差范数已经小于初始残差范数的千分之一，再继续加项只是过拟合训练集，所以直接提前退出。
参数max_terms控制模型复杂度。实际调参时我通常把这个值从小往大试，同时盯着残差下降曲线：如果加第 5 项之后残差几乎不再下降，说明结构已经抓完了，再多加只是让模型记住噪声。这里的早停策略已经帮你挡掉了大半过度拟合的风险。
4. 单变量分类：从回归输出到类别判定的落地实现
4.1 从连续响应到类别：把 PPR 输出改造成分类器
分类部分复用了回归的方向寻优逻辑，只是在最后加了一层类别判定。对于二分类，正常思路是想直接最小化分类损失，但分类损失对方向参数求导不友好，所以项目里做了一个近似处理：先拟合一个光滑的连续函数逼近类别标签，再用阈值把连续输出切成类别。
def decision_function(self, X): score = np.zeros(X.shape[0]) for direction, ridge in self.terms_: score += ridge(np.dot(X, direction)) return score def predict(self, X): return (self.decision_function(X) > self.threshold).astype(int)decision_function的形式和回归预测完全一样，区别只在最后一步：predict用threshold把连续分数映射成 0/1。threshold默认是 0.5，但遇到类别不平衡的数据时，这个值往往需要挪到训练集的正样本分数中位数附近。
多分类的处理方式则是对每一类做一次一对多拟合，每个类别拥有自己的一组投影方向和岭函数。预测时取所有类别里决策分数最高那个作为输出。这种做法实现简单，而且能和 scikit-learn 的predict_proba风格对齐，缺点是会稍微加重训练开销。
4.2 评估与交叉验证：别拿准确率当唯一指标
单变量分类数据里最常见的场景是正负样本比例悬殊，此时准确率会产生严重误导：99% 的负样本，模型全预测负类也能拿到 0.99 的准确率。项目代码里的score方法默认用的是准确率，但我在实际用的时候会换成 ROC-AUC 或 F1。
def score(self, X, y): from sklearn.metrics import roc_auc_score try: return roc_auc_score(y, self.decision_function(X)) except ValueError: # 单类别输入时 AUC 无法计算，退回准确率 return super().score(X, y)AUC 不依赖固定阈值，能更真实地反映方向寻优的质量。如果你发现 AUC 很高但准确率一般，说明模型结构是能分开两个类别的，只是threshold选得不好，此时调整阈值而不是重训模型。如果 AUC 本身就在 0.5 附近晃，那说明投影方向没有抓到区分性的特征组合，应该加max_terms或者检查特征预处理。
4.3 超参数怎么调：方向数、重启次数与带宽的搭配
| 参数 | 取值范围建议 | 影响 |
|---|---|---|
| max_terms | 1~8 | 太小欠拟合，太大过拟合 |
| n_restarts | 10~30 | 越大越稳，训练越慢 |
| bandwidth | 0.2~0.4 | 太小震荡，太大趋于线性 |
| threshold | 按类别分布调 | 只影响分类判定，不影响方向质量 |
调参顺序也有讲究。先固定bandwidth=0.3和n_restarts=20，从max_terms=1开始往回看残差曲线；再调宽带；最后如果训练结果依然不稳定，才去加大重启次数。这个顺序能避免同时动多个参数导致不知道是谁在起作用。网格搜索可以用，但方向寻优本身有随机性，所以同一组参数不同 seed 跑出来方差会很大，务必在交叉验证里固定随机种子，否则你看到的效果差距其实只是噪声。
5. 避坑与工程化：测试、CI、文档的常见翻车现场
5.1 投影方向震荡，换个 seed 结果就翻脸
现象：同一份训练数据，连着跑三次，三次的验证集误差差出 5 个百分点，投影方向权重也完全对不上。
原因：方向寻优的目标函数非凸，局部最优点到处都是。单次启动的梯度下降几乎不可能保证收敛到全局最优，随机种子一变就落进另一个局部坑里。
解决：把n_restarts从默认值往上加，并且在交叉验证和网格搜索里固定随机种子。常见做法是给模块加一个random_state参数，所有随机方向初始化都通过这个种子生成，这样每次实验可复现。注意n_restarts提升带来的收益有上限，到 30 次以后基本平稳，没必要无限堆算力。
5.2 岭函数在数据边界处突然起飞
现象：训练集整体误差正常，但把测试集样本按投影值排序后，处于两端极值的样本误差异常大，回归曲线在边界出现翘尾巴。
原因：局部加权回归在边界点能找到的邻居少，窗口被迫扩大后拟合出的二次曲线外推能力弱；再加上强制最少 5 个邻居的兜底逻辑，边界点的邻居很可能来自单侧，曲线斜率被带偏。
解决：先检查输入特征是否做了标准化，投影值的分布是否严重偏斜。其次可以适当调大bandwidth，让边界点能拉到更远的邻居；如果问题依旧，把岭函数换成带约束的样条拟合，比如在两端施加零二阶导数约束，这是我在实际项目中把默认 smoother 换掉的主要原因。
5.3 覆盖率卡死在 70% 附近上不去
现象：CI 流水线里的覆盖率统计一直在 70% 上下浮动，新增测试也补不上去，点开报告发现_optimizer.py里的失败分支几乎没被走到。
原因：测试全在 happy path 上。多起点重启逻辑里「几次重启全部失败」的分支、岭函数里window.sum() < 5的兜底分支、以及tol早停分支，都只有特定数据形状才能触发。写测试的人没专门构造这些边界样本。
解决：针对每个if分支单独构造数据集。想触发兜底逻辑就用样本量小于 5 的小批量数据；想触发早停就用一条线性关系非常强的数据，让第一项拟合完残差就低于阈值；想测重启失败就用全零特征列，让方向寻优始终停在原点。那之后我的覆盖率一般能推到 85% 以上。
5.4 CI 构建昨天还绿今天红成一片
现象：上午提交时全量测试通过，下午新拉了一次依赖就大面积报错，错误信息集中在scipy.linalg和numpy.polyfit的接口上。
原因：项目依赖没有锁版本，scipy发新版后改了部分子模块的函数签名或弃用警告升级成错误。测试本身没变，环境变了。
解决：在 CI 流水线里把依赖版本固定到具体版本区间，同时按 Python 版本做矩阵测试。另一个容易漏的是文档构建：示例脚本examples/plot_ppr_demo.py会随文档构建执行，只要示例里用到的 API 变动，文档构建就会失败。检查顺序永远是先看报错堆栈在哪一步，再确认依赖锁没锁住，最后看是不是示例代码过期。
6. 二次开发指南：把投影追踪模块改造进自己的项目
6.1 把默认平滑器换成样条基
_smoother.py里的_lowess_smooth是独立函数，接口是(z, y, bandwidth) -> fitted。我可以直接写一个样条版本，把内核从局部加权回归换成 B 样条基展开，这样做的好处是边界行为可控，拟合速度也更快：
def spline_smooth(z, y, knots=10): from scipy.interpolate import BSpline t = np.quantile(z, np.linspace(0, 1, knots + 3))[1:-1] bspl = BSpline(t, y, 3, extrapolate=False) return bspl(z)在调用处把self._smooth = _lowess_smooth换成self._smooth = spline_smooth，其余逻辑完全不用动。这就是我把_fit_one_term里所有平滑调用收敛到一个内部方法上的意义——替换实现不需要碰方向寻优代码。
6.2 接到 scikit-learn Pipeline 里做特征抽象层
投影追踪模块可以作为特征抽象层塞进数据预处理流程。比如原始特征是 50 维，PPR 学出 5 个投影方向之后，把decision_function的输出当新特征喂给下游模型：
from sklearn.pipeline import Pipeline from sklearn.preprocessing import StandardScaler pipe = Pipeline([ ("scale", StandardScaler()), ("ppr", ProjectionPursuitClassifier(max_terms=4, n_restarts=15)), ])这套接法的好处是投影方向本身携带了变量组合信息，下游模型能在更低维度上做决策。之前有个模拟项目 X 里我照搬这个模块，第一版没设早停阈值，max_terms拉到 12，训练跑了半个多小时还不见收敛。后来我强制把random_state固定、把max_terms从小往大试、再盯着残差下降曲线选拐点，才把训练时间压到分钟级。从那以后我每次碰投影追踪类模型，都会先跑一遍残差曲线再谈调参，希望帮到你。
本文还有配套的精品资源，点击获取
