MATLABとPythonでの線形計画の矛盾する解

Question

大規模な線形プログラムを実行しようとしていて、以前のコードの一部をMATLABからPythonに変換しています。しかし、問題はMATLABとPythonが劇的に矛盾する答えを与えているということです.MATLABコードは最適な解決策を見つけますが、Pythonコードでは問題は実行不可能です。これは、ell_infinity回帰またはミニマックス回帰の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方法が機能しなかったとしてWiはまた、追加された行

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

でCVXOPTを試してみました。しかし、再び、これは返さMATLAB 012と対比するdual-infeasibleのフラグ（成功）およびdual-simplexアルゴリズムを使用した。

私がテストしたデータは、関連する一連の応答変数を持つ500 x 11のデータマトリックスです。

私は出力にこのような大きな違いをもたらしたものと、それが正しいものに興味があります。混乱を招くことの1つは、入力のサイズをフルランクよりも小さくすることは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 

（この問題に対する他の製剤もありますが、いくつかは二度保存しない[link]を参照してください。）