Widersprüchliche Lösungen in der linearen Programmierung in MATLAB und Python

Question

Ich versuche, ein großangelegtes lineares Programm auszuführen und übersetze meinen früheren Code aus MATLAB in Python. Das Problem ist jedoch, dass MATLAB und Python mir drastisch widersprüchliche Antworten geben: Der MATLAB-Code findet die optimale Lösung, aber der Python-Code sagt, dass das Problem nicht machbar ist. Dies ist eine LP-Modellierung der ell_infinity-Regression oder der Minimax-Regression. Ich bin ziemlich zuversichtlich in der Einrichtung beider Funktionen.

MATLAB-Code:

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-Code

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 

Da die scipy Methode nicht wI arbeiten auch versucht CVXOPT mit den zusätzlichen Linien

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

Aber auch in diesem zurück eine Flagge von dual-infeasible, die den MATLAB 012 kontrastiert(Erfolg) und dual-simplex Algorithmus verwendet.

Die Daten, die ich getestet habe, sind eine 500 mal 11 Datenmatrix mit einem zugehörigen Satz von Antwortvariablen.

Ich bin interessiert, was einen so signifikanten Unterschied in den Ausgaben verursacht hat und welche davon richtig ist. Eine Sache, die Verwirrung stiftet, ist, dass das Reduzieren der Größe der Eingabe auf einen Wert, der kleiner als der vollständige Rang ist, das MATLAB-Ausgabeflag nicht zu beeinflussen scheint, aber der Python-Code keine optimale Lösung findet (wie ich denke) .

Die Daten sind bei https://www.dropbox.com/s/g4qcj1q0ncljawq/dataset.csv?dl=0 und die Matrix A ist die erste 11 Spalten und der Zielvektor ist die letzte Spalte.

Answer 1

(1) Standardmäßige variable Grenzen sind fast immer Null und unendlich. Das heißt, die meisten Löser gehen von nicht negativen Variablen aus, sofern nicht anders angegeben. Matlab verwendet einen anderen Standard. Dort sind standardmäßig Variablen frei.

Für Ihr Modell bedeutet dies, dass Sie linprog explizit mitteilen müssen, dass die Variablen x frei sind. Die Variable t kann entweder frei oder nicht negativ sein.

So

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

Das Modell für die Simplex-Methode ein wenig groß ist (die Simplex-Implementierung ist ein wenig von einem Spielzeug-Algorithmus). Die neue Interior-Point-Methode hat damit kein Problem.

(2) Die Linie

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

sollte lesen wahrscheinlich

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

(3) Beachten Sie auch, dass Ihr Modell

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

falsch ist.Es sollte sein:

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

Dies gibt als LP-Eingang:

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

(Es gibt auch andere Formulierungen für dieses Problem, einige speichern nicht zweimal im See [link].)