淘宝猜价格中的数值优化：二次01规划退化成01线性规划

最新推荐文章于 2024-11-09 21:51:27 发布

_唯欣主义

最新推荐文章于 2024-11-09 21:51:27 发布

阅读量687

点赞数

分类专栏：数值优化代码练习文章标签： python

本文链接：https://blog.csdn.net/m0_50523097/article/details/128495618

版权

代码练习同时被 2 个专栏收录

9 篇文章 0 订阅

订阅专栏

数值优化

2 篇文章 0 订阅

订阅专栏

本文探讨了一个数值优化问题，即如何选择商品使得选中商品的总价格接近目标价格。通过数学建模，将问题转化为一个二次01规划，并进一步简化为线性规划。介绍了将二次项和交叉项转化为一次项的方法，最后提供了Python实现的线性规划代码，展示了系数矩阵的构造。该问题规模随着商品数量增加而增大，但数值优化算法的迭代特性使得求解仍可行。

摘要由CSDN通过智能技术生成

目标问题及运用场景

重拾数值优化学习的契机：淘宝猜价格，当我用图像和爬虫将相关商品的价格预估出来后，如何选择一件或多件商品使得被选中商品的总价格接近目标价格呢？

这是一个看似很简单的问题，我们先将其数学形式表达出来：
$\min_{x\in \{0,1\}^n}\|a^\intercal x-b \|_2^2~~~~~~b\in\mathbb{R}~~~~x,a\in\mathbb{R}^n \tag{1}$
解释一下，这里的变量 b 就是我们的目标价格，a 存储了每件商品的预估价格（我们自己预估出来的），x 是元素只能为 0 或 1 的列向量， $x_i=0$ 意味中我们选中第 i 件商品； $x_i=1$ 意味中我们不选中第 i 件商品。 $\|\cdot\|_2^2$ 用来衡量误差的经典式子（也可以尝试其他范数衡量）。

整个式子的意思是：找到一个向量 x（它说明了哪些商品选，哪些不选），使得选中的商品的总价格和目标价格差距最小。

式子化简与变形

从数值优化的角度上看，(1)式是一个约束的二次优化问题，我们首先对（1）式进行化简。
$\begin{aligned} \|a^\intercal x-b \|_2^2 &= (a^\intercal x-b)^\intercal (a^\intercal x-b)\\ &=(x^\intercal a-b)(a^\intercal x-b) \\&= x^\intercal (aa^\intercal )x+(-2ba)^\intercal x+b^2 \\ &=x^\intercal Qx+c^\intercal x+b^2 \end{aligned}$
其中 $Q=aa^\intercal, c=-2ba$ ，到这里我们可以发现式子里 $b^2$ 这一项是固定的。因此我们的（1）式等同于下面式子：
$\min_{x\in \{0,1\}^n}\|a^\intercal x-b \|_2^2\Leftrightarrow\min_{x\in \{0,1\}^n}x^\intercal Qx+c^\intercal x\tag{2}$
接下来这一步很重要，将 一个二次01规划化为一次01规划，我们先给出结论
$\min_{x\in \{0,1\}^n}x^\intercal Qx+c^\intercal x \Leftrightarrow \min_{x\in \{0,1\}^n}\sum_{1\le i<j\le n}2q_{ij}y_{ij} + \sum_{i=1}^n(c_i+q_{ii})x_i \tag{3}$
可以看到化简后（3）式没有二次项，这种变换是针对 0-1规划 设计出来的。思考一下，我们如果要把（2）式里的二次项去掉呢？具体来讲，里面的二次项还有两种形式：一种是平方项 $x_{i}^2$ ，另外一种是交叉项 $x_{i}x_j(i\ne j)$ 。0-1规划 的特殊性在平方项上展现的淋漓尽致，可以发现：
$x_{i}*x_i =x_i\tag{4}$
因此我们可以直接将平方项变成一次项！！！
而对于交叉项，就没有这么好的性质了，我们采用引入新变量的方式解决它，令
$y_{ij}=x_i*x_j~~(1\le i<j\le n)\tag{5}$
因为 x 是有约束的，所以 y 必须也有相应约束！！！
一种等价的约束是：
$\forall x_i,~x_j\in\{0,1\},~y_{ij}=x_i*x_j ~~~\mathbf{iff}~~~y_{ij}=\max\{x_i+x_j-1, 0\}\tag{6}$
读者可以自行验证，因此我们得到关于 y 的约束条件如下：
$y_{ij}\in\{0,1\}~~~~~ y_{ij}\le x_i~~~~~y_{ij}\le x_j~~~~~y_{ij}\ge x_i+x_j-1\tag{7}$
至此，我们证明了（1）到（2），（2）到（3）的等价性，由等价的 传递性，我们得到
$\min_{x\in \{0,1\}^n}\|a^\intercal x-b \|_2^2\Leftrightarrow \min_{x,y_{ij}\in \{0,1\}^n}\sum_{1\le i<j\le n}2q_{ij}y_{ij} + \sum_{i=1}^n(c_i+q_{ii})x_i \tag{8}$

问题规模

假设 n 个商品，则有 $\frac{n(n-1)}{2}$ 个 $y$ 和 $n$ 个 $x$ ，共有变量 $\frac{n(n-1)}{2} +n$ 个。相应地，0-1约束也会 $\frac{n(n-1)}{2} + n$ 个， $y$ 与 $x$ 之间的约束有 $\frac{n(n-1)}{2} * 3$ （由（7）知），共有约束 $2n^2-n$ 个。

也就是说一个形式简单的 n阶二次、带n个约束的优化最终转化为 $\frac{n(n-1)}{2} + n$ 阶、带 $2n^2-n$ 个约束的线性规划。

