From 12668d99bd872539ff74d82a966ed049e66acde9 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Thu, 4 Apr 2024 15:27:42 +0200
Subject: Some optimizations.
---
least_squares.cpp | 57 ++++++++++++++++++++++++++++++++++---------------------
1 file changed, 35 insertions(+), 22 deletions(-)
diff --git a/least_squares.cpp b/least_squares.cpp
index 648592f..2ede87c 100644
--- a/least_squares.cpp
+++ b/least_squares.cpp
@@ -73,7 +73,7 @@ public:
return {H, gradH, hessH};
}
- Real getHamiltonian(const Vector& x) const {
+ virtual Real getHamiltonian(const Vector& x) const {
Real H;
std::tie(H, std::ignore, std::ignore) = HdHddH(x);
return H;
@@ -129,6 +129,15 @@ public:
return {H, dH, ddH};
}
+
+ /* Unfortunately benchmarking indicates that ignorning entries of a returned tuple doesn't result
+ * in those execution paths getting optimized out. It is much more efficient to compute the
+ * energy alone when only the energy is needed.
+ */
+ Real getHamiltonian(const Vector& x) const override {
+ Vector V = b + (A + 0.5 * (J * x)) * x;
+ return 0.5 * V.squaredNorm();
+ }
};
class CubicModel : public Model {
@@ -217,30 +226,36 @@ Vector findMinimum(const Model& M, const Vector& x0, Real ε = 1e-12) {
return x;
}
-Vector metropolisStep(const Model& M, const Vector& x0, Real β, Rng& r, Real ε = 1) {
- Vector Δx(M.N);
+Vector metropolisSweep(const Model& M, const Vector& x0, Real β, Rng& r, unsigned sweeps = 1, Real ε = 1) {
+ Vector x = x0;
+ Real H = M.getHamiltonian(x);
- for (Real& Δxᵢ : Δx) {
- Δxᵢ = ε * r.variate<Real, std::normal_distribution>();
- }
+ for (unsigned j = 0; j < sweeps; j++) {
+ Real rate = 0;
- Vector xNew = normalize(x0 + Δx);
+ for (unsigned i = 0; i < M.N; i++) {
+ Vector xNew = x;
- Real Hold = M.getHamiltonian(x0);
- Real Hnew = M.getHamiltonian(xNew);
+ for (Real& xNewᵢ : xNew) {
+ xNewᵢ += ε * r.variate<Real, std::normal_distribution>();
+ }
- if (exp(-β * (Hnew - Hold)) > r.uniform<Real>(0.0, 0.1)) {
- return xNew;
- } else {
- return x0;
- }
-}
+ xNew = normalize(xNew);
-Vector metropolisSweep(const Model& M, const Vector& x0, Real β, Rng& r, Real ε = 1) {
- Vector x = x0;
+ Real Hnew = M.getHamiltonian(xNew);
+
+ if (exp(-β * (Hnew - H)) > r.uniform<Real>(0.0, 0.1)) {
+ x = xNew;
+ H = Hnew;
+ rate++;
+ }
+ }
- for (unsigned i = 0; i < M.N; i++) {
- x = metropolisStep(M, x, β, r, ε);
+ if (rate / M.N < 0.5) {
+ ε /= 1.5;
+ } else {
+ ε *= 1.5;
+ }
}
return x;
@@ -300,9 +315,7 @@ int main(int argc, char* argv[]) {
for (unsigned sample = 0; sample < samples; sample++) {
QuadraticModel leastSquares(N, M, r, σ, A, J);
- for (unsigned i = 0; i < sweeps; i++) {
- x = metropolisSweep(leastSquares, x, β, r);
- }
+ x = metropolisSweep(leastSquares, x, β, r, sweeps);
x = findMinimum(leastSquares, x);
std::cout << " " << leastSquares.getHamiltonian(x) / N;
}
--
cgit v1.2.3-70-g09d2