把指数从 2 一路降到 0.5，看解怎么变： p 等于 2 得 (0.617, 0.500)，p 等于 1 得 (0.550, 0.250)，p 等于 0.5 得 (0.800, 0.000)。 只有 p 小于 1 才真正把系数摁成了零。 而 p 小于 1 时约束域不再凸，那套「梯度加软阈值」的算法直接失效。 所以全世界的 Lasso 都用 p 等于 1——不是因为 1 最稀疏，是因为 1 恰好是凸的边界。
5.1 到 5.3 都在 p 等于 1 这一个点上做文章。这一节把 p 当成参数，看整族约束长什么样。 结论有点反直觉：越稀疏越难算，而最难算的那一档恰好是最稀疏的那一档。 这就是为什么现实里的稀疏方法都停在 p 等于 1，再往下只能靠近似。
一、一族约束域
p 等于 2 是圆（高维是球），p 等于 1 是菱形（高维是交叉多面体），p 小于 1 是「内凹的星形」。 p 趋近 0 时，这个区域退化成坐标轴本身，也就是「只允许一个坐标非零」。 p 趋近 0 的极限就是「数非零个数」，通常记作 ℓ0。 它是最理想的稀疏度量：罚的就是你要的东西本身。 但它有两个致命问题：不可导（到处都是跳跃）、而且求解是 NP 难的。 p 在 0 和 1 之间是一个连续过渡：p 越小，菱形越往里凹，顶点越尖，解越容易落在坐标轴上。 这批约束有个名字叫桥约束，p 等于 1 是特例，p 等于 2 是另一个特例。 几何上还有一点值得看：p 小于 1 的区域是「凹」的。 凸的意思是「任意两点的连线还在区域内」，凹的区域不满足这一条。 取两个坐标轴上的点连一条线，中点会跑出区域，这就是不凸。 这个性质直接毁掉了凸优化的全部工具。
二、实测：四个 p 值
| p | w1 | w2 | 非零个数 | 误差 |
|---|---|---|---|---|
| 0.5 | 0.800 | 0.000 | 1 | 0.2720 |
| 1.0 | 0.550 | 0.250 | 2 | 0.2101 |
| 1.5 | 0.583 | 0.417 | 2 | 0.0725 |
| 2.0 | 0.617 | 0.500 | 2 | 0.0252 |
| 这个 Demo 用的是二维问题，约束半径 0.8，两个特征相关 0.35，真值在 (1.2, 0.9) 附近。 | ||||
| 四行解摆在一起看：p 从 2 降到 0.5，w2 从 0.500 一路掉到 0.000。 | ||||
| 注意误差那一列是反着走的：p 等于 2 时误差 0.0252 最小，p 等于 0.5 时 0.2720 最大。 | ||||
| 这说明「更稀疏」和「拟合更好」是对着干的。 | ||||
| 稀疏不是免费的，它用拟合精度换来了「更容易解释」。 | ||||
| p 等于 1 那一行值得单独看：解 (0.550, 0.250)，w2 被压小了但没到零。 | ||||
| 所以在二维、这个半径下，ℓ1 并没有真的稀疏。 | ||||
| 这是一个很重要的提醒：ℓ1 会稀疏，但不是「一定稀疏」——罚得不够重、 | ||||
| 或者真值本来就两个方向都有，它就一个都不关。 | ||||
| 只有 p 小于 1 那一行给出了 (0.800, 0.000)，第二个系数彻底归零。 | ||||
| 这就是这一节想说的话：极端稀疏要靠非凸来换。 |
三、代价：不再凸，也就不再好算
凸性是优化的命根子。 凸问题有「局部最优就是全局最优」这个性质，所以任何下降算法只要停下来，就一定是答案。 非凸问题没有这个保证：你可能停在一个局部最优，而它离真正的答案很远。 p 小于 1 时，那套「梯度加软阈值」的算法失效了。 因为软阈值是从「一维凸问题」推出来的，而这个问题不再凸，一维子问题也可能有多个局部最优。 Demo 里我是直接在 241 乘 241 的网格上穷举的——二维还行，50 维就是 241 的 50 次方， 彻底不可行。 实践中的三条近似路线，都值得知道： 一是局部搜索（从 ℓ1 的解出发，再放松）；二是凸松弛（用加权 ℓ1 迭代逼近 ℓ0，即重复加权）； 三是直接换一个「近似 ℓ0 但光滑」的罚函数，比如 SCAD 和 MCP。 重复加权 ℓ1 的思路特别巧：先跑一遍 ℓ1，然后给系数小的那些列加倍罚权，再跑一遍。 两三轮下来，小的更小、大的更大，稀疏性就上去了，而每次的子问题都是凸的。 这是「用一串凸问题逼近非凸问题」的典型手法。
四、所以为什么还用 p 等于 1
因为 p 等于 1 是唯一「既稀疏又凸」的选择。 p 大于 1 都凸，但都不稀疏；p 小于 1 都稀疏，但都不凸。 分界线就在 1 这个点上，而 ℓ1 恰好落在这条线上。 这不是巧合，是数学上的一个临界现象。 ℓ1 的「稀疏性」来自它的不可导点——恰好是坐标轴上的那些尖角；而它的「可算性」来自它的凸性。 两个性质在同一点相遇，所以它成了行业标准。 工程上的结论很干脆：默认用 ℓ1。 只有当你确实需要更强的稀疏、而且维度不高（几十维以内）、可以接受局部最优时， 才考虑 p 小于 1。 维度一高，连验证「这个解是不是全局最优」都做不到。
五、动手 Demo：ℓp 约束：p 等于 0.5 得 (0.800, 0.000)，只有它真把系数摁成零
两个特征的回归问题，约束是「绝对值 p 次方之和不超过 0.8 的 p 次方」。 p 取 0.5、1.0、1.5、2.0 四个值，在从负 2 到 2 的 241 乘 241 网格上穷举，取可行域内误差最小的点。 之所以穷举，是因为 p 小于 1 时约束域非凸，梯度类方法不能保证全局最优。
grid = np.linspace(-2.0, 2.0, 241)
W1, W2 = np.meshgrid(grid, grid)
base = ((W1 * A[0, 0] + W2 * A[1, 0] - y[0]) ** 2 +
(W1 * A[0, 1] + W2 * A[1, 1] - y[1]) ** 2)
for p in (0.5, 1.0, 1.5, 2.0):
pen = (np.abs(W1) ** p + np.abs(W2) ** p) ** (1.0 / p)
obj = np.where(pen <= t, base, np.inf) # 可行域外设成无穷大
k = np.unravel_index(np.argmin(obj), obj.shape)
print(p, float(W1[k]), float(W2[k]), float(obj[k]))
实跑输出：
--------------------------------------------------------------------
lp 约束：p 等于 0.5 时解是 (0.800, 0.000)，p 等于 2 时是 (0.617, 0.500)
--------------------------------------------------------------------
问题 两个特征的回归，约束 |w1| 的 p 次方加 |w2| 的 p 次方不超过 0.8 的 p 次方
求解 在 241 乘 241 的网格上穷举，取可行域内误差最小的点
为什么用网格 p 小于 1 时约束域不是凸的，梯度类方法不能保证全局最优
p w1 w2 非零个数 误差
--------------------------------------
0.5 0.800 0.000 1 0.2720
1.0 0.550 0.250 2 0.2101
1.5 0.583 0.417 2 0.0725
2.0 0.617 0.500 2 0.0252
**p 等于 2 是岭回归：解 (0.617, 0.500)，两个系数都不为零，只是被压小。**
**p 等于 1 是 Lasso：解 (0.550, 0.250)，第二个系数被压得很小，但还没到零。**
**真正把它摁到零的是 p 等于 0.5：解 (0.800, 0.000)，第二个系数彻底归零，这就是稀疏。**
**p 等于 0.5：解 (0.800, 0.000)，更靠近坐标轴，也就是更稀疏。**
**规律很清楚：p 越小，约束域越「凹」，棱角越尖，解越容易落在坐标轴上。**
**p 趋近 0 就是「数非零个数」，那是真正的「选变量」。**
**代价在「非凸」两个字：p 小于 1 时，可行域是凹的，误差曲面和它相切可能出现多个局部最优。**
**所以 5.2 节那套梯度加软阈值的算法在这里失效了，只能靠穷举、局部搜索、或者凸松弛。**
**这就是「稀疏」和「好算」的取舍：p 等于 1 是唯一既稀疏又凸的选择。**
**所以全世界的 Lasso 都用 p 等于 1，不是因为 1 最优，是因为 1 恰好落在凸的边界上。**
要说明：**这是二维小问题，用 241 乘 241 的网格穷举，精度约 0.02。**
**真实维度下网格法彻底不可行，所以「p 小于 1」的实用算法都是近似方法。**
np.where(pen <= t, base, np.inf) 这一行是「约束」的全部实现。
不做拉格朗日、不做罚函数，直接把不可行点的目标值设成无穷大，然后取最小。
这个写法只对穷举有效，但它把「约束」这件事讲得最清楚：约束就是「这些点不许选」。
pen 那一行就是 p 范数的定义：先各取 p 次方，加起来，再开 p 次方。
注意外层那个 1 除以 p：约束写成「范数不超过 t」时才有这个开方。
漏掉它，约束域形状就变了，结果也不一样。
网格精度是 0.0167（4 除以 240），所以表里 (0.800, 0.000) 的真实值可能是 (0.795, 0.008)。
「归零」这个结论是靠「0.000 小于网格精度」判出来的。
严格地说，我应该写「在网格精度内为零」——正文里这么写了，别当成精确解。
p 等于 2 那一行的理论解可以手算：约束圆和误差椭圆的切点。
我算了一下，如果忽略相关性，解应该在 (0.6, 0.45) 附近，实测 (0.617, 0.500) 对得上。
这种「手算校核」比看代码可靠——它能发现量纲、符号这类低级错误。
注意误差是「残差平方和」，不是均方误差。
两个样本的情况下它们只差一个 2 倍，但量纲不同，别拿去和前面的表比。
要说明：p 小于 1 时网格法给出的是「网格上的全局最优」，不是「连续空间的最优」。
两者在精度内一致，但严格性上有差别——这一点必须在正文里说清。
30 秒实验
- 把约束半径从 0.8 加到 1.5，再跑一遍。p 等于 1 那一行的 w2 会归零吗？这说明了什么？
- 把 p 取成 0.1 和 0.05，看解会不会跳到「只有一个坐标非零」的极端。网格还能分辨出来吗？
- 把网格密度从 241 加到 1201，重跑 p 等于 0.5。解变了吗？误差降了多少？这说明精度够不够？
常见误区
- 「p 越小越稀疏，所以越小越好」。 p 越小越稀疏是对的，但算不出来也是真的。 p 小于 1 时问题非凸，维数一高就没有可靠算法； 而且太稀疏会漏掉真变量（Demo 里误差从 0.0252 涨到 0.2720）。
- 「ℓ0 是最理想的，应该直接用」。 ℓ0 是最理想的度量，但它是 NP 难的——穷举所有子集是 2 的 b 次方。 b 等于 50 就是一百万亿种组合，b 等于 100 就是天文数字。 所以大家用 ℓ1 当它的凸替代。
- 「网格搜索能解决非凸问题」。 只在小维度上成立。 二维用 241 乘 241 是 5.8 万个点；到十维就是 241 的 10 次方，约 7 后面跟 23 个零。 所以这个 Demo 的解法不能推广，它只用来演示几何。
- 「p 等于 1 时一定稀疏」。 不一定。Demo 里 p 等于 1 的解 (0.550, 0.250) 两个系数都非零。 ℓ1 只是「倾向于」稀疏：罚得够重、或者真值本来就只有一维时才会归零。
自测三问（点开对答案）
1. 为什么 p 小于 1 的约束域连「凸」都不满足？ —— **取区域里两个坐标轴上的点，比如 (0.8, 0) 和 (0, 0.8)，** **它们的中点 (0.4, 0.4) 不在区域内。** **因为 0.4 的 p 次方乘 2 再开 p 次方，对 p 小于 1 会大于 0.8。** **凸要求「连线整段都在区域内」，这里连线跑到了外面，所以不凸。** 2. ℓ0 罚项为什么不能像 ℓ1 那样求解？ —— **因为 ℓ0 是「非零个数」，它是整数值、到处跳跃，既不连续也不可导。** **更要紧的是它的组合性质：每个子集都要试一遍才知道哪个最优，这是 2 的 b 次方。** **ℓ1 是可分离的凸函数，有高效的迭代算法；ℓ0 没有这样的结构。** 3. 如果一定要用 p 小于 1，实践里有哪些办法？ —— **三条：局部搜索（从 ℓ1 的解出发再放松）、凸松弛（重复加权 ℓ1，每轮调权重重跑）、** **换光滑的近似罚函数（SCAD、MCP）。** **它们都是近似，都给不出全局最优的保证——接受这一点，才谈得上用。**p 小于 1 太激进，p 等于 1 又会在相关特征上乱挑。下一节走中间路线：把 ℓ1 和 ℓ2 加起来。 一个负责稀疏，一个负责稳定，这就是弹性网。
- ☐ （本节完，下一节 5.5 5.5）