Solving alternative systems using `Highs::basisSolve` and `Highs::basisTransposeSolve`

[chrhansk]

I an wondering about the documentation specifically of the functions

• Highs_getBasisSolve and
• Highs_getBasisTransposeSolve.

As far as I see, the type of LPs HiGHs expects to solve are given as

$$ \begin{aligned} \min_{x} & c^{T} x \ \text{s.t. } & L \leq Ax \leq U\ & l \leq x \leq u \end{aligned} $$

where $A$ is num_rows times num_cols, $L, U$ num_rows, and $l, u$ num_cols (for num_cols = Highs_getNumCol() and num_rows = Highs_getNumRow())

As far as I see it, a (row) basis of this problem should consist of num_cols of the inequalities and bounds, forming a regular matrix, which, when solved against their respective right hand sides should reproduce the associated basic solution. In particular, if $A$ is empty, the row basis should consist only of variable bounds, and the corresponding matrix should be the identity. This matrix should (?) coincide with the basis matrix B mentioned in the documentation.

I would like to use the functions above in order to solve B against a right-hand side, but I am confused regarding the dimensions in the documentation:

• The dimension of the right-hand side input vector b is not specified (I assumed it should be num_rows as well?)
• The output should contain num_rows entries.

I would have thought that input / output have a size equal to the size of the row basis (num_cols). I tried to understand the operations using an example:

#include "interfaces/highs_c_api.h"

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

int main() {
  // This example does exactly the same as the minimal example above,
  // but illustrates the full C API.  It first forms and solves the LP
  //
  // Min    f  =  x_0 +  x_1 + 3
  // s.t.                x_1 <= 7
  //        5 <=  x_0 + 2x_1 <= 15
  //        6 <= 3x_0 + 2x_1
  // 0 <= x_0 <= 4; 1 <= x_1
  //
  // It then solves it as a maximization, then as a MIP.
  //
  const int num_col = 2;
  const int num_row = 3;
  const int num_nz = 5;
  // Define the optimization sense and objective offset
  int sense = kHighsObjSenseMinimize;
  const double offset = 3;

  double* rhs = (double*) malloc(sizeof(double)*num_row);
  double* solution_vector = (double*) malloc(sizeof(double)*num_row);
  HighsInt* solution_index = (HighsInt*) malloc(sizeof(HighsInt)*num_row);
  HighsInt solution_num_nz = 0;

  // Define the column costs, lower bounds and upper bounds
  const double col_cost[2] = {1.0, 1.0};
  const double col_lower[2] = {0.0, 1.0};
  const double col_upper[2] = {4.0, 1.0e30};
  // Define the row lower bounds and upper bounds
  const double row_lower[3] = {-1.0e30, 5.0, 6.0};
  const double row_upper[3] = {7.0, 15.0, 1.0e30};
  // Define the constraint matrix column-wise
  const int a_format = kHighsMatrixFormatColwise;
  const int a_start[2] = {0, 2};
  const int a_index[5] = {1, 2, 0, 1, 2};
  const double a_value[5] = {1.0, 3.0, 1.0, 2.0, 2.0};

  int run_status;
  int model_status;
  double objective_function_value;
  int simplex_iteration_count;
  int64_t mip_node_count;
  int primal_solution_status;
  int dual_solution_status;
  int basis_validity;

  // Create a Highs instance
  void* highs = Highs_create();

  // Pass the LP to HiGHS
  run_status = Highs_passLp(highs, num_col, num_row, num_nz, a_format, sense, offset,
                            col_cost, col_lower, col_upper,
                            row_lower, row_upper,
                            a_start, a_index, a_value);
  assert(run_status == kHighsStatusOk);

  // Solve the incumbent model
  run_status = Highs_run(highs);
  // The run must be successful
  assert(run_status == kHighsStatusOk);
  // Get the model status - which must be optimal
  model_status = Highs_getModelStatus(highs);
  assert(model_status == kHighsModelStatusOptimal);

  printf("Run status = %d; Model status = %d\n", run_status, model_status);

  // Get scalar information about the solution
  Highs_getDoubleInfoValue(highs, "objective_function_value", &objective_function_value);
  Highs_getIntInfoValue(highs, "simplex_iteration_count", &simplex_iteration_count);
  Highs_getIntInfoValue(highs, "primal_solution_status", &primal_solution_status);
  Highs_getIntInfoValue(highs, "dual_solution_status", &dual_solution_status);
  Highs_getIntInfoValue(highs, "basis_validity", &basis_validity);

  // The primal and dual solution status values should indicate feasibility
  assert(primal_solution_status == kHighsSolutionStatusFeasible);
  assert(dual_solution_status == kHighsSolutionStatusFeasible);
  // The basis status should indicate validity
  assert(basis_validity == kHighsBasisValidityValid);

  double* col_value = (double*)malloc(sizeof(double) * num_col);
  double* col_dual = (double*)malloc(sizeof(double) * num_col);
  double* row_value = (double*)malloc(sizeof(double) * num_row);
  double* row_dual = (double*)malloc(sizeof(double) * num_row);

  int* col_basis_status = (int*)malloc(sizeof(int) * num_col);
  int* row_basis_status = (int*)malloc(sizeof(int) * num_row);

  // Get the primal and dual solution
  Highs_getSolution(highs, col_value, col_dual, row_value, row_dual);
  // Get the basis
  Highs_getBasis(highs, col_basis_status, row_basis_status);

  // Report the column primal and dual values, and basis status
  for (int i = 0; i < num_col; i++) {
    printf("Col%d = %lf; dual = %lf; status = %d; \n", i, col_value[i], col_dual[i], col_basis_status[i]);
  }
  // Report the row primal and dual values, and basis status
  for (int i = 0; i < num_row; i++) {
    printf("Row%d = %lf; dual = %lf; status = %d; \n", i, row_value[i], row_dual[i], row_basis_status[i]);
  }
  printf("Objective value = %g; Iteration count = %d\n", objective_function_value, simplex_iteration_count);

  for(int i = 0; i < num_row;++i)
  {
    rhs[i] = 10*(i + 1);
    solution_vector[i] = -1.;
    solution_index[i] = -1;
    printf("Rhs %d = %f\n", i, rhs[i]);
  }

  run_status = Highs_getBasisSolve(highs, rhs,
                                   solution_vector,
                                   &solution_num_nz,
                                   solution_index);

  assert(run_status == kHighsStatusOk);

  for(int j = 0; j < solution_num_nz; ++j)
  {
    printf("Solution %d = %f\n", solution_index[j], solution_vector[j]);
  }

  return 0;
}

