CGAL users discussion list

[Jean-Philippe Bauchet <>]

Hello,

I am encountering a problem with the module QP_Solver of CGAL 4.9.

I am currently trying to optimize a quadratic function f(X) = 1/2 X^T D X + C^T X under linear constraints AX <= B, but the exact representation mode of each variable x_i of X, makes the program return incorrect results when the number of variables in X (say N) increases, whereas other solvers like IPOPT still return accurate results.

Indeed, due to the use of an exact number type (CGAL::Gpmzf, CGAL::MP_Float) each variable x seems to be computed as x = p / q where p and q are floats : thus there exists m and e such that q = m * 10^e. The problem is, e decreases with the size of X. Beyond a certain threshold for this size, p and q get null and I finally obtain x = 0/0. More details can be found in the attached PDF file. A short runnable code illustrating the problem is also enclosed with this e-mail.

Is there any way to fix this problem ?

Thank you very much for your time, your help and your consideration.

Best regards,

Jean-Philippe Bauchet

Attachment: mu.dat
Description: Binary data

#include <iostream>
#include <string>
#include <cassert>
#include <utility>
#include <vector>
#include <CGAL/basic.h>
#include <CGAL/QP_models.h>
#include <CGAL/QP_functions.h>
#include <CGAL/Gmpzf.h>

#pragma warning(disable:4996)

typedef CGAL::Gmpzf ET;
typedef CGAL::Quadratic_program_from_iterators <double**, double*, CGAL::Const_oneset_iterator<CGAL::Comparison_result>, bool*, double*, bool*, double*, double**, double*> Program;
typedef CGAL::Quadratic_program_solution<ET> Solution;


bool allocate_data_matrix(const char* filename, double ***p_A, int *p_rows, int *p_cols) {
	FILE* A_file = fopen(filename, "r");
	if (A_file != NULL) {
		fscanf(A_file, "%i %i", p_rows, p_cols);
		int rows = *p_rows;
		int cols = *p_cols;
		*p_A = new double*[rows];
		double** A = *p_A;
		for (int i = 0 ; i < rows ; i++) {
			A[i] = new double[cols];
			for (int j = 0 ; j < cols ; j++) {
				double a_ij;
				fscanf(A_file, "%lf", &a_ij);
				A[i][j] = a_ij;
			}
		}
		fclose(A_file);
		return true;
	} else {
		return false;
	}
}


void desallocate_data_matrix(double **A, int rows) {
	for (int i = 0; i < rows; i++) delete[] A[i];
	delete[] A;
}


void allocate_bounds(bool** p_fl, bool** p_fu, double** p_l, double** p_u, int individuals, int variables, double M) {
	*p_fl = new bool[variables];
	*p_fu = new bool[variables];
	*p_l = new double[variables];
	*p_u = new double[variables];
	bool *fl = *p_fl, *fu = *p_fu;
	double *l = *p_l, *u = *p_u;
	for (int i = 0; i < individuals; i++) {
		fl[i] = true;
		fu[i] = true;
		l[i] = -M;
		u[i] = M;
	}
	for (int i = individuals; i < variables; i++) {
		fl[i] = false;
		fu[i] = false;
		l[i] = 0;
		u[i] = 0;
	}
}


void allocate_objective_function(double*** p_D, double** p_C, int individuals, int variables, double d_coef) {
	*p_D = new double*[variables];
	double** D = *p_D;
	for (int i = 0; i < variables; i++) {
		D[i] = new double[variables];
		if (i < individuals) {
			for (int j = 0; j < variables; j++) D[i][j] = (i == j ? d_coef : 0);
		} else {
			for (int j = 0; j < variables; j++) D[i][j] = 0;
		}
	}
	*p_C = new double[variables];
	double* C = *p_C;
	for (int i = 0; i < individuals; i++) C[i] = 0;
	for (int i = individuals; i < variables; i++) C[i] = 1;
}


