MATLAB 和 Python 中线性规划的冲突解决方案

Question

我正在尝试运行一个大型线性程序，并且正在将我以前的一些代码从 MATLAB 翻译成 Python。然而，问题是 MATLAB 和 Python 给出了截然不同的答案：MATLAB 代码找到了最优解，而 Python 代码却说问题不可行。这是 ell_infinity 回归或 minimax 回归的 LP 建模。我对这两个功能的设置相当有信心。

MATLAB 代码：

function [ x_opt, f_opt, exitflag, output] = ell_infinity_reg_solver( A, b )
% Solves the ell_infinity regression problem ||Ax - b||_inf.  That is finds
% the least t for which Ax - b < t.ones and Ax - b > -t.ones.
[n,d] = size(A) ;  
if n == 0
    f_opt  = 0 ;
    x_opt = zeros(d,1) ;
    return
end

% infinity norm
f = [ zeros(d,1); 1 ];
A = sparse([ A, -ones(n,1) ; -A, -ones(n,1) ]);
b = [ b; -b ];
[xt, f_opt, exitflag, output] = linprog(f,A,b);
x_opt = xt(1:d,:);


end

Python代码

from scipy.optimize import linprog
import numpy as np


def ell_infinity_regression_solver(A, b):
    """Compute the solution to the ell_infinity regression problem.
    x_hat = arg min_x ||Ax - b||_infty by a reformulating as LP:

    min t :
    Ax - b <= t * ones
    b - Ax >= - t*ones

    i.e. finds minimal radius ball enclosing distance between all data points
    and target values b.

    Input:
    A - array or datafram of dimension n by d
    b - target vector of size n by 1.

    Output:
    x_opt which solves the optimisation.
    f_opt the radius of the enclosing ball
    """
    n,d = A.shape

    # include a check on n and d to ensure dimensionality makes sense

    if n > 0:
        f = np.vstack((np.zeros((d,1)), [1])).ravel()
        print("shape of f {}".format(f.shape)) # should be d+1 length
        big_A = np.hstack((np.vstack((A,-A)), np.ones((2*n,1))))
        print("shape of big_A {}".format(big_A.shape)) # should be 2n*(d+1)
        #big_b = np.vstack((b,-b)) # should be length 2n
        big_b = b.append(-b) # only works for data frames -- switch to np array?
        print("Type of b {}".format(type(b)))
        print("Shape of b {}".format(b.shape))
        print("shape of big_b {}".format(big_b.shape))
        output = linprog(f, A_ub=big_A, b_ub=big_b)
        #print(output)

    else:
        f_opt = 0
        x_opt = np.zeros_like(b)


    return output

由于 scipy 方法不起作用，我还尝试 CVXOPT 添加行

c = cvxopt.matrix(f)
G = cvxopt.matrix(X)
h = cvxopt.matrix(Y)
output = cvxopt.solvers.lp(c, G, h, solver='glpk')

但是这又返回了一个 dual-infeasible 的标志，它与使用的 MATLAB exitflag=1（成功）和 dual-simplex 算法形成对比。

我测试的数据是一个 500 x 11 数据矩阵和一组相关的响应变量。

我很想知道是什么导致了输出中如此显着的差异，以及哪些是正确的。确实增加混乱的一件事是，将输入的大小减少到小于满秩的大小似乎不会影响 MATLAB 输出标志，但 Python 代码没有找到最佳解决方案（如我认为是预期的）。

数据在https://www.dropbox.com/s/g4qcj1q0ncljawq/dataset.csv?dl=0，矩阵A是前11列，目标向量是最后一列。

Answer 1

(1) 默认变量边界几乎总是零和无穷大。即，除非另有说明，否则大多数求解器都假定非负变量。 Matlab 使用不同的默认值。默认情况下变量是免费的。

对于您的模型，这意味着您需要明确告诉 linprog x 变量是免费的。 t 变量可以是自由的或非负的。

因此

output = linprog(f, A_ub=big_A, b_ub=big_b, bounds=(None,None),method='interior-point' )

模型对于单纯形法来说有点大（单纯形实现有点像玩具算法）。新的内点法用起来没问题

(2) 行

big_A = np.hstack((np.vstack((A,-A)), np.ones((2*n,1))))

应该可以阅读

big_A = np.hstack((np.vstack((A,-A)), -np.ones((2*n,1))))

(3) 另请注意您的模型

min t :
Ax - b <= t * ones
b - Ax >= - t*ones

不正确。应该是：

min t :
Ax - b <= t * ones
Ax - b >= - t*ones

这给出了 LP 输入：

min t :
 Ax - t <= b
-Ax - t  <= -b

(这个问题还有其他的表述，有的不存储A两次。参见[link])

MATLAB 和 Python 中线性规划的冲突解决方案

Conflicting solutions in linear programming in MATLAB and Python

python

matlab

optimization

linear-programming

cvxopt