From 79bdab00719c75a359d3a1ec444db13fda7fcf62 Mon Sep 17 00:00:00 2001 From: Miles Date: Fri, 1 Aug 2025 12:02:39 -0400 Subject: [PATCH 1/6] Fix map_eigvals for fermionic case and test --- src/lib/ITensorsExtensions/src/itensor.jl | 43 +++++++++++++++-- test/test_itensorsextensions.jl | 58 ++++++++++++++++++++++- 2 files changed, 94 insertions(+), 7 deletions(-) diff --git a/src/lib/ITensorsExtensions/src/itensor.jl b/src/lib/ITensorsExtensions/src/itensor.jl index 8c834830..32d0a5f5 100644 --- a/src/lib/ITensorsExtensions/src/itensor.jl +++ b/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...) + + # + 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...) + # + 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`. diff --git a/test/test_itensorsextensions.jl b/test/test_itensorsextensions.jl index cb4b5ed0..fb80fe85 100644 --- a/test/test_itensorsextensions.jl +++ b/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 From 6fa2bc6bcb08d615e56d1861d8d04bf68de3b347 Mon Sep 17 00:00:00 2001 From: Miles Date: Sat, 2 Aug 2025 12:35:27 -0400 Subject: [PATCH 2/6] Update version number --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index c3d97358..f4a4b689 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ITensorNetworks" uuid = "2919e153-833c-4bdc-8836-1ea460a35fc7" authors = ["Matthew Fishman , Joseph Tindall and contributors"] -version = "0.13.12" +version = "0.13.13" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" From 9cb594d5efa61f71df7ae4087cf72fa5ecfb1d11 Mon Sep 17 00:00:00 2001 From: Miles Date: Wed, 6 Aug 2025 13:48:49 -0400 Subject: [PATCH 3/6] Update map_eigvals to use newer pattern and make more robust --- src/lib/ITensorsExtensions/src/itensor.jl | 15 ++++++++++----- test/test_itensorsextensions.jl | 12 ++++++++++++ 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/src/lib/ITensorsExtensions/src/itensor.jl b/src/lib/ITensorsExtensions/src/itensor.jl index 32d0a5f5..e31cb157 100644 --- a/src/lib/ITensorsExtensions/src/itensor.jl +++ b/src/lib/ITensorsExtensions/src/itensor.jl @@ -5,6 +5,7 @@ using ITensors: Index, commonind, dag, + dir, hasqns, inds, isdiag, @@ -12,6 +13,7 @@ using ITensors: map_diag, noncommonind, noprime, + permute, replaceind, replaceinds, sim, @@ -63,17 +65,20 @@ function eigendecomp(A::ITensor, linds, rinds; ishermitian=false, kwargs...) end function map_eigvals(f::Function, A::ITensor, Linds, Rinds; kws...) + Linds = isa(Linds, Index) ? [Linds] : collect(Linds) + Rinds = isa(Rinds, Index) ? [Rinds] : collect(Rinds) # - auto_fermion_enabled = ITensors.using_auto_fermion() - if auto_fermion_enabled + fermionic_itensor = + ITensors.using_auto_fermion() && ITensors.has_fermionic_subspaces(inds(A)) + if fermionic_itensor # 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) + 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)==Out, Rinds) && all(j->dir(j)==In, Linds) + elseif all(j->dir(j)==ITensors.Out, Rinds) && all(j->dir(j)==ITensors.In, Linds) ordered_inds = [Rinds..., reverse(Linds)...] else error( @@ -96,7 +101,7 @@ function map_eigvals(f::Function, A::ITensor, Linds, Rinds; kws...) end # - if auto_fermion_enabled + if fermionic_itensor # Ensure expA indices in "matrix" form before re-enabling fermion system mapped_A = permute(mapped_A, ordered_inds) ITensors.enable_auto_fermion() diff --git a/test/test_itensorsextensions.jl b/test/test_itensorsextensions.jl index fb80fe85..88532f43 100644 --- a/test/test_itensorsextensions.jl +++ b/test/test_itensorsextensions.jl @@ -79,6 +79,7 @@ using Test: @test, @testset 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") @@ -90,9 +91,11 @@ using Test: @test, @testset 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") @@ -125,6 +128,15 @@ using Test: @test, @testset @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 From f8db6d93144e19eff6ed0a26229f67753e787fbf Mon Sep 17 00:00:00 2001 From: Miles Date: Thu, 7 Aug 2025 00:13:52 -0400 Subject: [PATCH 4/6] Add Linds, Rinds test to map_eigvals and factor out make_bosonic code --- src/lib/ITensorsExtensions/src/itensor.jl | 51 +++++++++++++---------- test/test_itensorsextensions.jl | 10 ++--- 2 files changed, 35 insertions(+), 26 deletions(-) diff --git a/src/lib/ITensorsExtensions/src/itensor.jl b/src/lib/ITensorsExtensions/src/itensor.jl index e31cb157..a89154af 100644 --- a/src/lib/ITensorsExtensions/src/itensor.jl +++ b/src/lib/ITensorsExtensions/src/itensor.jl @@ -64,32 +64,41 @@ function eigendecomp(A::ITensor, linds, rinds; ishermitian=false, kwargs...) return Ul, D, dag(U) end +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) + 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 check_input(::typeof(map_eigvals), f, A, Linds, Rinds) + all(x -> x isa Index, Linds) || error("Left indices must be a collection of Index") + all(x -> x isa Index, Rinds) || error("Right indices must be a collection of Index") +end + function map_eigvals(f::Function, A::ITensor, Linds, Rinds; kws...) - Linds = isa(Linds, Index) ? [Linds] : collect(Linds) - Rinds = isa(Rinds, Index) ? [Rinds] : collect(Rinds) + check_input(map_eigvals, f, A, Linds, Rinds) # fermionic_itensor = ITensors.using_auto_fermion() && ITensors.has_fermionic_subspaces(inds(A)) if fermionic_itensor - # 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)==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 - A = permute(A, ordered_inds) - # A^n now sign free, ok to temporarily disable fermion system + A = make_bosonic(A::ITensor, Linds, Rinds) + ordered_inds = inds(A) ITensors.disable_auto_fermion() end @@ -102,7 +111,7 @@ function map_eigvals(f::Function, A::ITensor, Linds, Rinds; kws...) # if fermionic_itensor - # Ensure expA indices in "matrix" form before re-enabling fermion system + # Ensure indices in "matrix" form before re-enabling fermion system mapped_A = permute(mapped_A, ordered_inds) ITensors.enable_auto_fermion() end diff --git a/test/test_itensorsextensions.jl b/test/test_itensorsextensions.jl index 88532f43..39ab8005 100644 --- a/test/test_itensorsextensions.jl +++ b/test/test_itensorsextensions.jl @@ -60,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) @@ -109,7 +109,7 @@ using Test: @test, @testset T = apply(T, swapprime(dag(T), 0=>1)) # Matrix test - sqrtM = map_eigvals(sqrt, M, s1', dag(s1); ishermitian=true) + sqrtM = map_eigvals(sqrt, M, [s1'], [dag(s1)]; ishermitian=true) @test M ≈ apply(sqrtM, sqrtM) ## Tensor test @@ -133,7 +133,7 @@ using Test: @test, @testset 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) + sqrtT = map_eigvals(sqrt, T, [b'], [dag(b)]; ishermitian=true) @test T ≈ apply(sqrtT, sqrtT) ITensors.disable_auto_fermion() From a42ace50763b368a0ecde21b7a9f19999ae2ff8f Mon Sep 17 00:00:00 2001 From: Miles Date: Thu, 7 Aug 2025 09:58:09 -0400 Subject: [PATCH 5/6] Move function to fix compilation error --- src/lib/ITensorsExtensions/src/itensor.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/lib/ITensorsExtensions/src/itensor.jl b/src/lib/ITensorsExtensions/src/itensor.jl index a89154af..46df13c5 100644 --- a/src/lib/ITensorsExtensions/src/itensor.jl +++ b/src/lib/ITensorsExtensions/src/itensor.jl @@ -85,11 +85,6 @@ function make_bosonic(A::ITensor, Linds, Rinds) return permute(A, ordered_inds) end -function check_input(::typeof(map_eigvals), f, A, Linds, Rinds) - all(x -> x isa Index, Linds) || error("Left indices must be a collection of Index") - all(x -> x isa Index, Rinds) || error("Right indices must be a collection of Index") -end - function map_eigvals(f::Function, A::ITensor, Linds, Rinds; kws...) check_input(map_eigvals, f, A, Linds, Rinds) @@ -119,6 +114,11 @@ function map_eigvals(f::Function, A::ITensor, Linds, Rinds; kws...) return mapped_A end +function check_input(::typeof(map_eigvals), f, A, Linds, Rinds) + all(x -> x isa Index, Linds) || error("Left indices must be a collection of Index") + all(x -> x isa Index, Rinds) || error("Right indices must be a collection of Index") +end + # Analagous to `denseblocks`. # Extract the diagonal entries into a diagonal tensor. function diagblocks(D::Tensor) From e99019f97344e41b8c713fae5393ee2d4ccf0cc6 Mon Sep 17 00:00:00 2001 From: Miles Date: Thu, 7 Aug 2025 11:49:34 -0400 Subject: [PATCH 6/6] Call indices to normalize Linds, Rinds inputs --- src/lib/ITensorsExtensions/src/itensor.jl | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/lib/ITensorsExtensions/src/itensor.jl b/src/lib/ITensorsExtensions/src/itensor.jl index 46df13c5..cc03865f 100644 --- a/src/lib/ITensorsExtensions/src/itensor.jl +++ b/src/lib/ITensorsExtensions/src/itensor.jl @@ -7,6 +7,7 @@ using ITensors: dag, dir, hasqns, + indices, inds, isdiag, itensor, @@ -69,6 +70,8 @@ function make_bosonic(A::ITensor, Linds, Rinds) # 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) @@ -86,8 +89,6 @@ function make_bosonic(A::ITensor, Linds, Rinds) end function map_eigvals(f::Function, A::ITensor, Linds, Rinds; kws...) - check_input(map_eigvals, f, A, Linds, Rinds) - # fermionic_itensor = ITensors.using_auto_fermion() && ITensors.has_fermionic_subspaces(inds(A)) @@ -114,11 +115,6 @@ function map_eigvals(f::Function, A::ITensor, Linds, Rinds; kws...) return mapped_A end -function check_input(::typeof(map_eigvals), f, A, Linds, Rinds) - all(x -> x isa Index, Linds) || error("Left indices must be a collection of Index") - all(x -> x isa Index, Rinds) || error("Right indices must be a collection of Index") -end - # Analagous to `denseblocks`. # Extract the diagonal entries into a diagonal tensor. function diagblocks(D::Tensor)