使用 Scipy 查找成本矩阵与附加约束矩阵的最佳对

Question

我正在研究一个优化问题，其中有一组需要分配的对，每个对都有相关的成本。

到目前为止我一直在使用 scipy.optimize.linear_sum_assignment。

我现在想通过尝试使用附加约束矩阵来控制所选对来增加问题的复杂性。我仍然希望最小化所选对的成本，但现在有一个限制，即所选对相对于第二个矩阵的平均值保持接近某个阈值。

这是一个使用 NumPy 数组的示例：

import numpy as np

cost_matrix = np.array([[3.1, 3.0], [3.0, 3.1]])

secondary_matrix = np.array([[1.0, 10.0], [10.0, 1.0]])

# Define the threshold for the average value in secondary_matrix
average_threshold = 6.0

在这个玩具示例中，

cost_matrix

是一个方阵（很少出现这种情况），其中cost[i][j]表示配对pairs[i]和pairs[j]的成本。

secondary_matrix

是相同形状的矩阵，其中importance_values[i][j]对应于配对pairs[i]和pairs[j]的次要成本。 average_threshold 设置所选对的最大允许平均辅助成本值。

仅使用

scipy.optimize.linear

sum

assignment

会给出 [0,1] 和 [1,0] 中的对，而在控制辅助矩阵时会将对角线条目作为最佳对。这就是我想要在优化中实现的目标。

我正在寻求有关可以有效处理此类问题的合适优化方法或库的指导。到目前为止，我已经实现了这种方法：

import numpy as np
from scipy.optimize import linear_sum_assignment, minimize


def match_minimize(cost_matrix, constraint_matrix, average_threshold):
    def combined_objective(x):
        row_indices, col_indices = linear_sum_assignment(cost_matrix + x.reshape(cost_matrix.shape))
        total_cost_matrix = cost_matrix[row_indices, col_indices].sum()
        selected_pairs_average = np.mean(secondary_matrix[row_indices, col_indices])
        return total_cost_matrix + max(0, selected_pairs_average - average_threshold)

    # Initial guess for the optimization variable x
    x0 = np.zeros_like(cost_matrix).flatten()

    # Solve the optimization problem
    result = minimize(combined_objective, x0, method="COBYLA")

    # Extract the optimal assignments
    row_indices, col_indices = linear_sum_assignment(cost_matrix + result.x.reshape(cost_matrix.shape))
    return row_indices, col_indices

在这种方法中，当辅助矩阵的选定值的平均值高于平均阈值时，我会对成本添加惩罚

它似乎适用于小矩阵，但对于我的问题，我需要处理具有数千个元素的矩阵，当我有一个包含 > 255 个列元素的矩阵时，我会遇到 python 段错误。

我是优化的初学者，所以不确定我是否走在正确的道路上。

我的 match_minimize 方法看起来合理吗？有没有办法让它适用于更大的矩阵？或者我应该使用不同的方法？

非常感谢所有建议和建议。谢谢

Answer 1

非常简单的 LP 构造；尽管您需要一个权重参数来决定平均误差相对于主要分配成本的重要性。尝试将下面的

mean_err_weight

更改为 0.1 到 100 之间，观察差异。

import numpy as np
from scipy.optimize import milp, LinearConstraint
import scipy.sparse


def mean_assignment(
    asn_cost: np.ndarray,
    mean_cost: np.ndarray,
    target_mean: float,
    mean_err_weight: float = 1,
) -> tuple[
    np.ndarray,  # assignments
    float,  # mean
    float,  # error from mean
]:
    """
    Perform linear sum assignment with a secondary objective of minimizing deviation from a mean;
    in other words,
    minimize sum(asn_cost[assign]) + mean_err_weight*(sum(mean_cost[assign])/n - target_mean)

    :param asn_cost: Assigment costs, 2D
    :param mean_cost: Mean costs; must be same shape as assignment costs
    :param target_mean: Number to which the mean of mean_cost[assign] will be pulled
    :param mean_err_weight: Optimization objective weight for mean cost
    :return: Assignment array of same shape as asn_cost; resulting mean of assigned mean_cost; and
        resulting error from mean
    """

    if asn_cost.shape != mean_cost.shape:
        raise ValueError('Shape mismatch')

    # Ensure a long matrix during optimization
    transpose = asn_cost.shape[0] < asn_cost.shape[1]
    if transpose:
        asn_cost = asn_cost.T
        mean_cost = mean_cost.T

    m, n = asn_cost.shape

    # Variables:
    #    m*n assignments, binary
    #    mean deviance, continuous

    cost_coeff = np.concatenate((
        asn_cost.ravel(),    # Assignment costs
        (mean_err_weight,),  # Cost of error from mean
    ))
    integrality = np.ones(m*n + 1)  # Assignments are binary
    integrality[-1] = 0             # Mean error is continuous

    # There must be at most one assignment per row (Kronecker)
    row_excl = LinearConstraint(
        A=scipy.sparse.hstack(
            (
                scipy.sparse.kron(scipy.sparse.eye(m), np.ones((1, n))),
                scipy.sparse.csc_array((m, 1)),
            ),
            format='csc',
        ),
        lb=np.zeros(m),
        ub=np.ones(m),
    )

    # There must be exactly one assignment per column
    col_excl = LinearConstraint(
        A=scipy.sparse.hstack(
            (scipy.sparse.eye(n),)*m +
            (scipy.sparse.csc_array((n, 1)),),
            format='csc',
        ),
        lb=np.ones(n),
        ub=np.ones(n),
    )

    # The error from mean is taken as the absolute
    # err >= +sum(secondary[assign])/sum(assign) - target
    # err >= -sum(secondary[assign])/sum(assign) + target
    # cost_matrix*assign + ... - n*err <= n*target
    # cost_matrix*assign + ... + n*err >= n*target
    mean_err_abs = LinearConstraint(
        A=scipy.sparse.bmat([
            [mean_cost.ravel(), [-n]],
            [mean_cost.ravel(), [+n]],
        ], format='csc'),
        lb=[-np.inf, n*target_mean],
        ub=[n*target_mean, +np.inf],
    )

    result = milp(
        c=cost_coeff, integrality=integrality,
        constraints=(row_excl, col_excl, mean_err_abs),
    )
    if not result.success:
        raise ValueError(result.message)
    assignments = result.x[:-1].reshape((m, n)).astype(int)
    mean_err = result.x[-1]
    mean = mean_cost[assignments.astype(bool)].mean()
    if transpose:
        assignments = assignments.T
    return assignments, mean, mean_err


def test() -> None:
    assign, mean, mean_err = mean_assignment(
        asn_cost=np.array([
            [5, 2],
            [3, 4],
            [6, 7],
        ]),
        mean_cost=np.array([
            [1., 10.],
            [10., 1.],
            [2.,  2.],
        ]),
        target_mean=6.,
        mean_err_weight=1,  # change to 0.5 to see cost tradeoff
    )

    print(assign)
    print('mean =', mean)
    print('mean err =', mean_err)


if __name__ == '__main__':
    test()