Method SolveL1Regression

src/colmap/estimators/rotation_averaging_impl.cc:655–717

}

bool RotationAveragingSolver::SolveL1Regression(
    RotationAveragingProblem& problem) {
  LeastAbsoluteDeviationSolver::Options l1_solver_options;
  l1_solver_options.max_num_iterations = 10;
  l1_solver_options.solver_type =
      LeastAbsoluteDeviationSolver::Options::SolverType::SupernodalCholmodLLT;
  l1_solver_options.ridge_regularization = options_.ridge_regularization;

  LeastAbsoluteDeviationSolver l1_solver(l1_solver_options,
                                         problem.ConstraintMatrix());
  if (!l1_solver.Valid()) {
    LOG(ERROR) << "L1 regression linear solver factorization failed";
    return false;
  }

  double prev_norm = 0;
  double curr_norm = 0;

  Eigen::VectorXd step(problem.NumParameters());

  int iteration = 0;

  for (; iteration < options_.max_num_l1_iterations; iteration++) {
    VLOG(2) << "L1 ADMM iteration: " << iteration;

    problem.ComputeResiduals();

    prev_norm = curr_norm;

    step.setZero();
    if (!l1_solver.Solve(problem.Residuals(), &step) ||
        step.array().isNaN().any()) {
      LOG(ERROR) << "L1 regression solve failed (iteration " << iteration
                 << ")";
      return false;
    }

    if (VLOG_IS_ON(2))
      LOG(INFO) << "residual:"
                << (problem.ConstraintMatrix() * step - problem.Residuals())
                       .array()
                       .abs()
                       .sum();

    curr_norm = step.norm();
    problem.UpdateState(step);

    // Check convergence.
    constexpr double kEps = 1e-12;
    if (problem.AverageStepSize(step) <
            options_.l1_step_convergence_threshold ||
        std::abs(prev_norm - curr_norm) < kEps) {
      if (std::abs(prev_norm - curr_norm) < kEps)
        LOG(INFO) << "std::abs(prev_norm - curr_norm) < " << kEps;
      iteration++;
      break;
    }
    l1_solver_options.max_num_iterations =