寻找给定函数下方最大凸函数的高效Python实现方案

我有长度均为n的浮点型参数向量X和值向量Y，需要构造向量Yc，使得(X, Yc)构成逐点满足Yc ≤ Y的最大凸函数图像。朴素算法时间复杂度为O(n²)，希望找到更高效的实现方式。从数学角度看这属于凸包问题，但现有Python库大多聚焦2D/3D凸包，不清楚最快获取Yc的方法。

已尝试用scipy.spatial.ConvexHull计算凸包，但不知道如何提取下边界对应的Yc。

要得到满足要求的Yc，本质是计算Y关于X的下凸包（lower convex hull），再通过线性插值得到所有X点对应的Yc。下凸包的时间复杂度可以做到O(n log n)，远优于朴素算法。

核心思路 #

1. 确保X有序（若原X无序，先排序X和对应Y）；
2. 用Andrew算法计算下凸包，得到构成最大凸函数的顶点；
3. 对所有X点进行线性插值，得到每个X对应的Yc，保证Yc ≤ Y且是最大凸函数值。

完整代码 #

import numpy as np
import matplotlib.pyplot as plt

# 原始数据
X = np.linspace(1, 5, 9)
Y = np.array([2, 3, 1, 4, 6, 2, 1, 4, 3])

# 步骤1：确保X升序（若原X无序，先排序）
sorted_indices = np.argsort(X)
X_sorted = X[sorted_indices]
Y_sorted = Y[sorted_indices]
points = np.column_stack((X_sorted, Y_sorted))

# 步骤2：Andrew算法计算下凸包（O(n log n)复杂度）
def lower_convex_hull(points):
    # 按X坐标排序，X相同则按Y排序
    points_sorted = points[np.lexsort((points[:, 1], points[:, 0]))]
    lower = []
    for p in points_sorted:
        # 维护下凸包：移除不满足凸性的点
        while len(lower) >= 2:
            a, b = lower[-2], lower[-1]
            # 计算叉积判断是否为左转（保证凸性）
            cross = (b[0] - a[0]) * (p[1] - b[1]) - (b[1] - a[1]) * (p[0] - b[0])
            if cross <= 0:
                lower.pop()
            else:
                break
        lower.append(p)
    return np.array(lower)

lower_hull = lower_convex_hull(points)
lower_X, lower_Y = lower_hull[:, 0], lower_hull[:, 1]

# 步骤3：线性插值得到所有X对应的Yc
Yc_sorted = np.interp(X_sorted, lower_X, lower_Y)

# 恢复到原X的顺序（若之前排序过）
Yc = np.zeros_like(Y)
Yc[sorted_indices] = Yc_sorted

# 可视化结果
plt.figure(figsize=(10, 6))
plt.plot(X, Y, 'o', label='原始数据点')
plt.plot(X, Yc, 'r-', linewidth=2, label='最大凸函数Yc')
plt.plot(lower_X, lower_Y, 'g--', marker='s', label='下凸包顶点')
plt.legend()
plt.xlabel('X')
plt.ylabel('Y')
plt.title('满足Yc ≤ Y的最大凸函数')
plt.show()

print("原始Y数组:", Y)
print("构造的Yc数组:", Yc)

说明 #

• Andrew算法直接计算下凸包，避免了筛选完整凸包的额外步骤，效率更高；
• np.interp保证每个X点的Yc是下凸包分段线性函数的对应值，满足Yc ≤ Y且是最大的凸函数；
• 若原X无序，排序后处理再恢复顺序即可，不影响最终结果。

内容的提问来源于stack exchange，提问作者SBF