Fix and refactor truncation tests

Author: lkdvos

src/tensors/factorizations/truncation.jl

@@ -36,6 +36,7 @@ truncspace(space::ElementarySpace) = TruncationSpace(space)
 function truncate!(::typeof(svd_trunc!),
                    (U, S, Vᴴ)::Tuple{AbstractTensorMap,AbstractTensorMap,AbstractTensorMap},
                    strategy::TruncationStrategy)
+    strategy == notrunc() && return (U, S, Vᴴ)
     ind = findtruncated_sorted(diagview(S), strategy)
     V_truncated = spacetype(S)(c => length(I) for (c, I) in ind)
 
@@ -66,6 +67,7 @@ end
 function truncate!(::typeof(left_null!),
                    (U, S)::Tuple{AbstractTensorMap,AbstractTensorMap},
                    strategy::MatrixAlgebraKit.TruncationStrategy)
+    strategy == notrunc() && return (U, S)
     extended_S = SectorDict(c => vcat(diagview(b),
                                       zeros(eltype(b), max(0, size(b, 2) - size(b, 1))))
                             for (c, b) in blocks(S))
@@ -82,6 +84,7 @@ for f! in (:eig_trunc!, :eigh_trunc!)
     @eval function truncate!(::typeof($f!),
                              (D, V)::Tuple{AbstractTensorMap,AbstractTensorMap},
                              strategy::TruncationStrategy)
+        strategy == notrunc() && return (D, V)
         ind = findtruncated(diagview(D), strategy)
         V_truncated = spacetype(D)(c => length(I) for (c, I) in ind)
 
