Work on Stokes lines, and new method involving parametric fits.

4 files changed, 284 insertions, 249 deletions

diff --git a/dynamics.hpp b/dynamics.hpp
index 10b1be2..561714e 100644
--- a/dynamics.hpp
+++ b/dynamics.hpp
@@ -113,3 +113,20 @@ std::tuple<Real, Vector<Scalar>> metropolis(const Tensor<Scalar, p>& J, const Ve
 
   return {E, z};
 }
+
+template <class Real, class Scalar, int p, class Distribution, class Generator>
+Vector<Scalar> randomSaddle(const Tensor<Scalar, p>& J, Distribution d, Generator& r, Real ε) {
+  Vector<Scalar> zSaddle;
+  bool foundSaddle = false;
+
+  while (!foundSaddle) {
+    Vector<Scalar> z0 = normalize(randomVector<Scalar>(J.dimension(0), d, r.engine()));
+
+    try {
+      zSaddle = findSaddle(J, z0, ε);
+      foundSaddle = true;
+    } catch (std::exception& e) {}
+  }
+
+  return zSaddle;
+}
diff --git a/langevin.cpp b/langevin.cpp
index dc7c5cf..fe26332 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -1,8 +1,10 @@
 #include <getopt.h>
+#include <limits>
 #include <stereographic.hpp>
 #include <unordered_map>
 #include <list>
 
+#include "Eigen/src/Eigenvalues/ComplexEigenSolver.h"
 #include "complex_normal.hpp"
 #include "p-spin.hpp"
 #include "dynamics.hpp"
