Fix map_eigvals for fermionic case and test

Author: emstoudenmire

src/lib/ITensorsExtensions/src/itensor.jl

@@ -1,5 +1,6 @@
 using LinearAlgebra: LinearAlgebra, eigen, pinv
 using ITensors:
+  ITensors,
   ITensor,
   Index,
   commonind,
@@ -58,18 +59,50 @@ function eigendecomp(A::ITensor, linds, rinds; ishermitian=false, kwargs...)
   D, U = eigen(A, linds, rinds; ishermitian, kwargs...)
   ul, ur = noncommonind(D, U), commonind(D, U)
   Ul = replaceinds(U, vcat(rinds, ur), vcat(linds, ul))
-
   return Ul, D, dag(U)
 end
 
-function map_eigvals(f::Function, A::ITensor, inds...; ishermitian=false, kwargs...)
+function map_eigvals(f::Function, A::ITensor, Linds, Rinds; kws...)
+
+  # <fermions>
+  auto_fermion_enabled = ITensors.using_auto_fermion()
+  if auto_fermion_enabled
+    # If fermionic, bring indices into i',j',..,dag(j),dag(i)
+    # ordering with Out indices coming before In indices
+    # Resulting tensor acts like a normal matrix (no extra signs
+    # when taking powers A^n)
+    if all(j->dir(j)==Out, Linds) && all(j->dir(j)==In, Rinds)
+      ordered_inds = [Linds..., reverse(Rinds)...]
+    elseif all(j->dir(j)==Out, Rinds) && all(j->dir(j)==In, Linds)
+      ordered_inds = [Rinds..., reverse(Linds)...]
+    else
+      error(
+        "For fermionic exp, Linds and Rinds must have same directions within each set. Got dir.(Linds)=",
+        dir.(Linds),
+        ", dir.(Rinds)=",
+        dir.(Rinds),
+      )
+    end
+    A = permute(A, ordered_inds)
+    # A^n now sign free, ok to temporarily disable fermion system
+    ITensors.disable_auto_fermion()
+  end
+
   if isdiag(A)
-    return map_diag(f, A)
+    mapped_A = map_diag(f, A)
+  else
+    Ul, D, Ur = eigendecomp(A, Linds, Rinds; kws...)
+    mapped_A = Ul * map_diag(f, D) * Ur
   end
 
-  Ul, D, Ur = eigendecomp(A, inds...; ishermitian, kwargs...)
+  # <fermions>
+  if auto_fermion_enabled
+    # Ensure expA indices in "matrix" form before re-enabling fermion system
+    mapped_A = permute(mapped_A, ordered_inds)
+    ITensors.enable_auto_fermion()
+  end
 
-  return Ul * map_diag(f, D) * Ur
+  return mapped_A
 end
 
 # Analagous to `denseblocks`.

test/test_itensorsextensions.jl

@@ -4,17 +4,22 @@ using ITensors:
   ITensor,
   Index,
   QN,
+  apply,
   dag,
   delta,
   inds,
+  mapprime,
   noprime,
+  norm,
   op,
+  permute,
   prime,
   random_itensor,
   replaceind,
   replaceinds,
-  sim
-using ITensorNetworks.ITensorsExtensions: map_eigvals
+  sim,
+  swapprime
+using ITensorNetworks.ITensorsExtensions: eigendecomp, map_eigvals
 using StableRNGs: StableRNG
 using Test: @test, @testset
 @testset "ITensorsExtensions" begin
@@ -72,5 +77,54 @@ using Test: @test, @testset
     I = replaceind(inv_sqrtP * sqrtP, new_ind, i)
     @test I ≈ op("I", i)
   end
+
+  @testset "Fermionic eigendecomp" begin
+    s1 = Index([QN("Nf", 0, -1)=>2, QN("Nf", 1, -1)=>2], "Site,Fermion,n=1")
+    s2 = Index([QN("Nf", 0, -1)=>2, QN("Nf", 1, -1)=>2], "Site,Fermion,n=2")
+
+    # Make a random Hermitian matrix-like 4th order ITensor
+    T = random_itensor(s1', s2', dag(s2), dag(s1))
+    T = apply(T, swapprime(dag(T), 0=>1))
+    @test T ≈ swapprime(dag(T), 0=>1) # check Hermitian
+
+    Ul, D, Ur = eigendecomp(T, [s1', s2'], [dag(s1), dag(s2)]; ishermitian=true)
+
+    @test Ul*D*Ur ≈ T
+  end
+
+  @testset "Fermionic map eigvals tests" begin
+    s1 = Index([QN("Nf", 0, -1)=>2, QN("Nf", 1, -1)=>2], "Site,Fermion,n=1")
+    s2 = Index([QN("Nf", 0, -1)=>2, QN("Nf", 1, -1)=>2], "Site,Fermion,n=2")
+
+    # Make a random Hermitian matrix ITensor
+    M = random_itensor(s1', dag(s1))
+    #M = mapprime(prime(M)*swapprime(dag(M),0=>1),2=>1)
+    M = apply(M, swapprime(dag(M), 0=>1))
+
+    # Make a random Hermitian matrix-like 4th order ITensor
+    T = random_itensor(s1', s2', dag(s2), dag(s1))
+    T = apply(T, swapprime(dag(T), 0=>1))
+
+    # Matrix test
+    sqrtM = map_eigvals(sqrt, M, s1', dag(s1); ishermitian=true)
+    @test M ≈ apply(sqrtM, sqrtM)
+
+    ## Tensor test
+    sqrtT = map_eigvals(sqrt, T, [s1', s2'], [dag(s1), dag(s2)]; ishermitian=true)
+    @test T ≈ apply(sqrtT, sqrtT)
+
+    # Permute and test again
+    T = permute(T, dag(s2), s2', dag(s1), s1')
+    sqrtT = map_eigvals(sqrt, T, [s1', s2'], [dag(s1), dag(s2)]; ishermitian=true)
+    @test T ≈ apply(sqrtT, sqrtT)
+
+    ## Explicitly passing indices in different, valid orders
+    sqrtT = map_eigvals(sqrt, T, [s2', s1'], [dag(s2), dag(s1)]; ishermitian=true)
+    @test T ≈ apply(sqrtT, sqrtT)
+    sqrtT = map_eigvals(sqrt, T, [dag(s2), dag(s1)], [s2', s1'], ; ishermitian=true)
+    @test T ≈ apply(sqrtT, sqrtT)
+    sqrtT = map_eigvals(sqrt, T, [dag(s1), dag(s2)], [s1', s2'], ; ishermitian=true)
+    @test T ≈ apply(sqrtT, sqrtT)
+  end
 end
 end