Skip to content

Commit dbc6c63

Browse files
author
Teseo Schneider
committed
Merge branch 'main' into custom-fallback
2 parents e6ada8c + b08ab23 commit dbc6c63

22 files changed

+249
-140
lines changed

‎CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -173,7 +173,7 @@ if(POLYSOLVE_LARGE_INDEX)
173173
endif()
174174

175175
target_compile_definitions(polysolve_linear PRIVATE POLYSOLVE_LINEAR_SPEC="${PROJECT_SOURCE_DIR}/linear-solver-spec.json")
176-
target_compile_definitions(polysolve PRIVATE POLYSOLVE_NON_LINEAR_SPEC="${PROJECT_SOURCE_DIR}/non-linear-solver-spec.json")
176+
target_compile_definitions(polysolve PRIVATE POLYSOLVE_NON_LINEAR_SPEC="${PROJECT_SOURCE_DIR}/nonlinear-solver-spec.json")
177177
target_compile_definitions(polysolve_linear PUBLIC POLYSOLVE_JSON_SPEC_DIR="${PROJECT_SOURCE_DIR}")
178178

179179

‎non-linear-solver-spec.json‎ renamed to ‎nonlinear-solver-spec.json‎

Lines changed: 23 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,4 @@
1-
[
2-
{
1+
[{
32
"pointer": "/",
43
"default": null,
54
"type": "object",
@@ -512,17 +511,19 @@
512511
"max_step_size_iter_final",
513512
"default_init_step_size",
514513
"step_ratio",
515-
"Armijo"
514+
"Armijo",
515+
"RobustArmijo"
516516
],
517517
"doc": "Settings for line-search in the nonlinear solver"
518518
},
519519
{
520520
"pointer": "/line_search/method",
521-
"default": "Backtracking",
521+
"default": "RobustArmijo",
522522
"type": "string",
523523
"options": [
524524
"Armijo",
525525
"ArmijoAlt",
526+
"RobustArmijo",
526527
"Backtracking",
527528
"MoreThuente",
528529
"None"
@@ -582,10 +583,27 @@
582583
},
583584
{
584585
"pointer": "/line_search/Armijo/c",
585-
"default": 0.5,
586+
"default": 1e-4,
586587
"type": "float",
588+
"min_value": 0,
587589
"doc": "Armijo c parameter."
588590
},
591+
{
592+
"pointer": "/line_search/RobustArmijo",
593+
"default": null,
594+
"type": "object",
595+
"optional": [
596+
"delta_relative_tolerance"
597+
],
598+
"doc": "Options for RobustArmijo."
599+
},
600+
{
601+
"pointer": "/line_search/RobustArmijo/delta_relative_tolerance",
602+
"default": 0.1,
603+
"type": "float",
604+
"min_value": 0,
605+
"doc": "Relative tolerance on E to switch to approximate."
606+
},
589607
{
590608
"pointer": "/box_constraints",
591609
"type": "object",

‎src/polysolve/nonlinear/PostStepData.cpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,9 +3,10 @@
33
namespace polysolve::nonlinear
44
{
55
PostStepData::PostStepData(const int iter_num,
6+
const json &solver_info,
67
const Eigen::VectorXd &x,
78
const Eigen::VectorXd &grad)
8-
: iter_num(iter_num), x(x), grad(grad)
9+
: iter_num(iter_num), solver_info(solver_info), x(x), grad(grad)
910
{
1011
}
1112
} // namespace polysolve::nonlinear

‎src/polysolve/nonlinear/PostStepData.hpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,10 +9,12 @@ namespace polysolve::nonlinear
99
{
1010
public:
1111
PostStepData(const int iter_num,
12+
const json &solver_info,
1213
const Eigen::VectorXd &x,
1314
const Eigen::VectorXd &grad);
1415

1516
const int iter_num;
17+
const json &solver_info;
1618
const Eigen::VectorXd &x;
1719
const Eigen::VectorXd &grad;
1820
};

‎src/polysolve/nonlinear/Solver.cpp‎

Lines changed: 5 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -267,12 +267,10 @@ namespace polysolve::nonlinear
267267
objFunc.solution_changed(x);
268268
}
269269

270-
objFunc.post_step(PostStepData(this->m_current.iterations, x, grad));
271-
272270
const auto g_norm_tol = this->m_stop.gradNorm;
273271
this->m_stop.gradNorm = first_grad_norm_tol;
274272

275-
StopWatch stop_watch("non-linear solver", this->total_time, m_logger);
273+
StopWatch stop_watch("nonlinear solver", this->total_time, m_logger);
276274
stop_watch.start();
277275

278276
m_logger.debug(
@@ -283,6 +281,7 @@ namespace polysolve::nonlinear
283281
this->m_stop.fDelta, this->m_stop.gradNorm, this->m_stop.xDelta);
284282

285283
update_solver_info(objFunc.value(x));
284+
objFunc.post_step(PostStepData(this->m_current.iterations, solver_info, x, grad));
286285

287286
int f_delta_step_cnt = 0;
288287
double f_delta = 0;
@@ -446,15 +445,16 @@ namespace polysolve::nonlinear
446445
// -----------
447446
const double step = (rate * delta_x).norm();
448447

448+
update_solver_info(energy);
449+
objFunc.post_step(PostStepData(this->m_current.iterations, solver_info, x, grad));
450+
449451
if (objFunc.stop(x))
450452
{
451453
this->m_status = cppoptlib::Status::UserDefined;
452454
m_error_code = ErrorCode::SUCCESS;
453455
m_logger.debug("[{}][{}] Objective decided to stop", descent_strategy_name(), m_line_search->name());
454456
}
455457

456-
objFunc.post_step(PostStepData(this->m_current.iterations, x, grad));
457-
458458
if (f_delta < this->m_stop.fDelta)
459459
f_delta_step_cnt++;
460460
else
@@ -471,8 +471,6 @@ namespace polysolve::nonlinear
471471
if (++this->m_current.iterations >= this->m_stop.iterations)
472472
this->m_status = cppoptlib::Status::IterationLimit;
473473

474-
update_solver_info(energy);
475-
476474
// reset the tolerance, since in the first iter it might be smaller
477475
this->m_stop.gradNorm = g_norm_tol;
478476
} while (objFunc.callback(this->m_current, x) && (this->m_status == cppoptlib::Status::Continue));

‎src/polysolve/nonlinear/descent_strategies/Newton.hpp‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -113,7 +113,10 @@ namespace polysolve::nonlinear
113113
const double characteristic_length,
114114
spdlog::logger &logger);
115115

116-
std::string name() const override { return internal_name() + "RegularizedNewton"; }
116+
std::string name() const override
117+
{
118+
return fmt::format("{}RegularizedNewton (reg_weight={:g})", internal_name(), reg_weight);
119+
}
117120

118121
void reset(const int ndof) override;
119122
bool handle_error() override;

‎src/polysolve/nonlinear/line_search/Armijo.cpp‎

Lines changed: 16 additions & 53 deletions
Original file line numberDiff line numberDiff line change
@@ -2,10 +2,6 @@
22

33
#include "Armijo.hpp"
44

5-
#include <polysolve/Utils.hpp>
6-
7-
#include <spdlog/spdlog.h>
8-
95
namespace polysolve::nonlinear::line_search
106
{
117
Armijo::Armijo(const json &params, spdlog::logger &logger)
@@ -14,58 +10,25 @@ namespace polysolve::nonlinear::line_search
1410
c = params["line_search"]["Armijo"]["c"];
1511
}
1612

17-
double Armijo::compute_descent_step_size(
18-
const TVector &x,
19-
const TVector &searchDir,
13+
void Armijo::init_compute_descent_step_size(
14+
const TVector &delta_x,
15+
const TVector &old_grad)
16+
{
17+
armijo_criteria = c * delta_x.dot(old_grad);
18+
assert(armijo_criteria <= 0);
19+
}
20+
21+
bool Armijo::criteria(
22+
const TVector &delta_x,
2023
Problem &objFunc,
2124
const bool use_grad_norm,
22-
const double old_energy_in,
23-
const double starting_step_size)
25+
const double old_energy,
26+
const TVector &old_grad,
27+
const TVector &new_x,
28+
const double new_energy,
29+
const double step_size) const
2430
{
25-
TVector grad(x.rows());
26-
double f_in = old_energy_in;
27-
double alpha = starting_step_size;
28-
29-
bool valid;
30-
31-
TVector x1 = x + alpha * searchDir;
32-
{
33-
POLYSOLVE_SCOPED_STOPWATCH("constraint set update in LS", constraint_set_update_time, m_logger);
34-
objFunc.solution_changed(x1);
35-
}
36-
37-
objFunc.gradient(x, grad);
38-
39-
double f = use_grad_norm ? grad.squaredNorm() : objFunc.value(x1);
40-
const double cache = c * grad.dot(searchDir);
41-
valid = objFunc.is_step_valid(x, x1);
42-
43-
while ((std::isinf(f) || std::isnan(f) || f > f_in + alpha * cache || !valid) && alpha > current_min_step_size() && cur_iter <= current_max_step_size_iter())
44-
{
45-
alpha *= step_ratio;
46-
x1 = x + alpha * searchDir;
47-
48-
{
49-
POLYSOLVE_SCOPED_STOPWATCH("constraint set update in LS", constraint_set_update_time, m_logger);
50-
objFunc.solution_changed(x1);
51-
}
52-
53-
if (use_grad_norm)
54-
{
55-
objFunc.gradient(x1, grad);
56-
f = grad.squaredNorm();
57-
}
58-
else
59-
f = objFunc.value(x1);
60-
61-
valid = objFunc.is_step_valid(x, x1);
62-
63-
m_logger.trace("ls it: {} f: {} (f_in + alpha * Cache): {} invalid: {} ", cur_iter, f, f_in + alpha * cache, !valid);
64-
65-
cur_iter++;
66-
}
67-
68-
return alpha;
31+
return new_energy <= old_energy + step_size * armijo_criteria;
6932
}
7033

7134
}; // namespace polysolve::nonlinear::line_search

