如何对ndarray特定轴实现凸包矢量化与插值的高效处理？

我需要处理一个形状为(2000000,19)的浮点型2D ndarray（命名为arr），以及一个形状为(19,)的浮点型1D ndarray（命名为w）。当前已实现插值凸包的处理逻辑，但因循环遍历每个样本导致速度极慢，希望找到矢量化替代方案，实现类似result = interpolated_hull(arr, w, axis=0)的无循环接口。

现有代码实现 #

依赖导入与样本数据 #

import numpy as np
from scipy.interpolate import interp1d

# Sample data
arr = np.array([[49.38639913, 49.76769437, 49.66370476, 49.49905455, 49.15242251,
        48.0518658 , 45.998071  , 45.31347273, 45.29614113, 45.25281212,
        45.0448329 , 44.61154286, 43.72763117, 42.38443203, 41.17121991,
        40.48662165, 40.35663463, 39.88001558, 39.55938095],
       [47.97387359, 47.86121818, 47.69656797, 47.18528571, 46.70000087,
        45.39146494, 43.50232035, 43.18168571, 43.82295498, 43.62364156,
        43.31167273, 42.88704848, 42.37576623, 41.0585645 , 40.37396623,
        39.09142771, 38.79679048, 38.51948485, 38.52815065]])
w = np.array([2.1017, 2.1197, 2.1374, 2.1548, 2.172 , 2.1893, 2.2068, 2.2254,
       2.2417, 2.2592, 2.2756, 2.2928, 2.3097, 2.326 , 2.3421, 2.3588,
       2.3745, 2.3903, 2.4064])

上凸包计算函数 #

def upper_andrews_hull(points: np.ndarray):
    """
    Computes the upper half of the convex hull of a set of 2D points.
    :param points: an iterable sequence of (x, y) pairs representing the points.
    """
    # 2D cross product of OA and OB vectors, i.e. z-component of their 3D cross product.
    # Returns a positive value, if OAB makes a counter-clockwise turn,
    # negative for clockwise turn, and zero if the points are collinear.
    def cross(o, a, b):
        return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])

    # Reverse the points so that we can pop from the end
    points = np.flip(points, axis=0)
    # Build upper hull
    upper = []
    for p in points:
        while len(upper) >= 2 and cross(upper[-2], upper[-1], p) <= 0:
            upper.pop()
        upper.append(p)
    # Reverse the upper hull
    upper = np.flip(np.array(upper), axis=0)

    return upper

主处理逻辑（速度瓶颈） #

result = np.empty(arr.shape)
for i in range(arr.shape[0]):
    # Create points, using w as x values, and arr as y values
    points = np.stack((w, arr[i,:]), axis=1)
    # Calculate the convex hull around the points
    hull = upper_andrews_hull(points)
    # Interpolate the hull
    interp_function = interp1d(*hull.T)
    # Store interpolation's result to have the same x references as original points
    result[i,:] = interp_function(w)

当前核心问题 #

• 循环遍历arr的200万行样本，计算效率极低
• 每个样本的凸包点数不固定，无法用统一形状的ndarray批量存储所有凸包
• 需要实现全矢量化的处理流程，彻底替换显式循环

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