Fix map_eigvals for fermionic case and test (#250)

Co-authored-by: Matt Fishman <[email protected]>

Author: emstoudenmire

src/lib/ITensorsExtensions/src/itensor.jl

@@ -1,16 +1,20 @@
 using LinearAlgebra: LinearAlgebra, eigen, pinv
 using ITensors:
+  ITensors,
   ITensor,
   Index,
   commonind,
   dag,
+  dir,
   hasqns,
+  indices,
   inds,
   isdiag,
   itensor,
   map_diag,
   noncommonind,
   noprime,
+  permute,
   replaceind,
   replaceinds,
   sim,
@@ -58,18 +62,57 @@ 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 make_bosonic(A::ITensor, Linds, Rinds)
+  # 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)
+  Linds = indices(Linds)
+  Rinds = indices(Rinds)
+  if all(j->dir(j)==ITensors.Out, Linds) && all(j->dir(j)==ITensors.In, Rinds)
+    ordered_inds = [Linds..., reverse(Rinds)...]
+  elseif all(j->dir(j)==ITensors.Out, Rinds) && all(j->dir(j)==ITensors.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
+  # permuted A^n will be sign free, ok to temporarily disable fermion system
+  return permute(A, ordered_inds)
+end
+
+function map_eigvals(f::Function, A::ITensor, Linds, Rinds; kws...)
+  # <fermions>
+  fermionic_itensor =
+    ITensors.using_auto_fermion() && ITensors.has_fermionic_subspaces(inds(A))
+  if fermionic_itensor
+    A = make_bosonic(A::ITensor, Linds, Rinds)
+    ordered_inds = inds(A)
+    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 fermionic_itensor
+    # Ensure 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
@@ -55,9 +60,9 @@ using Test: @test, @testset
     rng = StableRNG(1234)
     A = random_itensor(rng, elt, i, j)
     P = A * prime(dag(A), i)
-    sqrtP = map_eigvals(sqrt, P, i, i'; ishermitian=true)
-    inv_P = dag(map_eigvals(inv, P, i, i'; ishermitian=true))
-    inv_sqrtP = dag(map_eigvals(inv ∘ sqrt, P, i, i'; ishermitian=true))
+    sqrtP = map_eigvals(sqrt, P, [i], [i']; ishermitian=true)
+    inv_P = dag(map_eigvals(inv, P, [i], [i']; ishermitian=true))
+    inv_sqrtP = dag(map_eigvals(inv ∘ sqrt, P, [i], [i']; ishermitian=true))
 
     new_ind = noprime(sim(i'))
     sqrtPdag = replaceind(dag(sqrtP), i', new_ind)
@@ -72,5 +77,66 @@ using Test: @test, @testset
     I = replaceind(inv_sqrtP * sqrtP, new_ind, i)
     @test I ≈ op("I", i)
   end
+
+  @testset "Fermionic eigendecomp" begin
+    ITensors.enable_auto_fermion()
+    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
+    ITensors.disable_auto_fermion()
+  end
+
+  @testset "Fermionic map eigvals tests" begin
+    ITensors.enable_auto_fermion()
+    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)
+
+    # Test bosonic index case while fermion system is enabled
+    b = Index([QN("Nb", 0)=>2, QN("Nb", 1)=>2])
+    T = random_itensor(b', dag(b))
+    T = apply(T, swapprime(dag(T), 0=>1))
+    sqrtT = map_eigvals(sqrt, T, [b'], [dag(b)]; ishermitian=true)
+    @test T ≈ apply(sqrtT, sqrtT)
+
+    ITensors.disable_auto_fermion()
+  end
 end
 end