‎src/polysolve/nonlinear/line_search/Armijo.hpp‎

Lines changed: 14 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,13 @@
11
#pragma once
22

3-
#include "LineSearch.hpp"
3+
#include "Backtracking.hpp"
44

55
namespace polysolve::nonlinear::line_search
66
{
7-
class Armijo : public LineSearch
7+
class Armijo : public Backtracking
88
{
99
public:
10-
using Superclass = LineSearch;
10+
using Superclass = Backtracking;
1111
using typename Superclass::Scalar;
1212
using typename Superclass::TVector;
1313

@@ -16,15 +16,21 @@ namespace polysolve::nonlinear::line_search
1616
virtual std::string name() override { return "Armijo"; }
1717

1818
protected:
19-
double compute_descent_step_size(
20-
const TVector &x,
19+
virtual void init_compute_descent_step_size(
20+
const TVector &delta_x,
21+
const TVector &old_grad) override;
22+
23+
virtual bool criteria(
2124
const TVector &delta_x,
2225
Problem &objFunc,
2326
const bool use_grad_norm,
24-
const double old_energy_in,
25-
const double starting_step_size) override;
27+
const double old_energy,
28+
const TVector &old_grad,
29+
const TVector &new_x,
30+
const double new_energy,
31+
const double step_size) const override;
2632

27-
private:
2833
double c;
34+
double armijo_criteria; ///< cached value: c * delta_x.dot(old_grad)
2935
};
3036
} // namespace polysolve::nonlinear::line_search

