From 09bd9058e00b609d8cee2baa8023710ba74e3deb Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Sat, 6 Nov 2021 22:03:34 +0100 Subject: Finished Levenburg-Marquardt for Stokes lines. --- langevin.cpp | 29 ++++++++++++- p-spin.hpp | 23 ++++++++++ stokes.hpp | 139 +++++++++++++++++++++++++++++++++-------------------------- 3 files changed, 128 insertions(+), 63 deletions(-) diff --git a/langevin.cpp b/langevin.cpp index d23cb33..adb5212 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -193,8 +193,35 @@ int main(int argc, char* argv[]) { } */ + /* + J(0,0,0) = Complex(2, 3); + J(1,1,1) = Complex(-2, 0.3); + J(0,1,1) = Complex(4, 0.2); + J(1,0,1) = Complex(4, 0.2); + J(1,1,0) = Complex(4, 0.2); + J(1,0,0) = Complex(0.7, 0.4); + J(0,1,0) = Complex(0.7, 0.4); + J(0,0,1) = Complex(0.7, 0.4); + + ComplexVector z0(2);; + z0 << Complex(0.8, 0.3), Complex(0.7, 0.2); + ComplexVector z1(2); + z1 << Complex(-0.5, 0.2), Complex(1.0, 1.0); + + Cord test(J, z0, z1, 2); + + test.gs[0](0) = Complex(0.2, 0.2); + test.gs[0](1) = Complex(0.4, 0.4); + test.gs[1](0) = Complex(0.1, 0.2); + test.gs[1](1) = Complex(0.3, 0.4); + + auto [dgs, ddgs] = test.gsGradHess(J, 0.7); + + std::cout << dgs << std::endl; + std::cout << ddgs << std::endl; + */ - Cord test(J, zSaddle, zSaddleNext, 2); + Cord test(J, zSaddle, zSaddleNext, 5); test.relaxNewton(J, 20, 1, 1e4); std::cout << test.z0.transpose() << std::endl; diff --git a/p-spin.hpp b/p-spin.hpp index 50da209..15d2525 100644 --- a/p-spin.hpp +++ b/p-spin.hpp @@ -120,6 +120,29 @@ Tensor ddzDot(const Vector& z, const Vector& dH) { + zT.conjugate().contract(zT.conjugate(), ei).contract(zT.conjugate(), ei) * (2.0 * dH.dot(z) / pow(z², 3)); } +template +Tensor dcdzDot(const Vector& z, const Vector& dH, const Matrix& ddH) { + Tensor zT = Eigen::TensorMap>(z.data(), {z.size()}); + Tensor dHT = Eigen::TensorMap>(dH.data(), {dH.size()}); + Tensor ddHT = Eigen::TensorMap>(ddH.data(), {dH.size(), dH.size()}); + + Matrix I = Matrix::Identity(z.size(), z.size()); + Tensor IT = Eigen::TensorMap>(I.data(), {z.size(), z.size()}); + + Eigen::array, 0> ei = {}; + + Scalar z² = z.squaredNorm(); + + return ddHT.conjugate().contract(zT.conjugate(), ei) / z² + +IT.contract(dHT.conjugate(), ei).shuffle((std::array){0,2,1}) / z² + -zT.contract(dHT.conjugate(), ei).contract(zT.conjugate(), ei) / pow(z², 2) + -ddHT.conjugate().contract(zT, ip00).contract(zT.conjugate(), ei).contract(zT.conjugate(), ei) / pow(z², 2) + -IT.contract(zT.conjugate(), ei) * (dH.dot(z) / pow(z², 2)) + -IT.contract(zT.conjugate(), ei).shuffle((std::array){0,2,1}) * (dH.dot(z) / pow(z², 2)) + +zT.contract(zT.conjugate(), ei).contract(zT.conjugate(), ei) * (2.0 * dH.dot(z) / pow(z², 3)) + ; +} + template Tensor ddzDotConjugate(const Vector& z, const Vector& dH, const Matrix& ddH, const Tensor& dddH) { Tensor zT = Eigen::TensorMap>(z.data(), {z.size()}); diff --git a/stokes.hpp b/stokes.hpp index 948a484..e14db52 100644 --- a/stokes.hpp +++ b/stokes.hpp @@ -344,6 +344,7 @@ public: Matrix dż = dzDot(z, dH); Matrix dżc = dzDotConjugate(z, dH, ddH); Tensor ddż = ddzDot(z, dH); + Tensor dcdż = dcdzDot(z, dH, ddH); Tensor ddżc = ddzDotConjugate(z, dH, ddH, dddH); Eigen::array, 1> ip20 = {Eigen::IndexPair(2, 0)}; @@ -352,7 +353,7 @@ public: Matrix ddżcdz = Eigen::Map>(ddżcdzT.data(), z.size(), z.size()); Tensor ddżdzT = ddż.contract(dzT.conjugate(), ip20); - Matrix ddżdz = Eigen::Map>(ddżcdzT.data(), z.size(), z.size()); + Matrix ddżdz = Eigen::Map>(ddżdzT.data(), z.size(), z.size()); Tensor ddżcżT = ddżc.contract(żT, ip20); Matrix ddżcż = Eigen::Map>(ddżcżT.data(), z.size(), z.size()); @@ -360,27 +361,35 @@ public: Tensor ddżżcT = ddż.contract(żT.conjugate(), ip20); Matrix ddżżc = Eigen::Map>(ddżżcT.data(), z.size(), z.size()); - Vector x(z.size() * gs.size()); + Tensor dcdżcdzT = dcdż.conjugate().contract(dzT, ip20); + Matrix dcdżcdz = Eigen::Map>(dcdżcdzT.data(), z.size(), z.size()); + + Tensor dcdżcżT = dcdż.conjugate().contract(żT, ip20); + Matrix dcdżcż = Eigen::Map>(dcdżcżT.data(), z.size(), z.size()); + + Vector x(2 * z.size() * gs.size()); Matrix M(x.size(), x.size()); for (unsigned i = 0; i < gs.size(); i++) { - Real fdg = gCoeff(i, t); - Real dfdg = dgCoeff(i, t); + Real fdgn = gCoeff(i, t); + Real dfdgn = dgCoeff(i, t); Vector dC = - 0.5 / ż.norm() / dz.norm() * ( - dfdg * ż.conjugate() + fdg * dżc * dz + fdg * dż * dz.conjugate() + dfdgn * ż.conjugate() + fdgn * dżc * dz + fdgn * dż * dz.conjugate() - real(dz.dot(ż)) * ( - dfdg * dz.conjugate() / dz.squaredNorm() + - fdg * (dżc * ż + dż * ż.conjugate()) / ż.squaredNorm() + dfdgn * dz.conjugate() / dz.squaredNorm() + + fdgn * (dżc * ż + dż * ż.conjugate()) / ż.squaredNorm() ) ); for (unsigned j = 0; j < dC.size(); j++) { x(z.size() * i + j) = dC(j); + x(z.size() * gs.size() + z.size() * i + j) = conj(dC(j)); } for (unsigned j = 0; j < gs.size(); j++) { Real fdgm = gCoeff(j, t); Real dfdgm = dgCoeff(j, t); + Scalar CC = - real(dz.dot(ż)) / ż.norm() / dz.norm(); Vector dCm = - 0.5 / ż.norm() / dz.norm() * ( dfdgm * ż.conjugate() + fdgm * dżc * dz + fdgm * dż * dz.conjugate() - real(dz.dot(ż)) * ( @@ -388,66 +397,80 @@ public: fdgm * (dżc * ż + dż * ż.conjugate()) / ż.squaredNorm() ) ); - Matrix ddC = - - 0.5 / ż.norm() / dz.norm() * ( - dfdg * fdgm * dżc.transpose() + dfdgm * fdg * dżc + - fdg * fdgm * ddżcdz + fdg * fdgm * ddżdz + + Matrix ddC = 0.5 / ż.norm() / dz.norm() * ( + - ( + dfdgn * fdgm * dżc.transpose() + dfdgm * fdgn * dżc + + fdgn * fdgm * ddżcdz + fdgn * fdgm * ddżdz ) - + 0.5 / dz.squaredNorm() * ( - dfdg * dz.conjugate() * dCm.transpose() + + - ż.norm() / dz.norm() * ( + dfdgn * dz.conjugate() * dCm.transpose() + dfdgm * dC * dz.adjoint() ) - + 0.5 / ż.squaredNorm() * ( - fdg * (dżc * ż + dż * ż.conjugate()) * dCm.transpose() + + - dz.norm() / ż.norm() * ( + fdgn * (dżc * ż + dż * ż.conjugate()) * dCm.transpose() + fdgm * dC * (dżc * ż + dż * ż.conjugate()).transpose() ) - + 0.5 / ż.norm() / dz.norm() * real(dz.dot(ż)) * fdg * fdgm * ( + - CC * dz.norm() / ż.norm() * fdgn * fdgm * ( ddżżc + ddżcż + dżc * dż.transpose() + dż * dżc.transpose() - ) / ż.squaredNorm() - - 0.5 / ż.norm() / dz.norm() * real(dz.dot(ż)) * dfdg * dfdgm * dz.conjugate() * dz.adjoint() / pow(dz.squaredNorm(), 2) - - 0.5 / ż.norm() / dz.norm() * real(dz.dot(ż)) / pow(ż.squaredNorm(), 2) - * fdg * (dżc * ż + dż * ż.conjugate()) + ) + - 0.5 * CC / dz.norm() / ż.norm() * ( + dfdgn * fdgm * dz.conjugate() * (dżc * ż + dż * ż.conjugate()).transpose() + + dfdgm * fdgn * (dżc * ż + dż * ż.conjugate()) * dz.adjoint() + ) + + 0.5 * CC * ż.norm() / pow(dz.norm(), 3) * dfdgn * dfdgm * dz.conjugate() * dz.adjoint() + + 0.5 * CC * dz.norm() / pow(ż.norm(), 3) + * fdgn * (dżc * ż + dż * ż.conjugate()) * fdgm * (dżc * ż + dż * ż.conjugate()).transpose() + ) + ; + + Matrix dcdC = 0.5 / ż.norm() / dz.norm() * ( + - ( + dfdgn * fdgm * dż + dfdgm * fdgn * dż.adjoint() + + fdgn * fdgm * (dcdżcdz + dcdżcdz.adjoint()) + ) + - ż.norm() / dz.norm() * ( + dfdgn * dz.conjugate() * dCm.adjoint() + + dfdgm * dC * dz.transpose() + ) + - dz.norm() / ż.norm() * ( + fdgn * (dżc * ż + dż * ż.conjugate()) * dCm.adjoint() + + fdgm * dC * (dżc * ż + dż * ż.conjugate()).adjoint() + ) + - CC * ż.norm() / dz.norm() * dfdgn * dfdgm * Matrix::Identity(z.size(), z.size()) + - CC * dz.norm() / ż.norm() * fdgn * fdgm * ( + dcdżcż + dcdżcż.adjoint() + dżc * dżc.adjoint() + dż * dż.adjoint() + ) + - 0.5 * CC / dz.norm() / ż.norm() * ( + dfdgn * fdgm * dz.conjugate() * (dżc * ż + dż * ż.conjugate()).adjoint() + + dfdgm * fdgn * (dżc * ż + dż * ż.conjugate()) * dz.transpose() + ) + + 0.5 * CC * ż.norm() / pow(dz.norm(), 3) * dfdgn * dfdgm * dz.conjugate() * dz.transpose() + + 0.5 * CC * dz.norm() / pow(ż.norm(), 3) * fdgn * fdgm + * (dżc * ż + dż * ż.conjugate()) + * (dżc * ż + dż * ż.conjugate()).adjoint() + ) ; for (unsigned k = 0; k < z.size(); k++) { for (unsigned l = 0; l < z.size(); l++) { - M(z.size() * i + k, z.size() * j + l) = ddC(k, l); + M(z.size() * i + k, z.size() * j + l) = dcdC(k, l); + M(z.size() * gs.size() + z.size() * i + k, z.size() * j + l) = conj(ddC(k, l)); + M(z.size() * i + k, z.size() * gs.size() + z.size() * j + l) = ddC(k, l); + M(z.size() * gs.size() + z.size() * i + k, z.size() * gs.size() + z.size() * j + l) = conj(dcdC(k, l)); } } } } - Cord cNew(*this); - Matrix Mf = Matrix::Zero(x.size(), x.size()); - for (unsigned i = 0; i < x.size(); i++) { - for (unsigned j = 0; j < x.size(); j++) { - Scalar hi = std::complex(0.03145, 0.08452); - Scalar hj = std::complex(0.09483, 0.02453); - Scalar c0 = cost(J, t); - cNew.gs[i / z.size()](i % z.size()) += hi; - Scalar ci = cNew.cost(J, t); - cNew.gs[j / z.size()](j % z.size()) += hj; - Scalar cij = cNew.cost(J, t); - cNew.gs[i / z.size()](i % z.size()) -= hi; - Scalar cj = cNew.cost(J, t); - cNew.gs[j / z.size()](j % z.size()) -= hj; - - Mf(i, j) = (cij - ci - cj + c0) / (hi * hj); - } - } - - std::cout << M << std::endl; - std::cout << Mf << std::endl; - getchar(); - return {x, M}; } template std::pair relaxStepNewton(const Tensor& J, unsigned nt, Real δ₀) { - Vector dgsTot = Vector::Zero(z0.size() * gs.size()); - Matrix HessTot = Matrix::Zero(z0.size() * gs.size(), z0.size() * gs.size()); + Vector dgsTot = Vector::Zero(2 * z0.size() * gs.size()); + Matrix HessTot = Matrix::Zero(2 * z0.size() * gs.size(), 2 * z0.size() * gs.size()); for (unsigned i = 0; i < nt; i++) { Real t = (i + 1.0) / (nt + 1.0); @@ -463,35 +486,27 @@ public: Real δ = δ₀; Real oldCost = totalCost(J, nt); + Real newCost = std::numeric_limits::infinity(); - Vector δg = HessTot.partialPivLu().solve(dgsTot); - - for (unsigned i = 0; i < gs.size(); i++) { - for (unsigned j = 0; j < z0.size(); j++) { - cNew.gs[i](j) = gs[i](j) - δg[z0.size() * i + j]; - } - } - - Real newCost = cNew.totalCost(J, nt); - if (newCost < oldCost) { - std::cout << "Newton step accepted!" << std::endl; - } + Matrix I = abs(HessTot.diagonal().array()).matrix().asDiagonal(); while (newCost > oldCost) { + Vector δg = (HessTot + δ * I).partialPivLu().solve(dgsTot); + for (unsigned i = 0; i < gs.size(); i++) { for (unsigned j = 0; j < z0.size(); j++) { - cNew.gs[i](j) = gs[i](j) - δ * conj(dgsTot[z0.size() * i + j]); + cNew.gs[i](j) = gs[i](j) - δg[z0.size() * gs.size() + z0.size() * i + j]; } } newCost = cNew.totalCost(J, nt); - δ /= 2; + δ *= 2; } gs = cNew.gs; - return {2*δ, stepSize}; + return {δ / 2, stepSize}; } template @@ -555,7 +570,7 @@ public: while (stepSize > 1e-7 && steps < maxSteps) { std::tie(δ, stepSize) = relaxStepNewton(J, nt, δ); std::cout << totalCost(J, nt) << " " << δ << " " << stepSize << std::endl; - δ *= 2; + δ /= 3; steps++; } } -- cgit v1.2.3-70-g09d2