Can one use quadratic programming to solve linear programs?

102 Views Asked by At

The quadratic programming objective function look like this:

$$\text{J}_{\text{min}} = \frac{1}{2}x^TQx + c^Tx$$

And the linear programming objective function look like this:

$$\text{J}_{\text{min}} = c^Tx$$

That means, if I say $Q = 0$, then I could use a quadratic programming method to optimize a linear programming problem.

Are there any drawbacks by doing so, rather than using Simplex method?

Because I find Hildreth's Method in C much easier than Simplex Method in C to use. It's just a faster method.

Here are two MATLAB-programs that involves linear programming and quadratic programming.

Linprog - Simpelx Method - Maximization & Minimization

% This simplex method has been written as it was C code
function [x, solution] = linprog(c, A, b, max_or_min)
  row_a = size(A, 1);
  column_a = size(A, 2);
  if(max_or_min == 0)
    % Maximization
    [x, solution] = opti(c, A, b, row_a, column_a, max_or_min);
  else
    % Minimization
    [x, solution] = opti(b, A', c, column_a, row_a, max_or_min);
  end
end
function [x, solution] = opti(c, A, b, row_a, column_a, max_or_min)
  % Clear the solution
  if(max_or_min == 0)
        x = zeros(column_a, 1);
    else
        x = zeros(row_a, 1);
  end
  % Create the tableau
  tableau = zeros(row_a + 1, column_a + row_a + 2);
  j = 1;
  for i = 1:row_a
    % First row
    tableau(i, 1:column_a) = A(i, 1:column_a);
    % Slack variable s
    j = column_a + i;
    tableau(i, j) = 1;
    % Add b vector
    tableau(i, column_a + row_a + 2) = b(i);
  end
  % Negative objective function
  tableau(row_a + 1, 1:column_a) = -c(1:column_a);
  % Slack variable for objective function
  tableau(row_a + 1, column_a + row_a + 1) = 1;
  % Do row operations
    entry = -1.0; % Need to start with a negative number because MATLAB don't have do-while! ;(
    pivotColumIndex = 0;
    pivotRowIndex = 0;
    pivot = 0.0;
    value1 = 0.0;
    value2 = 0.0;
    value3 = 0.0;
    smallest = 0.0;
    count = 0;
    solution = true;
  while(entry < 0) % Continue if we have still negative entries
    % Find our pivot column
    pivotColumIndex = 1;
    entry = 0.0;
    for i = 1:column_a + row_a + 2
      value1 = tableau(row_a + 1, i);
      if(value1 < entry)
        entry = value1;
        pivotColumIndex = i;
      end
    end
    % If the smallest entry is equal to 0 or larger than 0, break
    if(entry >= 0.0)
      break;
    end
    % If we found no solution
    if(count == 255)
        solution = false;
        return;
    end
    % Find our pivot row
    pivotRowIndex = 1;
    value1 = tableau(1, pivotColumIndex); % Value in pivot column
    value2 = tableau(1, column_a+row_a+2); % Value in the b vector
    smallest = value2/value1; % Initial smalles value1
    for i = 2:row_a
      value1 = tableau(i, pivotColumIndex); % Value in pivot column
      value2 = tableau(i, column_a+row_a+2); % Value in the b vector
      value3 = value2/value1;
      if(or(and(value3 > 0, value3 < smallest), smallest < 0))
        smallest = value3;
        pivotRowIndex = i;
      end
    end
    % We know where our pivot is. Turn the pivot into 1
    % 1/pivot * PIVOT_ROW -> PIVOT_ROW
    pivot = tableau(pivotRowIndex, pivotColumIndex); % Our pivot value
    for i = 1:column_a + row_a + 2
      value1 = tableau(pivotRowIndex, i); % Our row value at pivot row
      tableau(pivotRowIndex, i) = value1 * 1/pivot; % When value1 = pivot, then pivot will be 1
    end
    % Turn all other values in pivot column into 0. Jump over pivot row
        % -value1* PIVOT_ROW + ROW -> ROW
    for i = 1:row_a + 1
      if(i ~= pivotRowIndex)
        value1 = tableau(i, pivotColumIndex); %  This is at pivot column
        for j = 1:column_a+row_a+2
          value2 = tableau(pivotRowIndex, j); % This is at pivot row
          value3 = tableau(i, j); % This is at the row we want to be 0 at pivot column
          tableau(i, j) = -value1*value2 + value3;
        end
      end
    end
    % Count for the iteration
        count = count + 1;
  end
  % If max_or_min == 0 -> Maximization problem
  if(max_or_min == 0)
    % Now when we have shaped our tableau. Let's find the optimal solution. Sum the columns
    for i = 1:column_a
      value1 = 0; % Reset
      for j = 1:row_a + 1
        value1 = value1 + tableau(j, i); % Summary
        value2 = tableau(j, i); %  If this is 1 then we are on the selected
        % Check if we have a value that are very close to 1
        if(and(and(value1 < 1 + eps, value1 > 1 - eps), value2 > 1 - eps))
          x(i) = tableau(j, column_a+row_a+2);
        end
      end
    end
  else
    % Minimization (The Dual method) - Only take the bottom rows on the slack variables
    for i = 1:row_a
      x(i) = tableau(row_a+1, i + column_a); % We take only the bottom row at start index column_a
    end
  end
end

Quadprog - Hildreth's Method - Maximization & Minimization

function [x, solution] = quadprog(Q, c, A, b)
      % Assume that the solution is false
      solution = true;
      % Set number of iterations
      number_of_iterations = 255;
      % Same as in C code
      FLT_EPSILON = 1.19209290e-07;
      % Unconstrained solution
      x = -linsolve(Q, c);
      % Constraints difference
      K = b - A*x;
      % Check constraint violation
      if(sum(K <= 0) == 0)
        return; % No violation
      end
      % Create P
      P = linsolve(Q, A');
      % Create H = A*Q*A'
      H = A*P;
      % Solve lambda from H*lambda = -K, where lambda >= 0
      [m, n] = size(K);
      lambda = zeros(m, n);
      for km = 1:number_of_iterations
        lambda_p = lambda;
        % Use Gauss Seidel
        for i = 1:m
          w = -1.0/H(i,i)*(K(i) + H(i,:)*lambda - H(i,i)*lambda(i));
          lambda(i) = max(0, w);
        end
        % Check if the minimum convergence has been reached
        w = (lambda - lambda_p)'*(lambda - lambda_p);
        if (w < FLT_EPSILON)
          break;
        end
        % Check if the maximum iteration have been reached
        if(km == 255)
          solution = false;
          return;
        end
      end
      % Find the solution: x = -inv(Q)*c - inv(Q)*A'*lambda
      x = x - P*lambda;
    end