Skip to content

Commit

Permalink
Now checked Hessian column deletion rigorously in unit test
Browse files Browse the repository at this point in the history
  • Loading branch information
@jajhall
jajhall committed May 25, 2024
1 parent e0676e3 commit 8afbb10
Show file tree
Hide file tree
Showing 2 changed files with 56 additions and 49 deletions.
75 changes: 39 additions & 36 deletions check/TestQpSolver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -601,13 +601,15 @@ TEST_CASE("test-semi-definite2", "[qpsolver]") {
REQUIRE(fabs(solution.col_value[1] - 2) < double_equal_tolerance);
}

void hessianProduct(const HighsHessian& hessian, const std::vector<double>& arg, std::vector<double>& result) {
void hessianProduct(const HighsHessian& hessian, const std::vector<double>& arg,
std::vector<double>& result) {
HighsInt dim = hessian.dim_;
assert(HighsInt(arg.size()) == dim);
result.resize(dim);
for (HighsInt iCol = 0; iCol < dim; iCol++) {
double sum = 0;
for (HighsInt iEl = hessian.start_[iCol]; iEl < hessian.start_[iCol+1]; iEl++)
for (HighsInt iEl = hessian.start_[iCol]; iEl < hessian.start_[iCol + 1];
iEl++)
sum += hessian.value_[iEl] * arg[hessian.index_[iEl]];
result[iCol] = sum;
}
Expand Down Expand Up @@ -677,8 +679,7 @@ TEST_CASE("test-qp-modification", "[qpsolver]") {
std::vector<double> result1;
arg0.resize(dim);
HighsRandom random;
for (HighsInt iCol = 0; iCol < dim; iCol++)
arg0[iCol] = random.fraction();
for (HighsInt iCol = 0; iCol < dim; iCol++) arg0[iCol] = random.fraction();
HighsHessian hessian0 = incumbent_model.hessian_;
arg1 = arg0;

Expand All @@ -696,17 +697,17 @@ TEST_CASE("test-qp-modification", "[qpsolver]") {

dim--;
for (HighsInt iCol = delete_col; iCol < dim; iCol++)
arg1[iCol] = arg1[iCol+1];
arg1[iCol] = arg1[iCol + 1];
arg0[delete_col] = 0;
hessianProduct(hessian0, arg0, result0);
for (HighsInt iCol = delete_col; iCol < dim; iCol++)
result0[iCol] = result0[iCol+1];
result0[iCol] = result0[iCol + 1];

arg1.resize(dim);
hessianProduct(incumbent_model.hessian_, arg1, result1);
for (HighsInt iCol = 0; iCol < dim; iCol++)
REQUIRE(result0[iCol] == result1[iCol]);

// Deleting column 0 removes only nonzero from the Hessian, so problem is an
// LP
delete_col = 0;
Expand Down Expand Up @@ -741,15 +742,14 @@ TEST_CASE("test-qp-delete-col", "[qpsolver]") {
lp.a_matrix_.value_ = {1, 1, 1, 1, 1};
hessian.dim_ = 5;
hessian.start_ = {0, 4, 7, 10, 11, 12};
hessian.index_ = {0, 1, 3, 4, 1, 2, 4, 2, 3, 4, 3, 4};
hessian.index_ = {0, 1, 3, 4, 1, 2, 4, 2, 3, 4, 3, 4};
hessian.value_ = {11, 21, 41, 51, 22, 32, 52, 33, 43, 53, 44, 55};

Highs highs;
// highs.setOptionValue("output_flag", dev_run);
highs.setOptionValue("output_flag", dev_run);
const HighsModel& incumbent_model = highs.getModel();
REQUIRE(highs.passModel(model) == HighsStatus::kOk);
// if (dev_run)
incumbent_model.hessian_.print();
if (dev_run) incumbent_model.hessian_.print();

HighsInt dim = incumbent_model.lp_.num_col_;
std::vector<double> arg0;
Expand All @@ -758,15 +758,13 @@ TEST_CASE("test-qp-delete-col", "[qpsolver]") {
std::vector<double> result1;
arg0.resize(dim);
HighsRandom random;
for (HighsInt iCol = 0; iCol < dim; iCol++)
arg0[iCol] = random.fraction();
for (HighsInt iCol = 0; iCol < dim; iCol++) arg0[iCol] = random.fraction();
HighsHessian hessian0 = incumbent_model.hessian_;
arg1 = arg0;

std::vector<HighsInt> set = {1, 3};
REQUIRE(highs.deleteCols(2, set.data()) == HighsStatus::kOk);
// if (dev_run)
incumbent_model.hessian_.print();
if (dev_run) incumbent_model.hessian_.print();

arg1[1] = arg1[2];
arg1[2] = arg1[4];
Expand All @@ -782,7 +780,8 @@ TEST_CASE("test-qp-delete-col", "[qpsolver]") {
for (HighsInt iCol = 0; iCol < dim; iCol++)
REQUIRE(result0[iCol] == result1[iCol]);

dim = 10;
dim = 100;
lp.clear();
lp.num_col_ = dim;
lp.num_row_ = 1;

Expand All @@ -799,6 +798,8 @@ TEST_CASE("test-qp-delete-col", "[qpsolver]") {
lp.a_matrix_.index_[iCol] = iCol;
lp.a_matrix_.value_[iCol] = 1;
}
lp.a_matrix_.start_.push_back(dim);

lp.sense_ = ObjSense::kMinimize;
lp.offset_ = 0;
lp.row_lower_ = {-inf};
Expand All @@ -810,62 +811,64 @@ TEST_CASE("test-qp-delete-col", "[qpsolver]") {
std::vector<double> hessian_col(dim);
for (HighsInt iCol = 0; iCol < dim; iCol++) {
hessian_col.assign(dim, 0);
HighsInt kmax = std::max(1, (dim-iCol)/3);
hessian_col[iCol] = dim * iCol + (iCol+1);
if (iCol < dim-1) {
HighsInt kmax = std::max(1, (dim - iCol) / 3);
hessian_col[iCol] = dim * iCol + iCol;
if (iCol < dim - 1) {
for (HighsInt k = 0; k < kmax; k++) {
HighsInt iRow = iCol+1 + random.integer(dim-iCol-1);
assert(iRow >= 0);
assert(iRow < dim);
hessian_col[iRow] = dim * iCol + (iRow+1);
HighsInt iRow = iCol + 1 + random.integer(dim - iCol - 1);
assert(iRow >= 0);
assert(iRow < dim);
hessian_col[iRow] = dim * iCol + iRow;
}
}
for (HighsInt iRow = iCol; iRow < dim; iRow++) {
if (hessian_col[iRow]) {
hessian.index_.push_back(iRow);
hessian.value_.push_back(hessian_col[iRow]);
hessian.index_.push_back(iRow);
hessian.value_.push_back(hessian_col[iRow]);
}
}
hessian.start_.push_back(HighsInt(hessian.index_.size()));
}

REQUIRE(highs.passModel(model) == HighsStatus::kOk);
incumbent_model.hessian_.print();
if (dev_run) incumbent_model.hessian_.print();

arg0.resize(dim);
for (HighsInt iCol = 0; iCol < dim; iCol++)
arg0[iCol] = 1;//random.fraction();
for (HighsInt iCol = 0; iCol < dim; iCol++) arg0[iCol] = random.fraction();
hessian0 = incumbent_model.hessian_;
arg1 = arg0;

std::vector<HighsInt> mask;
mask.assign(dim, 0);

HighsInt kmax = std::max(1, dim/3);
HighsInt kmax = std::max(1, dim / 3);
for (HighsInt k = 0; k < kmax; k++) {
HighsInt iRow = random.integer(dim);
assert(iRow >= 0);
assert(iRow < dim);
mask[iRow] = 1;
}
highs.deleteCols(mask.data());
if (dev_run) incumbent_model.hessian_.print();

for (HighsInt iCol = 0; iCol < dim; iCol++) {
HighsInt iRow = mask[iCol];
if (iRow < 0) {
arg0[iCol] = 0;
} else {
arg1[iRow] = arg1[iCol];
}

}
hessianProduct(hessian0, arg0, result0);
for (HighsInt iCol = 0; iCol < dim; iCol++) {
HighsInt iRow = mask[iCol];
if (iRow >= 0) result0[iRow] = result0[iCol];
}
dim = incumbent_model.hessian_.dim_;
arg1.resize(dim);
hessianProduct(incumbent_model.hessian_, arg1, result1);
for (HighsInt iCol = 0; iCol < dim; iCol++)
REQUIRE(result0[iCol] == result1[iCol]);


for (HighsInt iCol = 0; iCol < dim; iCol++) {
REQUIRE(result0[iCol] == result1[iCol]);
}
}


30 changes: 17 additions & 13 deletions src/model/HighsHessian.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,9 @@ void HighsHessian::deleteCols(const HighsIndexCollection& index_collection) {
}
for (HighsInt iCol = keep_from_col; iCol <= keep_to_col; iCol++)
new_index[iCol] = new_dim++;
// When using a mask, to_k = this->dim_, but consecutive
// keep/delete entries are accumulated, so may not need all passes
if (keep_to_col >= this->dim_ - 1) break;
}
assert(new_dim < this->dim_);
// Now perform the pass that deletes rows/columns
Expand All @@ -82,41 +85,42 @@ void HighsHessian::deleteCols(const HighsIndexCollection& index_collection) {
new_dim = 0;
HighsInt new_num_nz = 0;
HighsInt new_num_entries = 0;
std::vector<HighsInt> save_start = this->start_;
for (HighsInt k = from_k; k <= to_k; k++) {
updateOutInIndex(index_collection, delete_from_col, delete_to_col,
keep_from_col, keep_to_col, current_set_entry);
if (k == from_k) {
// Account for the initial columns being kept
for (HighsInt iCol = 0; iCol < delete_from_col; iCol++) {
assert(new_index[iCol] >= 0);
for (HighsInt iEl = this->start_[iCol]; iEl < this->start_[iCol + 1];
for (HighsInt iEl = save_start[iCol]; iEl < save_start[iCol + 1];
iEl++) {
HighsInt iRow = new_index[this->index_[iEl]];
if (iRow < 0) continue;
this->index_[new_num_entries] = iRow;
this->value_[new_num_entries] = this->value_[iEl];
if (this->value_[new_num_entries]) new_num_nz++;
new_num_entries++;
this->index_[new_num_entries] = iRow;
this->value_[new_num_entries] = this->value_[iEl];
if (this->value_[new_num_entries]) new_num_nz++;
new_num_entries++;
}
new_dim++;
this->start_[new_dim] = new_num_entries;
new_dim++;
this->start_[new_dim] = new_num_entries;
}
assert(new_dim == delete_from_col);
}
for (HighsInt iCol = keep_from_col; iCol <= keep_to_col; iCol++) {
assert(new_index[iCol] >= 0);
for (HighsInt iEl = this->start_[iCol]; iEl < this->start_[iCol + 1];
iEl++) {
for (HighsInt iEl = save_start[iCol]; iEl < save_start[iCol + 1]; iEl++) {
HighsInt iRow = new_index[this->index_[iEl]];
if (iRow < 0) continue;
this->index_[new_num_entries] = iRow;
this->value_[new_num_entries] = this->value_[iEl];
if (this->value_[new_num_entries]) new_num_nz++;
new_num_entries++;
this->index_[new_num_entries] = iRow;
this->value_[new_num_entries] = this->value_[iEl];
if (this->value_[new_num_entries]) new_num_nz++;
new_num_entries++;
}
new_dim++;
this->start_[new_dim] = new_num_entries;
}
if (keep_to_col >= this->dim_ - 1) break;
}
assert(new_dim == check_new_dim);
this->dim_ = new_dim;
Expand Down

0 comments on commit 8afbb10

Please sign in to comment.