‎src/polysolve/nonlinear/line_search/Backtracking.cpp‎

Lines changed: 35 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -4,8 +4,6 @@
44

55
#include <spdlog/spdlog.h>
66

7-
#include <cfenv>
8-
97
namespace polysolve::nonlinear::line_search
108
{
119

@@ -20,58 +18,68 @@ namespace polysolve::nonlinear::line_search
2018
Problem &objFunc,
2119
const bool use_grad_norm,
2220
const double old_energy,
21+
const TVector &old_grad,
2322
const double starting_step_size)
2423
{
2524
double step_size = starting_step_size;
2625

27-
TVector grad(x.rows());
26+
init_compute_descent_step_size(delta_x, old_grad);
2827

29-
// Find step that reduces the energy
30-
double cur_energy = std::nan("");
31-
bool is_step_valid = false;
32-
while (step_size > current_min_step_size() && cur_iter < current_max_step_size_iter())
28+
for (; step_size > current_min_step_size() && cur_iter < current_max_step_size_iter(); step_size *= step_ratio, ++cur_iter)
3329
{
34-
TVector new_x = x + step_size * delta_x;
30+
const TVector new_x = x + step_size * delta_x;
3531

3632
try
3733
{
38-
POLYSOLVE_SCOPED_STOPWATCH("solution changed - constraint set update in LS", this->constraint_set_update_time, m_logger);
34+
POLYSOLVE_SCOPED_STOPWATCH("solution changed - constraint set update in LS", constraint_set_update_time, m_logger);
3935
objFunc.solution_changed(new_x);
4036
}
4137
catch (const std::runtime_error &e)
4238
{
4339
m_logger.warn("Failed to take step due to \"{}\", reduce step size...", e.what());
44-
45-
step_size *= step_ratio;
46-
this->cur_iter++;
4740
continue;
4841
}
4942

50-
if (use_grad_norm)
43+
if (!objFunc.is_step_valid(x, new_x))
5144
{
52-
objFunc.gradient(new_x, grad);
53-
cur_energy = grad.squaredNorm();
45+
continue;
5446
}
55-
else
56-
cur_energy = objFunc.value(new_x);
5747

58-
is_step_valid = objFunc.is_step_valid(x, new_x);
48+
const double new_energy = objFunc.value(new_x);
5949

60-
m_logger.trace("ls it: {} delta: {} invalid: {} ", this->cur_iter, (cur_energy - old_energy), !is_step_valid);
61-
62-
if (!std::isfinite(cur_energy) || cur_energy >= old_energy || !is_step_valid)
50+
if (!std::isfinite(new_energy))
6351
{
64-
step_size *= step_ratio;
65-
// max_step_size should return a collision free step
66-
// assert(objFunc.is_step_collision_free(x, new_x));
52+
continue;
6753
}
68-
else
54+
55+
m_logger.trace("ls it: {} ΔE: {}", cur_iter, new_energy - old_energy);
56+
57+
if (criteria(delta_x, objFunc, use_grad_norm, old_energy, old_grad, new_x, new_energy, step_size))
6958
{
70-
break;
59+
break; // found a good step size
7160
}
72-
this->cur_iter++;
7361
}
7462

7563
return step_size;
7664
}
65+
66+
bool Backtracking::criteria(
67+
const TVector &delta_x,
68+
Problem &objFunc,
69+
const bool use_grad_norm,
70+
const double old_energy,
71+
const TVector &old_grad,
72+
const TVector &new_x,
73+
const double new_energy,
74+
const double step_size) const
75+
{
76+
if (use_grad_norm)
77+
{
78+
TVector new_grad;
79+
objFunc.gradient(new_x, new_grad);
80+
return new_grad.norm() < old_grad.norm(); // TODO: cache old_grad.norm()
81+
}
82+
return new_energy < old_energy;
83+
}
84+
7785
} // namespace polysolve::nonlinear::line_search

0 commit comments

Comments
 (0)