SciPy优化中能否将返回布尔值的函数作为非线性约束使用

现有方案的问题 #

• 你当前用布尔返回值作为非线性约束的方式不可行：SLSQP等梯度-based优化器要求约束函数输出连续可微的数值，用来计算梯度更新优化方向，布尔值没有梯度信息，优化器无法迭代。而且你当前的NonlinearConstraint缺少上下界参数，本身就无法运行。
• 正如你所说，该问题的约束本质是线性的，完全不需要引入非线性约束，用线性约束求解效率和精度都会高很多。

最优实现思路 #

核心利用两个凸集和线性变换的性质：

1. 线性变换下，凸多面体的像等于其顶点像的凸包：因此只需要保证空间A的所有顶点经过映射后落在空间B中，就能保证A中所有点映射后都在B范围内。
2. 凸多面体可以用半空间不等式表示为 C * x ≤ d，所有属于该凸多面体的点都满足该不等式组。

实现步骤 #

1. 先对空间B的所有点求凸包，将其转成C * y ≤ d的半空间表示，可用scipy.spatial.ConvexHull直接实现。
2. 将约束C * (T * a_i) ≤ d（a_i为空间A的顶点）转换为对矩阵T拉直向量的线性约束，利用克罗内克积实现矩阵运算向量化。
3. 你的映射误差目标也可以转换为线性最小二乘问题，整个问题变成带线性约束的线性最小二乘问题，用scipy.optimize.lsq_linear即可高效求解。

参考实现代码 #

import numpy as np
from scipy.spatial import ConvexHull
from scipy.optimize import lsq_linear, LinearConstraint

def get_polytope_halfspace(points):
    """将点集的凸包转换为 C @ x <= d 的半空间表示"""
    hull = ConvexHull(points)
    C = hull.equations[:, :-1]
    d = -hull.equations[:, -1]
    return C, d

def make_T_optimal(local_points, target_points):
    n_dim = local_points.shape[1]
    # 构造目标空间B的半空间约束
    B_points = np.vstack([target_points, local_points])
    C, d = get_polytope_halfspace(B_points)
    num_constraints_per_point = C.shape[0]
    num_local_vertices = local_points.shape[0]

    # 构造所有顶点的线性约束矩阵
    A_constraint = []
    b_constraint = []
    for a_i in local_points:
        # 利用克罗内克积将 T@a_i 转换为拉直向量x的线性组合
        transform_mat = np.kron(a_i.T, np.eye(n_dim))
        constraint_part = C @ transform_mat
        A_constraint.append(constraint_part)
        b_constraint.append(d)
    
    A_constraint = np.vstack(A_constraint)
    b_constraint = np.hstack(b_constraint)
    linear_constraint = LinearConstraint(A_constraint, -np.inf, b_constraint)

    # 构造最小二乘目标：最小化第一个配对点的映射误差
    a0 = local_points[0]
    b0 = target_points[0]
    A_lsq = np.kron(a0.T, np.eye(n_dim))
    b_lsq = b0

    # 求解带约束的线性最小二乘
    x0 = np.eye(n_dim).flatten()
    res = lsq_linear(A_lsq, b_lsq, constraints=linear_constraint, x0=x0)

    # 转换回矩阵形式
    T_opt = res.x.reshape(n_dim, n_dim)
    return T_opt, res

扩展说明 #

• 如果你有多个配对点需要满足映射要求，只需要把每个配对点对应的A_lsq部分和b_lsq部分堆叠即可，目标会自动变成所有配对点的误差平方和最小。
• 如果需要配对点映射完全相等，可以将T@a_i = b_i作为等式约束加入LinearConstraint中。
• 原有check_point_within_polytope函数可以用来最后验证求解得到的T是否满足要求，不需要放到优化流程中。

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