summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJaron Kent-Dobias <jaron@kent-dobias.com>2021-11-06 22:03:34 +0100
committerJaron Kent-Dobias <jaron@kent-dobias.com>2021-11-06 22:03:34 +0100
commit09bd9058e00b609d8cee2baa8023710ba74e3deb (patch)
treef35fbfba5fcc2d8998dad3743928b500af1749e5
parentc60430fc1c5d90ae06d1fd019257474c8f395bef (diff)
downloadcode-09bd9058e00b609d8cee2baa8023710ba74e3deb.tar.gz
code-09bd9058e00b609d8cee2baa8023710ba74e3deb.tar.bz2
code-09bd9058e00b609d8cee2baa8023710ba74e3deb.zip
Finished Levenburg-Marquardt for Stokes lines.
-rw-r--r--langevin.cpp29
-rw-r--r--p-spin.hpp23
-rw-r--r--stokes.hpp139
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
@@ -121,6 +121,29 @@ Tensor<Scalar, 3> ddzDot(const Vector<Scalar>& z, const Vector<Scalar>& dH) {
}
template <class Scalar>
+Tensor<Scalar, 3> dcdzDot(const Vector<Scalar>& z, const Vector<Scalar>& dH, const Matrix<Scalar>& ddH) {
+ Tensor<Scalar, 1> zT = Eigen::TensorMap<Tensor<const Scalar, 1>>(z.data(), {z.size()});
+ Tensor<Scalar, 1> dHT = Eigen::TensorMap<Tensor<const Scalar, 1>>(dH.data(), {dH.size()});
+ Tensor<Scalar, 2> ddHT = Eigen::TensorMap<Tensor<const Scalar, 2>>(ddH.data(), {dH.size(), dH.size()});
+
+ Matrix<Scalar> I = Matrix<Real>::Identity(z.size(), z.size());
+ Tensor<Scalar, 2> IT = Eigen::TensorMap<Tensor<const Scalar, 2>>(I.data(), {z.size(), z.size()});
+
+ Eigen::array<Eigen::IndexPair<int>, 0> ei = {};
+
+ Scalar z² = z.squaredNorm();
+
+ return ddHT.conjugate().contract(zT.conjugate(), ei) / z²
+ +IT.contract(dHT.conjugate(), ei).shuffle((std::array<int, 3>){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<int, 3>){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 <class Scalar>
Tensor<Scalar, 3> ddzDotConjugate(const Vector<Scalar>& z, const Vector<Scalar>& dH, const Matrix<Scalar>& ddH, const Tensor<Scalar, 3>& dddH) {
Tensor<Scalar, 1> zT = Eigen::TensorMap<Tensor<const Scalar, 1>>(z.data(), {z.size()});
Tensor<Scalar, 1> dHT = Eigen::TensorMap<Tensor<const Scalar, 1>>(dH.data(), {dH.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<Scalar> dż = dzDot(z, dH);
Matrix<Scalar> dżc = dzDotConjugate(z, dH, ddH);
Tensor<Scalar, 3> ddż = ddzDot(z, dH);
+ Tensor<Scalar, 3> dcdż = dcdzDot(z, dH, ddH);
Tensor<Scalar, 3> ddżc = ddzDotConjugate(z, dH, ddH, dddH);
Eigen::array<Eigen::IndexPair<int>, 1> ip20 = {Eigen::IndexPair<int>(2, 0)};
@@ -352,7 +353,7 @@ public:
Matrix<Scalar> ddżcdz = Eigen::Map<Matrix<Scalar>>(ddżcdzT.data(), z.size(), z.size());
Tensor<Scalar, 2> ddżdzT = ddż.contract(dzT.conjugate(), ip20);
- Matrix<Scalar> ddżdz = Eigen::Map<Matrix<Scalar>>(ddżcdzT.data(), z.size(), z.size());
+ Matrix<Scalar> ddżdz = Eigen::Map<Matrix<Scalar>>(ddżdzT.data(), z.size(), z.size());
Tensor<Scalar, 2> ddżcżT = ddżc.contract(żT, ip20);
Matrix<Scalar> ddżcż = Eigen::Map<Matrix<Scalar>>(ddżcżT.data(), z.size(), z.size());
@@ -360,27 +361,35 @@ public:
Tensor<Scalar, 2> ddżżcT = ddż.contract(żT.conjugate(), ip20);
Matrix<Scalar> ddżżc = Eigen::Map<Matrix<Scalar>>(ddżżcT.data(), z.size(), z.size());
- Vector<Scalar> x(z.size() * gs.size());
+ Tensor<Scalar, 2> dcdżcdzT = dcdż.conjugate().contract(dzT, ip20);
+ Matrix<Scalar> dcdżcdz = Eigen::Map<Matrix<Scalar>>(dcdżcdzT.data(), z.size(), z.size());
+
+ Tensor<Scalar, 2> dcdżcżT = dcdż.conjugate().contract(żT, ip20);
+ Matrix<Scalar> dcdżcż = Eigen::Map<Matrix<Scalar>>(dcdżcżT.data(), z.size(), z.size());
+
+ Vector<Scalar> x(2 * z.size() * gs.size());
Matrix<Scalar> 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<Scalar> 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<Scalar> 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<Scalar> ddC =
- - 0.5 / ż.norm() / dz.norm() * (
- dfdg * fdgm * dżc.transpose() + dfdgm * fdg * dżc +
- fdg * fdgm * ddżcdz + fdg * fdgm * ddżdz
+
+ Matrix<Scalar> 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<Scalar> 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<Scalar>::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<Scalar> Mf = Matrix<Scalar>::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<Real>(0.03145, 0.08452);
- Scalar hj = std::complex<Real>(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 <int p>
std::pair<Real, Real> relaxStepNewton(const Tensor<Scalar, p>& J, unsigned nt, Real δ₀) {
- Vector<Scalar> dgsTot = Vector<Scalar>::Zero(z0.size() * gs.size());
- Matrix<Scalar> HessTot = Matrix<Scalar>::Zero(z0.size() * gs.size(), z0.size() * gs.size());
+ Vector<Scalar> dgsTot = Vector<Scalar>::Zero(2 * z0.size() * gs.size());
+ Matrix<Scalar> HessTot = Matrix<Scalar>::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<Real>::infinity();
- Vector<Scalar> δ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<Scalar> I = abs(HessTot.diagonal().array()).matrix().asDiagonal();
while (newCost > oldCost) {
+ Vector<Scalar> δ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 <int p>
@@ -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++;
}
}