带和约束的面积差快速优化方法咨询（含Scipy实现）

一、先优化损失函数的计算效率 #

当前损失函数的计算已用到numpy向量化，但可以通过预计算固定项进一步减少重复运算：

1. 预计算梯形权重：原代码的梯形求和逻辑可以转化为加权点积，提前计算好每个点的权重，避免每次迭代都做切片和加法操作。
2. 优化矩阵乘法形式：将数组转置后用矩阵乘法替代逐元素相乘再求和，代码更清晰且numpy的矩阵乘法效率更高。

优化后的损失函数代码：

import numpy as np

# 预计算所有固定项
xpoint_diffs = xpoints[1:] - xpoints[:-1]
# 计算梯形法对应的每个点权重
weights = np.zeros(len(at))
weights[0] = xpoint_diffs[0] / 2
weights[-1] = xpoint_diffs[-1] / 2
weights[1:-1] = (xpoint_diffs[:-1] + xpoint_diffs[1:]) / 2
# 转置数组，方便后续矩阵乘法（shape: 500×4）
A = np.array([a1, a2, a3, a4]).T
target = at

def lossfunc(x):
    # 计算线性组合与目标的绝对差，再加权求和
    abs_diff = np.abs(target - A @ x)
    return abs_diff @ weights

二、换用线性规划求解器（更高效的核心方案） #

你的问题属于带约束的L1最小化问题，变量仅4个，完全可以转化为线性规划（LP）问题求解，效率远高于scipy.minimize的通用非线性优化方法。

转化思路 #

目标是最小化sum(w_i * |t_i - (c₁a₁ᵢ + c₂a₂ᵢ + c₃a₃ᵢ + c₄a₄ᵢ)|)（w_i为预计算的梯形权重）。引入辅助变量z_i≥0，将绝对误差约束转化为两组线性不等式：

• t_i - (Aᵢ@c) ≤ z_i
• -(t_i - (Aᵢ@c)) ≤ z_i
  此时目标函数变为sum(w_i * z_i)，同时保留原约束：

1. c₁ + c₂ + c₃ + c₄ = 1
2. 0 ≤ cⱼ ≤ 1（j=1~4）

具体实现代码 #

from scipy.optimize import linprog

# 沿用之前预计算的A、weights、target

# 构建线性规划目标系数：前4位是c₁~c₄（系数为0），后500位是z₁~z₅₀₀（系数为weights）
c = np.concatenate([np.zeros(4), weights])

# 构建不等式约束矩阵与右端项
A_ub = np.zeros((2*500, 4+500))
b_ub = np.zeros(2*500)

for i in range(500):
    # 第一组约束：t_i - Aᵢ@c ≤ z_i → -Aᵢ@c - z_i ≤ -t_i
    A_ub[i, :4] = -A[i]
    A_ub[i, 4+i] = -1
    b_ub[i] = -target[i]
    # 第二组约束：Aᵢ@c - t_i ≤ z_i → Aᵢ@c - z_i ≤ t_i
    A_ub[500+i, :4] = A[i]
    A_ub[500+i, 4+i] = -1
    b_ub[500+i] = target[i]

# 构建等式约束：c₁+c₂+c₃+c₄=1
A_eq = np.zeros((1, 4+500))
A_eq[0, :4] = 1
b_eq = np.array([1.0])

# 变量边界：cⱼ∈[0,1]，z_i≥0
bounds = [(0,1)]*4 + [(0, None)]*500

# 调用线性规划求解器（highs是高效的现代LP求解器）
result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method='highs')

# 提取最优系数与最小损失
optimal_c = result.x[:4]
min_loss = result.fun

效率优势 #

线性规划求解器（如linprog的highs方法）是针对凸优化问题的专用算法，收敛速度远快于通用非线性优化器；且无需计算损失函数的梯度（避免数值梯度的额外计算开销），即使约束规模大，求解效率依然很高。

三、其他小优化 #

• 若精度允许，将数组类型从float64改为float32，减少内存占用并加快运算速度。
• 所有固定计算（如权重、转置数组）都提前完成，避免在优化迭代中重复计算。

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