@@ -84,11 +86,10 @@ int main(int argc, char* argv[]) {
   Real T = 1;    // temperature
   Real Rκ = 0;   // real part of distribution parameter
   Real Iκ = 0;   // imaginary part of distribution parameter
-
   // simulation parameters
   Real ε = 1e-15;
   Real εJ = 1e-5;
-  Real δ = 1e-2;   // threshold for determining saddle
+  Real δ = 1;   // threshold for determining saddle
   Real Δ = 1e-3;
   Real γ = 1e-1;   // step size
   unsigned t = 1000; // number of Langevin steps
@@ -140,7 +141,7 @@ int main(int argc, char* argv[]) {
 
   complex_normal_distribution<Real> d(0, 1, 0);
 
-  ComplexTensor J = generateCouplings<Complex, PSPIN_P>(N, complex_normal_distribution<Real>(0, σ, κ), r.engine());
+  ComplexTensor J = generateCouplings<Complex, p>(N, complex_normal_distribution<Real>(0, σ, κ), r.engine());
 
   std::function<Real(const ComplexTensor&, const ComplexVector&)> energyNormGrad = []
     (const ComplexTensor& J, const ComplexVector& z) {
@@ -149,106 +150,58 @@ int main(int argc, char* argv[]) {
       return W;
     };
 
-
-  std::function<Real(const ComplexTensor&, const ComplexVector&)> energyThres = [N, κ]
-    (const ComplexTensor& J, const ComplexVector& z) {
-      Real a = z.squaredNorm() / (Real)N;
-      /*
-      Complex ε = getHamiltonian(J, z) / (Real)N;
-      Real t = norm(ε) / getThresholdEnergyDensity(p, κ, ε, a);
-      */
-
-      Real amin = 1.1547;
-      Real amax = 1.46806;
-
-      /*
-      Real tmin = 1;
-      Real tmax = 1.17514;
-      */
-
-      Real E = 0;
-
-      if (a > amax) {
-        E += pow(a - amax, 2);
-      } else if (a < amin) {
-        E += pow(a - amin, 2);
-      }
-
-      /*
-      if (t > tmax) {
-        E += pow(t - tmax, 2);
-      } else if (t < tmin) {
-        E += pow(t - tmin, 2);
-      }
-      */
-
-      return E;
-    };
-
-  Real largestThreshold = 0;
-  ComplexVector saddlePastThreshold;
-  bool foundSaddle = false;
-
-#pragma omp parallel default(none) shared(largestThreshold, saddlePastThreshold, foundSaddle, M, J, energyThres, d, N, γ, ε, κ)
-  {
-    Rng rPrivate;
-    ComplexVector z = normalize(randomVector<Complex>(N, d, rPrivate.engine()));
-
-    while (largestThreshold < 1) { // Until we find a saddle past the threshold...
-      std::tie(std::ignore, z) = metropolis(J, z, energyThres, (Real)0.1, γ, M, d, rPrivate.engine());
-      try {
-        ComplexVector latestSaddle = findSaddle(J, z, ε);
-        Real threshold = getProportionOfThreshold(κ, J, latestSaddle);
-
-#pragma omp critical
-        {
-          if (threshold > largestThreshold + 1e-6) {
-            largestThreshold = threshold;
-            saddlePastThreshold = latestSaddle;
-            std::cerr << "Found saddle with threshold porportion " << largestThreshold << std::endl;;
-          }
-        }
-      } catch (std::exception& e) {}
-    }
-  }
-
-  std::cerr << "Found saddle with energy beyond threshold." << std::endl;
+  ComplexVector zSaddle = randomSaddle(J, d, r, ε);
+  std::cerr << "Found saddle." << std::endl;
 
   ComplexVector zSaddleNext;
-  Real closestSaddle = std::numeric_limits<Real>::infinity();
-  ComplexVector z = saddlePastThreshold;
-
-  while (closestSaddle > δ) { // Until we find two saddles sufficiently close...
-    std::tie(std::ignore, z) = metropolis(J, saddlePastThreshold, energyNormGrad, T, γ, 1e4, d, r.engine());
+  bool foundSaddle = false;
+  while (!foundSaddle) {
+    ComplexVector z0 = normalize(zSaddle + δ * randomVector<Complex>(N, d, r.engine()));
     try {
-      zSaddleNext = findSaddle(J, z, ε);
-      Real saddleDistance = (zSaddleNext - saddlePastThreshold).norm();
-      if (saddleDistance < closestSaddle && saddleDistance > 1e-2) {
-        closestSaddle = saddleDistance;
-        std::cerr << "Nearby saddle: found saddle at distance " << saddleDistance << std::endl;
+      zSaddleNext = findSaddle(J, z0, ε);
+      Real saddleDistance = (zSaddleNext - zSaddle).norm();
+      if (saddleDistance / N > 1e-2) {
+        foundSaddle = true;
       }
     } catch (std::exception& e) {}
   }
 
-  std::cerr << "Found sufficiently nearby saddles, perturbing J to equalize Im H." << std::endl;
+  auto [H1, dH1, ddH1] = hamGradHess(J, zSaddle);
+  auto [H2, dH2, ddH2] = hamGradHess(J, zSaddleNext);
 
-  ComplexVector z1 = saddlePastThreshold;
-  ComplexVector z2 = zSaddleNext;
+  Eigen::ComplexEigenSolver<ComplexMatrix> ces;
+  ces.compute(ddH1);
 
-  std::tie(J, z1, z2) = matchImaginaryEnergies(J, z1, z2, 1e-14, ε, r);
+  Real φ = atan2( H2.imag() - H1.imag(), H1.real() - H2.real());
+  std::cerr << (zSaddle - zSaddleNext).norm() / (Real)N << " " << φ << " " << H1 * exp(Complex(0, φ)) << " " << H2 * exp(Complex(0, φ)) << std::endl;
+  Real smallestNorm = std::numeric_limits<Real>::infinity();
+  for (Complex z : ces.eigenvalues()) {
+    if (norm(z) < smallestNorm) {
+      smallestNorm = norm(z);
+    }
+  }
+  std::cerr << smallestNorm << std::endl;
 
-  std::cerr << "Im H is now sufficently close, starting to relax rope." << std::endl;
+  J = exp(Complex(0, φ)) * J;
 
-  std::cerr << "Threshold proportions of saddles are " << getProportionOfThreshold(κ, J, z1) << " and " << getProportionOfThreshold(κ, J, z2) << std::endl;
+  /*
+  if (stokesLineTest<Real>(J, zSaddle, zSaddleNext, 10, 4)) {
+    std::cerr << "Found a Stokes line" << std::endl;
+//    stokesLineTestNew<Real>(J, zSaddle, zSaddleNext, 10, 3);
+  } else {
+    std::cerr << "Didn't find a Stokes line" << std::endl;
+  }
+  */
 
-  Rope stokes(10, z1, z2, J);
 
-  for (unsigned i = 0; i < 9; i++) {
-    stokes.relaxDiscreteGradient(J, 1e6, 0.1, 0);
+  Cord test(J, zSaddle, zSaddleNext, 5);
+  test.relax(J, 20, 1, 1e5);
 
-    std::cout << stokes.n() << " " << stokes.cost(J) << " " << stokes.length() << std::endl;
+  std::cout << test.z0.transpose() << std::endl;
+  std::cout << test.z1.transpose() << std::endl;
 
-    stokes = stokes.interpolate();
+  for (Vector<Complex>& g : test.gs) {
+    std::cout << g.transpose() << std::endl;
   }
 
   return 0;
diff --git a/p-spin.hpp b/p-spin.hpp
index b2bdf07..e5cd195 100644
--- a/p-spin.hpp
+++ b/p-spin.hpp
@@ -89,3 +89,18 @@ std::tuple<Real, Vector<Scalar>> WdW(const Tensor<Scalar, p>& J, const Vector<Sc
 
   return {W, dW};
 }
+
+template <class Scalar>
+Matrix<Scalar> dzDot(const Vector<Scalar>& z, const Vector<Scalar>& dH) {
+  Real z² = z.squaredNorm();
+  return (dH.conjugate() - (dH.dot(z) / z²) * z.conjugate()) * z.adjoint() / z²;
+}
+
+template <class Scalar>
+Matrix<Scalar> dzDotConjugate(const Vector<Scalar>& z, const Vector<Scalar>& dH, const Matrix<Scalar>& ddH) {
+  Real z² = z.squaredNorm();
+  return -ddH + (ddH * z.conjugate()) * z.transpose() / z²
+          + (z.dot(dH) / z²) * (
+              Matrix<Scalar>::Identity(ddH.rows(), ddH.cols()) - z.conjugate() * z.transpose() / z²
+            );
+}
diff --git a/stokes.hpp b/stokes.hpp
index cbe5437..b2694a9 100644
--- a/stokes.hpp
+++ b/stokes.hpp
@@ -4,6 +4,8 @@
 #include "complex_normal.hpp"
 #include "dynamics.hpp"
 
+#include <iostream>
+
 class ropeRelaxationStallException: public std::exception {
   virtual const char* what() const throw() {
     return "Gradient descent stalled.";
@@ -11,49 +13,6 @@ class ropeRelaxationStallException: public std::exception {
 };
 
 template <class Scalar>
-Vector<Scalar> variation(const Vector<Scalar>& z, const Vector<Scalar>& z´, const Vector<Scalar>& z´´, const Vector<Scalar>& dH, const Matrix<Scalar>& ddH) {
-  Real z² = z.squaredNorm();
-  Real z´² = z´.squaredNorm();
-
-  Vector<Scalar> ż = zDot(z, dH);
-  Real ż² = ż.squaredNorm();
-
-  Real Reż·z´ = real(ż.dot(z´));
-
-  Matrix<Scalar> dż = (dH.conjugate() - (dH.dot(z) / z²) * z.conjugate()) * z.adjoint() / z²;
-  Matrix<Scalar> dżc = -ddH + (ddH * z.conjugate()) * z.transpose() / z²
-    + (z.dot(dH) / z²) * (
-        Matrix<Scalar>::Identity(ddH.rows(), ddH.cols()) - z.conjugate() * z.transpose() / z²
-      );
-
-  Vector<Scalar> dLdz = - (
-        dżc * z´ + dż * z´.conjugate() - (dż * ż.conjugate() + dżc * ż) * Reż·z´ / ż²
-      ) / sqrt(ż² * z´²) / 2;
-
-  Vector<Scalar> ż´ = -(ddH * z´).conjugate() + ((ddH * z´).dot(z) / z²) * z.conjugate() + (
-        dH.dot(z) * z´.conjugate() + dH.dot(z´) * z.conjugate() - (
-          dH.dot(z) * (z´.dot(z) + z.dot(z´)) / z²
-        ) * z.conjugate()
-      ) / z²;
-
-  Real dReż·z´ = real(ż.dot(z´´) + ż´.dot(z´));
-
-  Vector<Scalar> ddtdLdz´ = - (
-      (
-        ż´.conjugate() - (
-          Reż·z´ * z´´.conjugate() + dReż·z´ * z´.conjugate()
-          - (Reż·z´ / z´²) * (z´´.dot(z´) + z´.dot(z´´)) * z´.conjugate() 
-          ) / z´²
-        )
-      - 0.5 * (
-        (ż.dot(ż´) + ż´.dot(ż)) / ż² + (z´´.dot(z´) + z´.dot(z´´)) / z´²
-        ) * (ż.conjugate() - Reż·z´ / z´² * z´.conjugate())
-      ) / sqrt(ż² * z´²) / 2;
-
-  return dLdz - ddtdLdz´;
-}
-
-template <class Scalar>
 class Rope {
   public:
     std::vector<Vector<Scalar>> z;
@@ -93,65 +52,10 @@ class Rope {
       return l;
     }
 
-    template <int p>
-    Real error(const Tensor<Scalar, p>& J) const {
-      Scalar H0, HN;
-      std::tie(H0, std::ignore, std::ignore) = hamGradHess(J, z.front());
-      std::tie(HN, std::ignore, std::ignore) = hamGradHess(J, z.back());
-
-      Real ImH = imag((H0 + HN) / 2.0);
-
-      Real err = 0;
-
-      for (unsigned i = 1; i < z.size() - 1; i++) {
-        Scalar Hi;
-        std::tie(Hi, std::ignore, std::ignore) = hamGradHess(J, z[i]);
-
-        err += pow(imag(Hi) - ImH, 2);
-      }
-
-      return sqrt(err);
-    }
-
     Vector<Scalar> dz(unsigned i) const {
       return z[i + 1] - z[i - 1];
     }
 
-    Vector<Scalar> ddz(unsigned i) const {
-      return 4.0 * (z[i + 1] + z[i - 1] - 2.0 * z[i]);
-    }
-
-    template <int p>
-    std::vector<Vector<Scalar>> generateGradientδz(const Tensor<Scalar, p>& J) const {
-      std::vector<Vector<Scalar>> δz(z.size());
-
-#pragma omp parallel for
-      for (unsigned i = 1; i < z.size() - 1; i++) {
-        Vector<Scalar> dH;
-        Matrix<Scalar> ddH;
-        std::tie(std::ignore, dH, ddH) = hamGradHess(J, z[i]);
-
-        δz[i] = variation(z[i], dz(i), ddz(i), dH, ddH);
-      }
-
-      for (unsigned i = 1; i < z.size() - 1; i++) {
-        δz[i] = δz[i].conjugate() - (δz[i].dot(z[i]) / z[i].squaredNorm()) * z[i].conjugate();
-//        δz[i] = δz[i] - ((δz[i].conjugate().dot(dz(i))) / dz(i).squaredNorm()) * dz(i).conjugate();
-      }
-
-      // We return a δz with average norm of one.
-      Real mag = 0;
-      for (unsigned i = 1; i < z.size() - 1; i++) {
-        mag += δz[i].norm();
-      }
-
-      for (unsigned i = 1; i < z.size() - 1; i++) {
-        δz[i] /= mag / n();
-      }
-
-      return δz;
-    }
-
     template <int p>
     std::vector<Vector<Scalar>> generateDiscreteGradientδz(const Tensor<Scalar, p>& J, Real γ) const {
       std::vector<Vector<Scalar>> δz(z.size());
@@ -185,29 +89,38 @@ class Rope {
           dC +=  0.5 * (ż[i + 1].conjugate() - dz(i + 1).conjugate() * real(ż[i + 1].dot(dz(i + 1))) / dz(i + 1).squaredNorm()) / (ż[i + 1].norm() * dz(i + 1).norm());
         }
 
-        dC += - γ * (z[i - 1] + z[i + 1]).conjugate();
+        dC += γ * (2 * z[i] - z[i - 1] - z[i + 1]).conjugate();
 
-        δz[i] = dC;
-      }
+        δz[i] = dC.conjugate();
 
-      Real size = 0;
-      for (unsigned i = 1; i < z.size() - 1; i++) {
-        δz[i] = δz[i].conjugate() - (δz[i].dot(z[i]) / z[i].squaredNorm()) * z[i].conjugate();
+        δz[i] -=  z[i].conjugate() *  z[i].conjugate().dot(δz[i]) /  z²;
       }
 
       return δz;
     }
 
-    template<class Gen>
-    std::vector<Vector<Scalar>> generateRandomδz(Gen& r) const {
-      std::vector<Vector<Scalar>> δz(z.size());
+    void spread() {
+      Real l = length();
+
+      Real a = 0;
+      unsigned pos = 0;
+
+      std::vector<Vector<Scalar>> zNew = z;
 
-      complex_normal_distribution<> d(0, 1, 0);
       for (unsigned i = 1; i < z.size() - 1; i++) {
-        δz[i] = randomVector<Scalar>(z[0].size(), d, r);
+        Real b = i * l / (z.size() - 1);
+
+        while (b > a) {
+          pos++;
+          a += (z[pos] - z[pos - 1]).norm();
+        }
+
+        Vector<Scalar> δz = z[pos] - z[pos - 1];
+
+        zNew[i] = normalize(z[pos] - (a - b) / δz.norm() * δz);
       }
 
-      return δz;
+      z = zNew;
     }
 
     template<int p>
@@ -222,8 +135,6 @@ class Rope {
           rNew.z[i] = normalize(z[i] - (δ * Δl) * δz[i]);
         }
 
-        rNew.spread();
-
         if (rNew.cost(J, γ) < cost(J, γ)) {
           break;
         } else {
@@ -235,71 +146,33 @@ class Rope {
         }
       }
 
+//      rNew.spread();
+
       z = rNew.z;
 
       return δ;
     }
 
-    void spread() {
-      Real l = length();
-
-      Real a = 0;
-      unsigned pos = 0;
-
-      std::vector<Vector<Scalar>> zNew = z;
-
-      for (unsigned i = 1; i < z.size() - 1; i++) {
-        Real b = i * l / (z.size() - 1);
-
-        while (b > a) {
-          pos++;
-          a += (z[pos] - z[pos - 1]).norm();
-        } 
-
-        Vector<Scalar> δz = z[pos] - z[pos - 1];
-
-        zNew[i] = normalize(z[pos] - (a - b) / δz.norm() * δz);
-      }
-
-      z = zNew;
-    }
-
-    template <int p>
-    void relaxGradient(const Tensor<Scalar, p>& J, unsigned N, Real δ0) {
-      Real δ = δ0;
-      try {
-        for (unsigned i = 0; i < N; i++) {
-          std::vector<Vector<Scalar>> δz = generateGradientδz(J);
-          δ = 1.1 * perturb(J, δ, δz);
-        }
-      } catch (std::exception& e) {
-      }
-    }
-
     template <int p>
     void relaxDiscreteGradient(const Tensor<Scalar, p>& J, unsigned N, Real δ0, Real γ) {
       Real δ = δ0;
       try {
         for (unsigned i = 0; i < N; i++) {
           std::vector<Vector<Scalar>> δz = generateDiscreteGradientδz(J, γ);
-          δ = 1.1 * perturb(J, δ, δz, γ);
+          double stepSize = 0;
+          for (const Vector<Scalar>& v : δz) {
+            stepSize += v.norm();
+          }
+          if (stepSize / δz.size() < 1e-6) {
+            break;
+          }
+          std::cout << cost(J) << " " << stepSize / δz.size() << std::endl;
+          δ = 2 * perturb(J, δ, δz, γ);
         }
       } catch (std::exception& e) {
       }
     }
 
-    template <int p, class Gen>
-    void relaxRandom(const Tensor<Scalar, p>& J, unsigned N, Real δ0, Gen& r) {
-      Real δ = δ0;
-        for (unsigned i = 0; i < N; i++) {
-          try {
-            std::vector<Vector<Scalar>> δz = generateRandomδz(r);
-            δ = 1.1 * perturb(J, δ, δz);
-          } catch (std::exception& e) {
-        }
-      }
-    }
-
     template <int p>
     Real cost(const Tensor<Scalar, p>& J, Real γ = 0) const {
       Real c = 0;
@@ -334,3 +207,180 @@ class Rope {
       return r;
     }
 };
+
+template <class Real, class Scalar, int p>
+bool stokesLineTest(const Tensor<Scalar, p>& J, const Vector<Scalar>& z1, const Vector<Scalar>& z2, unsigned n0, unsigned steps) {
+  Rope stokes(n0, z1, z2, J);
+
+  Real oldCost = stokes.cost(J);
+
+  for (unsigned i = 0; i < steps; i++) {
+    stokes.relaxDiscreteGradient(J, 1e6, 1, pow(2, steps));
+
+    Real newCost = stokes.cost(J);
+
+    if (newCost > oldCost) {
+      return false;
+    }
+
+    oldCost = newCost;
+
+    stokes = stokes.interpolate();
+  }
+  return true;
+}
+
+template <class Scalar>
+class Cord {
+public:
+  std::vector<Vector<Scalar>> gs;
+  Vector<Scalar> z0;
+  Vector<Scalar> z1;
+
+  template <int p>
+  Cord(const Tensor<Scalar, p>& J, const Vector<Scalar>& z2, const Vector<Scalar>& z3, unsigned ng) : gs(ng, Vector<Scalar>::Zero(z2.size())) {
+    Scalar H2 = getHamiltonian(J, z2);
+    Scalar H3 = getHamiltonian(J, z3);
+
+    if (real(H2) > real(H3)) {
+      z0 = z2;
+      z1 = z3;
+    } else {
+      z0 = z3;
+      z1 = z2;
+    }
+  }
+
+  Real gCoeff(unsigned i, Real t) const {
+    return (1 - t) * t * pow(t, i);
+  }
+
+  Real dgCoeff(unsigned i, Real t) const {
+    return (i + 1) * (1 - t) * pow(t, i) - pow(t, i + 1);
+  }
+
+  Vector<Scalar> f(Real t) const {
+    Vector<Scalar> z = (1 - t) * z0 + t * z1;
+
+    for (unsigned i = 0; i < gs.size(); i++) {
+      z += gCoeff(i, t) * gs[i];
+    }
+
+    return z;
+  }
+
+  Vector<Scalar> df(Real t) const {
+    Vector<Scalar> z = z1 - z0;
+
+    for (unsigned i = 0; i < gs.size(); i++) {
+      z += dgCoeff(i, t) * gs[i];
+    }
+
+    return z;
+  }
+
+  template <int p>
+  Real cost(const Tensor<Scalar, p>& J, Real t) const {
+    Vector<Scalar> z = f(t);
+    Scalar H;
+    Vector<Scalar> dH;
+    std::tie(H, dH, std::ignore) = hamGradHess(J, z);
+    Vector<Scalar> ż = zDot(z, dH);
+    Vector<Scalar> dz = df(t);
+
+    return 1 - real(ż.dot(dz)) / ż.norm() / dz.norm();
+  }
+
+  template <int p>
+  Real totalCost(const Tensor<Scalar, p>& J, unsigned nt) const {
+    Real tc = 0;
+
+    for (unsigned i = 0; i < nt; i++) {
+      Real t = (i + 1.0) / (nt + 1.0);
+      tc += cost(J, t);
+    }
+
+    return tc;
+  }
+
+  template <int p>
+  std::vector<Vector<Scalar>> dgs(const Tensor<Scalar, p>& J, Real t) const {
+    Vector<Scalar> z = f(t);
+    auto [H, dH, ddH] = hamGradHess(J, z);
+    Vector<Scalar> ż = zDot(z, dH);
+    Vector<Scalar> dz = df(t);
+    Matrix<Scalar> dż = dzDot(z, dH);
+    Matrix<Scalar> dżc = dzDotConjugate(z, dH, ddH);
+
+    std::vector<Vector<Scalar>> x;
+    x.reserve(gs.size());
+
+    for (unsigned i = 0; i < gs.size(); i++) {
+      Real fdg = gCoeff(i, t);
+      Real dfdg = dgCoeff(i, t);
+      Vector<Scalar> dC = - 0.5 / ż.norm() / dz.norm() * (
+          dfdg * ż.conjugate() + fdg * dżc * dz + fdg * dż * dz.conjugate()
+          - real(dz.dot(ż)) * (
+              dfdg * dz.conjugate() / dz.squaredNorm() +
+              fdg * (dżc * ż + dż * ż.conjugate()) / ż.squaredNorm()
+            )
+        );
+
+      x.push_back(dC.conjugate());
+    }
+
+    return x;
+  }
+
+  template <int p>
+  std::pair<Real, Real> relaxStep(const Tensor<Scalar, p>& J, unsigned nt, Real δ₀) {
+    std::vector<Vector<Scalar>> dgsTot(gs.size(), Vector<Scalar>::Zero(z0.size()));
+
+    for (unsigned i = 0; i < nt; i++) {
+      Real t = (i + 1.0) / (nt + 1.0);
+      std::vector<Vector<Scalar>> dgsI = dgs(J, t);
+
+      for (unsigned j = 0; j < gs.size(); j++) {
+        dgsTot[j] += dgsI[j] / nt;
+      }
+    }
+
+    Real stepSize = 0;
+    for (const Vector<Scalar>& dgi : dgsTot) {
+      stepSize += dgi.squaredNorm();
+    }
+    stepSize = sqrt(stepSize);
+
+    Cord cNew(*this);
+
+    Real δ = δ₀;
+    Real oldCost = totalCost(J, nt);
+    Real newCost = std::numeric_limits<Real>::infinity();
+
+    while (newCost > oldCost) {
+      for (unsigned i = 0; i < gs.size(); i++) {
+        cNew.gs[i] = gs[i] - δ * dgsTot[i];
+      }
+
+      newCost = cNew.totalCost(J, nt);
+
+      δ /= 2;
+    }
+
+    gs = cNew.gs;
+
+    return {2 * δ, stepSize};
+  }
+
+  template <int p>
+  void relax(const Tensor<Scalar, p>& J, unsigned nt, Real δ₀, unsigned maxSteps) {
+    Real δ = δ₀;
+    Real stepSize = std::numeric_limits<Real>::infinity();
+    unsigned steps = 0;
+    while (stepSize > 1e-7 && steps < maxSteps) {
+      std::tie(δ, stepSize) = relaxStep(J, nt, δ);
+      std::cout << totalCost(J, nt) << " " << δ << " " << stepSize << std::endl;
+      steps++;
+    }
+  }
+};