void allocate_constraints(double*** p_A, double **p_B, int individuals, std::vector< std::pair<int, int> > & potentials, int variables, double** mu, double** targets) {
	*p_A = new double*[variables];
	double** A = *p_A;
	for (int i = 0; i < variables; i++) {
		A[i] = new double[2 * potentials.size()];
		for (int j = 0; j < 2 * potentials.size(); j++) A[i][j] = 0;
	}
	for (int p = 0; p < potentials.size(); p++) {
		int i = potentials[p].first;
		int j = potentials[p].second;
		A[i][2 * p + 0] = -2 * mu[i][j];
		A[j][2 * p + 0] = 2 * mu[i][j];
		A[individuals + p][2 * p + 0] = -1;
		A[i][2 * p + 1] = 2 * mu[i][j];
		A[j][2 * p + 1] = -2 * mu[i][j];
		A[individuals + p][2 * p + 1] = -1;
	}
	*p_B = new double[2 * potentials.size()];
	double* B = *p_B;
	for (int p = 0; p < potentials.size(); p++) {
		int i = potentials[p].first;
		int j = potentials[p].second;
		B[2 * p + 0] = -2 * mu[i][j] * targets[i][j];
		B[2 * p + 1] = 2 * mu[i][j] * targets[i][j];
	}
}


void desallocate_constraints(double** A, double* B, int variables) {
	for (int i = 0; i < variables; i++) {
		delete[] A[i];
	}
	delete[] A;
	delete[] B;
}


void desallocate_bounds(bool* fl, bool* fu, double* l, double* u) {
	delete[] fl;
	delete[] fu;
	delete[] l;
	delete[] u;
}


void desallocate_objective_function(double** D, double* C, int variables) {
	for (int i = 0; i < variables; i++) delete[] D[i];
	delete[] D;
	delete[] C;
}


int main()
{
	// Constant parameters, shouldn't be modified
	const double K = 0.9; // Scale factor
	const double M = 7;   // Range domain for the variables

	// Reads input data from files
	double **mu, **targets;
	int mu_rows, mu_cols, targets_rows, targets_cols;
	bool read_mu = allocate_data_matrix("mu.dat", &mu, &mu_rows, &mu_cols);
	bool read_targets = allocate_data_matrix("targets.dat", &targets, &targets_rows, &targets_cols);
	if (!read_mu || !read_targets) {
		std::cout << "Didn't read at least one input file : mu.dat, targets.dat" << std::endl;
		return 1;
	}

	// Checks
	assert(mu_rows == mu_cols);
	assert(mu_rows == targets_rows);
	assert(targets_rows == targets_cols);

	// Kept individuals : maximum value for this variable is 16 (mu_rows = 16)
	int individuals = 12;
	assert(individuals <= mu_rows);

	// Finds interactions
	std::vector< std::pair<int, int> > potentials;
	for (int i = 0 ; i < individuals ; i++) {
		for (int j = i + 1 ; j < individuals; j++) {
			if (mu[i][j] > 0) {
				potentials.push_back(std::make_pair(i, j));
			}
		}
	}
	int variables = individuals + potentials.size();
	std::cout << "Individuals : " << individuals << std::endl;
	std::cout << "Potentials  : " << potentials.size() << std::endl;
	std::cout << "Variables   : " << variables << std::endl << std::endl;

	// Defines the matrices describing our quadratic program
	bool *fl, *fu;
	double *l, *u;
	allocate_bounds(&fl, &fu, &l, &u, individuals, variables, M);

	double **D, *C;
	double d_coef = 2 * (1 - K) / (M * M);
	allocate_objective_function(&D, &C, individuals, variables, d_coef);

	double **A, *B;
	allocate_constraints(&A, &B, individuals, potentials, variables, mu, targets);

	CGAL::Const_oneset_iterator<CGAL::Comparison_result> r(CGAL::SMALLER);

	// Program options
	CGAL::Quadratic_program_options options;
	options.set_pricing_strategy(CGAL::QP_FILTERED_DANTZIG);
	options.set_auto_validation(true);

	// Solves
	Program qp(variables, 2 * int(potentials.size()), A, B, r, fl, l, fu, u, D, C);
	Solution S = CGAL::solve_quadratic_program(qp, ET(), options);
	
	std::cout << S << std::endl;

	// Frees memory
	desallocate_bounds(fl, fu, l, u);
	desallocate_objective_function(D, C, variables);
	desallocate_constraints(A, B, variables);
	desallocate_data_matrix(mu, mu_rows);
	desallocate_data_matrix(targets, targets_rows);

	return 0;
}

Attachment: targets.dat
Description: Binary data

Attachment: CGAL_QP.pdf
Description: Adobe PDF document