@@ -136,10 +139,14 @@ function _findnexttruncvalue(S, truncdim::SectorDict{I,Int}; by=identity,
 end
 
 # implementations
+function findtruncated_sorted(S::SectorDict, strategy::TruncationStrategy)
+    return findtruncated(S, strategy)
+end
+
 function findtruncated_sorted(S::SectorDict, strategy::TruncationKeepAbove)
     atol = rtol_to_atol(S, strategy.p, strategy.atol, strategy.rtol)
     findtrunc = Base.Fix2(findtruncated_sorted, truncbelow(atol))
-    return SectorDict(c => findtrunc(d) for (c, d) in Sd)
+    return SectorDict(c => findtrunc(d) for (c, d) in S)
 end
 function findtruncated(S::SectorDict, strategy::TruncationKeepAbove)
     atol = rtol_to_atol(S, strategy.p, strategy.atol, strategy.rtol)
@@ -195,13 +202,10 @@ function findtruncated(Sd::SectorDict, strategy::TruncationKeepSorted)
         _, cmin = next
         truncdim[cmin] -= 1
         totaldim -= dim(cmin)
-        if totaldim < strategy.howmany
-            # truncdim[cmin] += 1
-            break
-        end
         if truncdim[cmin] == 0
             delete!(truncdim, cmin)
         end
+        totaldim <= strategy.howmany && break
     end
     return SectorDict(c => permutations[c][Base.OneTo(d)] for (c, d) in truncdim)
 end

test/factorizations.jl

@@ -204,6 +204,54 @@ for V in spacelist
             end
         end
 
+        @testset "truncated SVD" begin
+            for T in eltypes,
+                t in (randn(T, W, W), randn(T, W, W)',
+                      randn(T, W, V1), randn(T, V1, W),
+                      randn(T, W, V1)', randn(T, V1, W)',
+                      DiagonalTensorMap(randn(T, reduceddim(V1)), V1))
+
+                @constinferred normalize!(t)
+
+                U, S, Vᴴ = @constinferred svd_trunc(t; trunc=notrunc())
+                @test U * S * Vᴴ ≈ t
+                @test isisometry(U)
+                @test isisometry(Vᴴ; side=:right)
+
+                trunc = truncrank(dim(domain(S)) ÷ 2)
+                U1, S1, Vᴴ1 = @constinferred svd_trunc(t; trunc)
+                @test t * Vᴴ1' ≈ U1 * S1
+                @test isisometry(U1)
+                @test isisometry(Vᴴ1; side=:right)
+                @test dim(domain(S1)) <= trunc.howmany
+
+                λ = minimum(minimum, values(LinearAlgebra.diag(S1)))
+                trunc = trunctol(λ - 10eps(λ))
+                U2, S2, Vᴴ2 = @constinferred svd_trunc(t; trunc)
+                @test t * Vᴴ2' ≈ U2 * S2
+                @test isisometry(U2)
+                @test isisometry(Vᴴ2; side=:right)
+                @test minimum(minimum, values(LinearAlgebra.diag(S1))) >= λ
+                @test U2 ≈ U1
+                @test S2 ≈ S1
+                @test Vᴴ2 ≈ Vᴴ1
+
+                trunc = truncspace(space(S2, 1))
+                U3, S3, Vᴴ3 = @constinferred svd_trunc(t; trunc)
+                @test t * Vᴴ3' ≈ U3 * S3
+                @test isisometry(U3)
+                @test isisometry(Vᴴ3; side=:right)
+                @test space(S3, 1) ≾ space(S2, 1)
+
+                trunc = truncerr(0.5)
+                U4, S4, Vᴴ4 = @constinferred svd_trunc(t; trunc)
+                @test t * Vᴴ4' ≈ U4 * S4
+                @test isisometry(U4)
+                @test isisometry(Vᴴ4; side=:right)
+                @test norm(t - U4 * S4 * Vᴴ4) <= 0.5
+            end
+        end
+
         @testset "Eigenvalue decomposition" begin
             for T in eltypes,
                 t in

test/tensors.jl

@@ -432,50 +432,6 @@ for V in spacelist
             @test LinearAlgebra.isdiag(D)
             @test LinearAlgebra.diag(D) == d
         end
-        @timedtestset "Tensor truncation" begin
-            for T in (Float32, ComplexF64)
-                for p in (1, 2, 3, Inf)
-                    # Test both a normal tensor and an adjoint one.
-                    ts = (randn(T, V1 ⊗ V2 ⊗ V3, V4 ⊗ V5),
-                          randn(T, V4 ⊗ V5, V1 ⊗ V2 ⊗ V3)')
-                    for t in ts
-                        U₀, S₀, V₀, = tsvd(t)
-                        t = rmul!(t, 1 / norm(S₀, p))
-                        U, S, V = @constinferred tsvd(t; trunc=truncerr(5e-1, p))
-                        ϵ = TensorKit._norm(LinearAlgebra.svdvals(U * S * V - t), p,
-                                            zero(scalartype(S)))
-                        p == 2 && @test ϵ < 5e-1
-                        # @show p, ϵ
-                        # @show domain(S)
-                        # @test min(space(S,1), space(S₀,1)) != space(S₀,1)
-                        U′, S′, V′ = tsvd(t; trunc=truncerr(ϵ + 10eps(ϵ), p))
-                        ϵ′ = TensorKit._norm(LinearAlgebra.svdvals(U′ * S′ * V′ - t), p,
-                                             zero(scalartype(S)))
-
-                        @test (U, S, V, ϵ) == (U′, S′, V′, ϵ′)
-                        U′, S′, V′ = tsvd(t; trunc=truncdim(ceil(Int, dim(domain(S)))))
-                        ϵ′ = TensorKit._norm(LinearAlgebra.svdvals(U′ * S′ * V′ - t), p,
-                                             zero(scalartype(S)))
-                        @test (U, S, V, ϵ) == (U′, S′, V′, ϵ′)
-                        U′, S′, V′ = tsvd(t; trunc=truncspace(space(S, 1)))
-                        ϵ′ = TensorKit._norm(LinearAlgebra.svdvals(U′ * S′ * V′ - t), p,
-                                             zero(scalartype(S)))
-                        @test (U, S, V, ϵ) == (U′, S′, V′, ϵ′)
-                        # results with truncationcutoff cannot be compared because they don't take degeneracy into account, and thus truncate differently
-                        U, S, V = tsvd(t; trunc=truncbelow(1 / dim(domain(S₀))))
-                        ϵ = TensorKit._norm(LinearAlgebra.svdvals(U * S * V - t), p,
-                                            zero(scalartype(S)))
-                        # @show p, ϵ
-                        # @show domain(S)
-                        # @test min(space(S,1), space(S₀,1)) != space(S₀,1)
-                        U′, S′, V′ = tsvd(t; trunc=truncspace(space(S, 1)))
-                        ϵ′ = TensorKit._norm(LinearAlgebra.svdvals(U′ * S′ * V′ - t), p,
-                                             zero(scalartype(S)))
-                        @test (U, S, V, ϵ) == (U′, S′, V′, ϵ′)
-                    end
-                end
-            end
-        end
         if BraidingStyle(I) isa Bosonic && hasfusiontensor(I)
             @timedtestset "Tensor functions" begin
                 W = V1 ⊗ V2