The row basis consists of constraints 2 & 3 at their lower bounds, producing the solution (.5, 2.25). So I would have expected that solving the basis would solve against the corresponding matrix

$$ \begin{pmatrix} 1 & 2 \ 3 & 2 \end{pmatrix}, $$

which would have inputs and outputs in R^2. This does not seem to be the case however as the input and output is in R^3.

Could you explain to me the definition of the basis matrix and the semantics of the two functions?

[jajhall]

There's only ever one "basis" in the HiGHS simplex solver, and it is of dimension num_rows. This is true of any simplex solver except the weird stuff in Soplex and (hence) CPLEX. Hence I didn't feel that it was necessary to specify the dimension of the RHS and solution in these calls.

Given an LP in the format above, HiGHS introduces (implicit) slack variables so that it has an equation system $Ax+y=0$, where $-U\le y\le -L$. A basis is given by a set of num_rows (basic) variables whose values are solved using the equations once the remaining num_cols (nonbasic) variables have been fixed - normally at bounds. Hence the only known "invertible representation" is of the matrix of dimension num_rows formed by the columns of $[A~ I]$ corresponding to the basic variables.

In the example, what you refer to as the "row basis" corresponds to the nonbasic slacks. It has a "kernel" role in Soplex, but isn't identified specifically in HiGHS.

When num_rows >= num_cols, it is the matrix $B_0$ in the basis matrix

$$\left[ \begin{matrix} B_0&\ B_1&I \end{matrix}\right] $$

Hence, with HiGHS, it is still possible to solve systems involving this matrix $B_0$ by augmenting the RHS of dimension num_cols with zero values - whether the regular or transposed system is to be solved.

When num_rows < num_cols, the HiGHS basis matrix takes the role of $B_0$ in the above, and the systems involving your matrix can be solved, since the non-trivial part of the solve involves the HiGHS basis matrix.