python实现

由于最终是用一次形式，故不需要使用库 qpsolvers，只需使用线性规划的库 pulp，另外还需使用库 numpy 进行矩阵、向量的变换。

import pulp as pu
import numpy as np

# a是商品预估价格所形成的列表，b是目标价格
def quad_01_programmer(a: list, b):
    n = len(a)
    a = np.array(a).reshape(-1, 1)
    Q = a * a.transpose()
    c = (-2) * b * a
    # print(Q)

    # row = int(3 * n * (n - 1) / 2)
    col = int(n * (n - 1) / 2 + n)

    var_xy = pu.LpVariable.dicts("z", range(col), cat="Binary")
    # print(var_xy)
    A = [2 * Q[i, j] for i in range(n - 1) for j in range(i + 1, n)]
    c = c.tolist()
    c = [int(x) for item in c for x in item]
    for i in range(n):
        c[i] += Q[i][i]
    A = A + c
    # print(A)

    pro = pu.LpProblem('', sense=pu.LpMinimize)
    # pro += pu.lpSum(A[i] * var_xy.values(i) for i in range(col))
    pro += pu.lpDot(var_xy.values(), A)

    bound = n * (n-1) >> 1
    posi  = 0
    for i in range(n - 1):
        for j in range(i + 1, n):
            pro += (var_xy[posi] - var_xy[bound + i] <= 0)
            posi += 1

    posi = 0
    for i in range(n - 1):
        for j in range(i + 1, n):
            pro += (var_xy[posi] - var_xy[bound + j] <= 0)
            posi += 1

    posi = 0
    for i in range(n - 1):
        for j in range(i + 1, n):
            pro += (var_xy[bound + i] + var_xy[bound + j] - var_xy[posi] <= 1)
            posi += 1

    pro.solve()
    solution = []
    # for i in range(n):
    #     print(pu.value(var_xy[bound + i]))
    for i in range(n):
        if pu.value(var_xy[bound + i]) == 1:
            solution.append(i)

    return solution

系数矩阵

虽然规模看起来很大，但数值优化算法一般是 迭代算法，系数矩阵的 稀疏程度 会影响到算法的速度。为此以下代码，展示了在 matlab 里面写时它的系数矩阵长什么样，并给出例子。

import pulp as pu
import numpy as np

def quad_01_looklike(a: list, b):
    n = len(a)
    a = np.array(a).reshape(-1, 1)
    Q = a * a.transpose()
    c = (-2) * b * a

    row = int(3 * n * (n - 1) / 2)
    col = int(n * (n - 1) / 2 + n)

    # print(row)
    # print(col)

    # variables = [pu.LpVariable('X%d' % i, cat='Binary') for i in range(row)]
    B = np.zeros(row).reshape(-1, 1)
    coef_object = [0 for _ in range(col)]
    posi = 0

    for i in range(n - 1):
        for j in range(i + 1, n):
            coef_object[posi] = Q[i, j] * 2
            posi += 1

    for i in range(n):
        coef_object[posi] = int(c[i])
        posi += 1

    # print(coef_object)

    A = np.zeros([row, col])
    posi = 0
    bound = n * (n - 1) >> 1
    for i in range(n - 1):
        for j in range(i + 1, n):
            A[posi, posi] = 1
            A[posi, bound + i] = -1
            posi += 1

    posi2 = 0
    for i in range(n - 1):
        for j in range(i + 1, n):
            A[posi, posi2] = 1
            A[posi, bound + j - 1] = -1
            posi += 1
            posi2 += 1

    posi2 = 0
    for i in range(n - 1):
        for j in range(i + 1, n):
            A[posi, posi2] = -1
            A[posi, bound + i] = 1
            A[posi, bound + j] = 1
            posi += 1
            posi2 += 1

    print(A)
    # prob = pu.LpProblem(sense=pu.LpMinimize)
    # prob += abs(b- sum([coef_object[i] * variables[i] for i in range(row)]))


a = [1, 2, 3, 4]
b = 6
quad_01_looklike(a, b)

得到结果如下：

[[ 1.  0.  0.  0.  0.  0. -1.  0.  0.  0.]
 [ 0.  1.  0.  0.  0.  0. -1.  0.  0.  0.]
 [ 0.  0.  1.  0.  0.  0. -1.  0.  0.  0.]
 [ 0.  0.  0.  1.  0.  0.  0. -1.  0.  0.]
 [ 0.  0.  0.  0.  1.  0.  0. -1.  0.  0.]
 [ 0.  0.  0.  0.  0.  1.  0.  0. -1.  0.]
 [ 1.  0.  0.  0.  0.  0. -1.  0.  0.  0.]
 [ 0.  1.  0.  0.  0.  0.  0. -1.  0.  0.]
 [ 0.  0.  1.  0.  0.  0.  0.  0. -1.  0.]
 [ 0.  0.  0.  1.  0.  0.  0. -1.  0.  0.]
 [ 0.  0.  0.  0.  1.  0.  0.  0. -1.  0.]
 [ 0.  0.  0.  0.  0.  1.  0.  0. -1.  0.]
 [-1.  0.  0.  0.  0.  0.  1.  1.  0.  0.]
 [ 0. -1.  0.  0.  0.  0.  1.  0.  1.  0.]
 [ 0.  0. -1.  0.  0.  0.  1.  0.  0.  1.]
 [ 0.  0.  0. -1.  0.  0.  0.  1.  1.  0.]
 [ 0.  0.  0.  0. -1.  0.  0.  1.  0.  1.]
 [ 0.  0.  0.  0.  0. -1.  0.  0.  1.  1.]]

Process finished with exit code 0