matrix sqrt and log

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

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

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/n - 1, i/n - 1> = |i|^2/n^2 + |1|^2 - (<i,1> + <1,i>)/n
-    # d/dn |x|^2 = -2 |i|^2/n^3 + (<i,1> + <1,i>)/n^2 = 0 -> 2|i|^2 == n (<i,1> + <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

lib/gpt/core/matrix/sqrt.py

@@ -0,0 +1,45 @@
+#
+#    GPT - Grid Python Toolkit
+#    Copyright (C) 2020-22  Christoph Lehner ([email protected], 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

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

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])