Function minimizeComponentIca

highs/presolve/ICrashUtil.cpp:166–228

}

double minimizeComponentIca(const int col, const double mu,
                            const std::vector<double>& lambda,
                            const HighsLp& lp, double& objective,
                            std::vector<double>& residual, HighsSolution& sol) {
  // Minimize quadratic for column col.

  // Formulas for a and b when minimizing for x_j
  // a = (1/(2*mu)) * sum_i a_ij^2
  // b = -(1/(2*mu)) sum_i (2 * a_ij * (sum_{k!=j} a_ik * x_k - b_i)) + c_j (\)
  //     + sum_i a_ij * lambda_i
  // b / 2 = -(1/(2*mu)) sum_i (2 * a_ij
  double a = 0.0;
  double b = 0.0;

  for (int k = lp.a_matrix_.start_[col]; k < lp.a_matrix_.start_[col + 1];
       k++) {
    int row = lp.a_matrix_.index_[k];
    a += lp.a_matrix_.value_[k] * lp.a_matrix_.value_[k];
    // matlab but with b = b / 2
    double bracket =
        -residual[row] - lp.a_matrix_.value_[k] * sol.col_value[col];
    bracket += lambda[row];
    // clp minimizing for delta_x
    // double bracket_clp = - residual_[row];
    b += lp.a_matrix_.value_[k] * bracket;
  }

  a = (0.5 / mu) * a;
  b = (0.5 / mu) * b + 0.5 * lp.col_cost_[col];

  double theta = -b / a;
  double delta_x = 0;

  // matlab
  double new_x;
  if (theta > 0)
    new_x = std::min(theta, lp.col_upper_[col]);
  else
    new_x = std::max(theta, lp.col_lower_[col]);
  delta_x = new_x - sol.col_value[col];

  // clp minimizing for delta_x
  // if (theta > 0)
  //   delta_x = std::min(theta, lp_.col_upper_[col] - col_value_[col]);
  // else
  //   delta_x = std::max(theta, lp_.col_lower_[col] - col_value_[col]);

  sol.col_value[col] += delta_x;

  // std::cout << "col " << col << ": " << delta_x << std::endl;

  // Update objective, row_value, residual after each component update.
  objective += lp.col_cost_[col] * delta_x;
  for (int k = lp.a_matrix_.start_[col]; k < lp.a_matrix_.start_[col + 1];
       k++) {
    int row = lp.a_matrix_.index_[k];
    residual[row] -= lp.a_matrix_.value_[k] * delta_x;
    sol.row_value[row] += lp.a_matrix_.value_[k] * delta_x;