diff options
-rw-r--r-- | topology.cpp | 58 |
1 files changed, 31 insertions, 27 deletions
diff --git a/topology.cpp b/topology.cpp index 6650ea6..1d46649 100644 --- a/topology.cpp +++ b/topology.cpp @@ -103,15 +103,40 @@ public: return {∂L, ∂∂L}; } - Vector newtonMethod(const Vector& v0, Real γ = 1, Real ε = 1e-12) const { + Vector newtonMethod(const Vector& v0, Real γ = 1, Real ε = 1e-12, unsigned maxSteps = 1e2) const { Vector v = v0; + unsigned steps = 0; Vector ∂L; Matrix ∂∂L; while (std::tie(∂L, ∂∂L) = ∂L_∂∂L(v), ∂L.squaredNorm() > ε) { + if (v.tail(M + 1).squaredNorm() > 1 / ε || steps > N * maxSteps) { + return Vector::Zero(N + M + 1); + } + v -= γ * ∂∂L.partialPivLu().solve(∂L); - v.head(N) = normalize(v.head(N)); // might as well stay on the sphere + + v.head(N) = normalize(v.head(N)); // Stay on the sphere + steps++; + } + + return v; + } + + Vector randomStationaryPoint(Rng& r, Real γ = 1, Real ε = 1e-12, unsigned maxTries = 1e1) const { + Vector v = Vector::Zero(N + M + 1); + unsigned tries = 0; + + while (v.squaredNorm() < ε && tries < maxTries) { + Vector v0(N + M + 1); + for (Real& v0ᵢ : v0) { + v0ᵢ = r.variate<Real, std::normal_distribution>(); + } + v0.head(N) = normalize(v0.head(N)); + + v = newtonMethod(v0, γ, ε); + tries++; } return v; @@ -235,30 +260,15 @@ int main(int argc, char* argv[]) { std::cout << std::setprecision(15); if (quadratic) { - Vector x = Vector::Zero(N); - x(0) = sqrt(N); - for (unsigned sample = 0; sample < samples; sample++) { - Vector v0(N + M + 1); - v0.head(N) = x; - for (Real& v0ᵢ : v0.tail(M + 1)) { - v0ᵢ = r.variate<Real, std::normal_distribution>(); - } ExtensiveQuadratic S(N, M, E, r); - Vector v = S.newtonMethod(v0, γ); - std::cout << S.overlap(v) << " " << S.cost(v) << std::endl; + Vector v = S.randomStationaryPoint(r, γ); + std::cout << S.overlap(v) << std::endl; } } else if (!count) { - Vector x = Vector::Zero(N); - x(0) = sqrt(N); - for (unsigned sample = 0; sample < samples; sample++) { - Vector v0(N + 2); - v0 << x, - r.variate<Real, std::normal_distribution>(), - r.variate<Real, std::normal_distribution>(); Spherical3Spin S(N, E, r); - Vector v = S.newtonMethod(v0, γ); + Vector v = S.randomStationaryPoint(r, γ); std::cout << S.overlap(v) << std::endl; } } else { @@ -266,13 +276,7 @@ int main(int argc, char* argv[]) { std::list<Vector> points; for (unsigned sample = 0; sample < samples; sample++) { - Vector v0(N + 2); - for (Real& v0i : v0) { - v0i = r.variate<Real, std::normal_distribution>(); - } - v0.head(N) = normalize(v0.head(N)); - - Vector v = S.newtonMethod(v0, γ); + Vector v = S.randomStationaryPoint(r, γ); if (v(N) == v(N)) { bool inList = false; for (const Vector& u : points) { |