From a20d6c31acbcca2b952cac7553fb7744dfaaf728 Mon Sep 17 00:00:00 2001 From: Christoph Lehner Date: Mon, 11 Dec 2023 12:16:10 +0100 Subject: [PATCH] matrix sqrt and log --- lib/gpt/core/matrix/__init__.py | 3 +- lib/gpt/core/matrix/inv.py | 11 +++++-- lib/gpt/core/matrix/log.py | 53 +++++++++++++++++---------------- lib/gpt/core/matrix/sqrt.py | 45 ++++++++++++++++++++++++++++ tests/core/matrix.py | 18 +++++++---- tests/qcd/gauge.py | 20 +++++++++++++ 6 files changed, 115 insertions(+), 35 deletions(-) create mode 100644 lib/gpt/core/matrix/sqrt.py diff --git a/lib/gpt/core/matrix/__init__.py b/lib/gpt/core/matrix/__init__.py index 1f21ab503..f39d24cc4 100644 --- a/lib/gpt/core/matrix/__init__.py +++ b/lib/gpt/core/matrix/__init__.py @@ -16,7 +16,8 @@ # with this program; if not, write to the Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # +from gpt.core.matrix.inv import inv +from gpt.core.matrix.sqrt import sqrt from gpt.core.matrix.exp import exp from gpt.core.matrix.log import log -from gpt.core.matrix.inv import inv from gpt.core.matrix.det import det diff --git a/lib/gpt/core/matrix/inv.py b/lib/gpt/core/matrix/inv.py index 8c2f899f1..3f175e1f9 100644 --- a/lib/gpt/core/matrix/inv.py +++ b/lib/gpt/core/matrix/inv.py @@ -22,7 +22,14 @@ def inv(A): A = gpt.eval(A) assert isinstance(A, gpt.lattice) - A_inv = gpt.lattice(A) + to_list = gpt.util.to_list - cgpt.invert_matrix(to_list(A_inv), to_list(A)) + + Al = to_list(A) + + if Al[0].otype.shape == (1,): + return gpt.component.inv(A) + + A_inv = gpt.lattice(A) + cgpt.invert_matrix(to_list(A_inv), Al) return A_inv diff --git a/lib/gpt/core/matrix/log.py b/lib/gpt/core/matrix/log.py index caae21e41..91e0bde51 100644 --- a/lib/gpt/core/matrix/log.py +++ b/lib/gpt/core/matrix/log.py @@ -19,32 +19,33 @@ import gpt, numpy -def log(i, convergence_threshold=0.5): - i = gpt.eval(i) - # i = n*(1 + x), log(i) = log(n) + log(1+x) - # x = i/n - 1, |x|^2 = = |i|^2/n^2 + |1|^2 - ( + <1,i>)/n - # d/dn |x|^2 = -2 |i|^2/n^3 + ( + <1,i>)/n^2 = 0 -> 2|i|^2 == n ( + <1,i>) - if i.grid.precision != gpt.double: - x = gpt.convert(i, gpt.double) - else: - x = gpt.copy(i) - lI = gpt.identity(x.new()) - n = gpt.norm2(x) / gpt.inner_product(x, lI).real - x /= n - x -= lI - n2 = gpt.norm2(x) ** 0.5 / x.grid.gsites - order = 8 * int(16 / (-numpy.log10(n2))) - # n2 = (gpt.norm2(x) / x.grid.gsites / x.otype.shape[0]) ** 0.5 - # order = int(16 / (-numpy.log10(n2))) - assert n2 < convergence_threshold +def log(i, convergence_threshold=0.0001): + x = gpt.eval(i) + + if x.grid.precision is gpt.single: + # perform intermediate calculation with enhanced precision + return gpt.convert(log(gpt.convert(x, gpt.double)), gpt.single) + + Id = gpt.identity(x) + + nrm = gpt.norm2(Id) + + # log(x^{1/2^n}) = 1/2^n log(x) + scale = 1.0 + for n in range(20): + x = gpt.matrix.sqrt(x) + + eps2 = gpt.norm2(x - Id) / nrm + scale *= 2.0 + if eps2 < convergence_threshold: + break + + x = gpt(x - Id) o = gpt.copy(x) - xn = gpt.copy(x) - for j in range(2, order + 1): - xn @= xn * x + xn = x + for j in range(2, 8): + xn = gpt(xn * x) o -= xn * ((-1.0) ** j / j) - o += lI * numpy.log(n) - if i.grid.precision != gpt.double: - r = gpt.lattice(i) - gpt.convert(r, o) - o = r + + o = gpt(scale * o) return o diff --git a/lib/gpt/core/matrix/sqrt.py b/lib/gpt/core/matrix/sqrt.py new file mode 100644 index 000000000..4c4cbc8e4 --- /dev/null +++ b/lib/gpt/core/matrix/sqrt.py @@ -0,0 +1,45 @@ +# +# GPT - Grid Python Toolkit +# Copyright (C) 2020-22 Christoph Lehner (christoph.lehner@ur.de, https://github.com/lehner/gpt) +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. +# +import gpt + + +def sqrt(A): + # Denman and Beavers 1976 + Mk = A + Xk = gpt.identity(A) + norm = gpt.norm2(Xk) + for k in range(100): + Xkp1 = gpt(0.5 * Xk + 0.5 * gpt.matrix.inv(Mk)) + Mkp1 = gpt(0.5 * Mk + 0.5 * gpt.matrix.inv(Xk)) + + eps = (gpt.norm2(Mk - Mkp1) / norm) ** 0.5 + Xk = Xkp1 + Mk = Mkp1 + + if eps < 1e3 * A.grid.precision.eps: + # gpt.message(f"Converged after {k} iterations") + # one final Newton iteration (including the original A to avoid error accumulation) + Mk = gpt(0.5 * Mk + 0.5 * gpt.matrix.inv(Mk) * A) + + # compute error + # gpt.message("err",gpt.norm2(Mk * Mk - A) / gpt.norm2(A)) + return Mk + + gpt.message("Warning: sqrt not converged") + return Mk diff --git a/tests/core/matrix.py b/tests/core/matrix.py index 4ec27f185..bb69f7a75 100755 --- a/tests/core/matrix.py +++ b/tests/core/matrix.py @@ -70,26 +70,32 @@ def mod0p1(x): g.message( f""" - Test log,exp,det,tr for {grid.precision.__name__} + Test sqrt,log,exp,det,tr for {grid.precision.__name__} """ ) for dtype in [g.mspincolor, g.mcolor, g.mspin, lambda grid: g.mcomplex(grid, 8)]: rng = g.random("test") m = rng.cnormal(dtype(grid)) + + sqrt_m = g.matrix.sqrt(m) + eps2 = g.norm2(sqrt_m * sqrt_m - m) / g.norm2(m) + g.message(f"test sqrt(M)*sqrt(M) = M for {m.otype.__name__}: {eps2}") + assert eps2 < eps**2 + minv = g.matrix.inv(m) eye = g.identity(m) - eps2 = g.norm2(m * minv - eye) / (12 * grid.fsites) + eps2 = g.norm2(m * minv - eye) / g.norm2(eye) g.message(f"test M*M^-1 = 1 for {m.otype.__name__}: {eps2}") assert eps2 < eps**2 - # make logarithm well defined - m @= eye + 0.01 * m m2 = g.matrix.exp(g.matrix.log(m)) eps2 = g.norm2(m - m2) / g.norm2(m) g.message(f"exp(log(m)) == m: {eps2}") - assert eps2 < eps**2.0 + assert eps2 < eps**2.0 * 1e3 + # make inverse well defined + m @= eye + 0.01 * m eps2 = g.norm2(g.matrix.log(g.matrix.det(g.matrix.exp(m))) - g.trace(m)) / g.norm2(m) g.message(f"log(det(exp(m))) == tr(m): {eps2}") - assert eps2 < eps**2.0 + assert eps2 < eps**2.0 * 1e3 diff --git a/tests/qcd/gauge.py b/tests/qcd/gauge.py index b8709c536..c789a808b 100755 --- a/tests/qcd/gauge.py +++ b/tests/qcd/gauge.py @@ -222,3 +222,23 @@ eps1 = g.norm2(f_test.gradient(V1, V1)) ** 0.5 / f_test(V1) g.message(f"df/f after {tag} gauge fix: {eps1}, improvement: {eps1/eps0}") assert eps1 / eps0 < expected_improvement + + +# test temporal gauge +def identity(U, mu=3): + # U'(x) = V(x) U_mu(x) Vdag(x+mu) + # V = [1, U[0], U[0] U[1], U[0] U[1] U[2], ...] + U_n = g.separate(U[mu], mu) + V_n = [g.identity(U_n[0])] + N = len(U_n) + for n in range(N - 1): + V_n.append(g(V_n[n] * U_n[n])) + return g.merge(V_n, mu) + + +V = identity(U) + +Up = g.qcd.gauge.transformed(U, V) + +for t in range(16): + print(t, Up[3][0, 0, 0, t])