summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJaron Kent-Dobias <jaron@kent-dobias.com>2024-04-12 11:41:22 +0200
committerJaron Kent-Dobias <jaron@kent-dobias.com>2024-04-12 11:41:22 +0200
commitc0636131197d2d48de5b42b42411e16c55429d7d (patch)
tree6827a359c94bc6973a106c9720327f186b91039a
parent4bdde3b2c1b6ea5e828a987211437f30587356ae (diff)
downloadcode-c0636131197d2d48de5b42b42411e16c55429d7d.tar.gz
code-c0636131197d2d48de5b42b42411e16c55429d7d.tar.bz2
code-c0636131197d2d48de5b42b42411e16c55429d7d.zip
Some small optimizations, and better judging of the endpoint for minimization.
-rw-r--r--least_squares.cpp38
1 files changed, 16 insertions, 22 deletions
diff --git a/least_squares.cpp b/least_squares.cpp
index 7b9c2c9..f69f117 100644
--- a/least_squares.cpp
+++ b/least_squares.cpp
@@ -40,7 +40,7 @@ Vector normalize(const Vector& x) {
}
Vector makeTangent(const Vector& v, const Vector& x) {
- return v - (v.dot(x) / x.squaredNorm()) * x;
+ return v - (v.dot(x) / x.size()) * x;
}
Real HFromV(const Vector& V) {
@@ -167,31 +167,24 @@ Vector gradientDescent(const QuadraticModel& M, const Vector& x0, Real ε = 1e-7
return x;
}
-Vector findMinimum(const QuadraticModel& M, const Vector& x0, Real ε = 1e-12) {
+Vector findMinimum(const QuadraticModel& M, const Vector& x0, Real ε = 1e-5) {
Vector x = x0;
Real λ = 100;
auto [H, g, m] = M.hamGradHess(x0);
- while (g.norm() / M.N > ε && λ * ε < 1) {
- while (λ * ε < 1) {
- Vector dx = (m - λ * (Matrix)m.diagonal().cwiseAbs().asDiagonal()).partialPivLu().solve(g);
- dx -= x.dot(dx) * x / M.N;
- Vector xNew = normalize(x - dx);
+ while (λ * ε < 1) {
+ Vector dx = (m - λ * (Matrix)m.diagonal().cwiseAbs().asDiagonal()).partialPivLu().solve(g);
+ Vector xNew = normalize(x - makeTangent(dx, x));
+ Real HNew = M.getHamiltonian(xNew);
- auto [HNew, gNew, mNew] = M.hamGradHess(xNew);
+ if (HNew > H) {
+ x = xNew;
+ std::tie(H, g, m) = M.hamGradHess(xNew);
- if (HNew > H) {
- x = xNew;
- H = HNew;
- g = gNew;
- m = mNew;
-
- λ /= 2;
- break;
- } else {
- λ *= 1.5;
- }
+ λ /= 2;
+ } else {
+ λ *= 1.5;
}
}
@@ -301,10 +294,11 @@ int main(int argc, char* argv[]) {
for (unsigned sample = 0; sample < samples; sample++) {
QuadraticModel leastSquares(N, M, r, σ, A, J);
- x = metropolisSweep(leastSquares, x, β, r, sweeps);
- std::cout << " " << leastSquares.getHamiltonian(x) / N;
+ if (β != 0) {
+ x = metropolisSweep(leastSquares, x, β, r, sweeps);
+ }
x = findMinimum(leastSquares, x);
- std::cout << " " << leastSquares.getHamiltonian(x) / N;
+ std::cout << " " << leastSquares.getHamiltonian(x) / N << std::flush;
}
std::cout << std::endl;