add projection tests

Author: lkdvos

src/TensorKit.jl

@@ -137,7 +137,7 @@ using LinearAlgebra: norm, dot, normalize, normalize!, tr,
     adjoint, adjoint!, transpose, transpose!,
     lu, pinv, sylvester,
     eigen, eigen!, svd, svd!,
-    isposdef, isposdef!, ishermitian, rank, cond,
+    isposdef, isposdef!, rank, cond,
     Diagonal, Hermitian
 using MatrixAlgebraKit
 

src/factorizations/matrixalgebrakit.jl

@@ -46,7 +46,7 @@ for f! in (
 end
 
 # Handle these separately because single output instead of tuple
-for f! in (:qr_null!, :lq_null!)
+for f! in (:qr_null!, :lq_null!, :project_hermitian!, :project_antihermitian!, :project_isometric!)
     @eval function MAK.$f!(t::AbstractTensorMap, N, alg::AbstractAlgorithm)
         MAK.check_input($f!, t, N, alg)
 
@@ -444,14 +444,24 @@ end
 
 # Projections
 # -----------
-function MAK.check_input(::typeof(project_hermitian!), tsrc::AbstractTensorMap, tdst::AbstractTensorMap)
+function MAK.check_input(::typeof(project_hermitian!), tsrc::AbstractTensorMap, tdst::AbstractTensorMap, ::AbstractAlgorithm)
     domain(tsrc) == codomain(tsrc) || throw(ArgumentError("Hermitian projection requires square input tensor"))
     tsrc === tdst || @check_space(tdst, space(tsrc))
     return nothing
 end
 
-MAK.check_input(::typeof(project_antihermitian!), tsrc::AbstractTensorMap, tdst::AbstractTensorMap) =
-    MAK.check_input(project_hermitian!, tsrc, tdst)
+MAK.check_input(::typeof(project_antihermitian!), tsrc::AbstractTensorMap, tdst::AbstractTensorMap, alg::AbstractAlgorithm) =
+    MAK.check_input(project_hermitian!, tsrc, tdst, alg)
 
-MAK.initialize_output(::typeof(project_hermitian!), tsrc::AbstractTensorMap) = tsrc
-MAK.initialize_output(::typeof(project_antihermitian!), tsrc::AbstractTensorMap) = tsrc
+function MAK.check_input(::typeof(project_isometric!), t::AbstractTensorMap, W::AbstractTensorMap, alg::AbstractAlgorithm)
+    codomain(t) ≿ domain(t) || throw(ArgumentError("Isometric projection requires `codomain(t) ≿ domain(t)`"))
+    @check_space W space(t)
+    @check_scalar(W, t)
+
+    return nothing
+end
+
+
+MAK.initialize_output(::typeof(project_hermitian!), tsrc::AbstractTensorMap, ::AbstractAlgorithm) = tsrc
+MAK.initialize_output(::typeof(project_antihermitian!), tsrc::AbstractTensorMap, ::AbstractAlgorithm) = tsrc
+MAK.initialize_output(::typeof(project_isometric!), tsrc::AbstractTensorMap, ::AbstractAlgorithm) = similar(tsrc)

test/tensors/factorizations.jl

@@ -389,5 +389,62 @@ for V in spacelist
                 @test cond(t) ≈ λmax / λmin
             end
         end
+
+        @testset "Hermitian projections" begin
+            for T in eltypes,
+                    t in (
+                        rand(T, V1, V1), rand(T, W, W), rand(T, W, W)',
+                        DiagonalTensorMap(rand(T, reduceddim(V1)), V1),
+                    )
+                normalize!(t)
+                noisefactor = eps(real(T))^(3 / 4)
+
+                th = (t + t') / 2
+                ta = (t - t') / 2
+                tc = copy(t)
+
+                th′ = @constinferred project_hermitian(t)
+                @test ishermitian(th′)
+                @test th′ ≈ th
+                @test t == tc
+                th_approx = th + noisefactor * ta
+                @test !ishermitian(th_approx) || (T <: Real && t isa DiagonalTensorMap)
+                @test ishermitian(th_approx; atol = 10 * noisefactor)
+
+                ta′ = project_antihermitian(t)
+                @test isantihermitian(ta′)
+                @test ta′ ≈ ta
+                @test t == tc
+                ta_approx = ta + noisefactor * th
+                @test !isantihermitian(ta_approx)
+                @test isantihermitian(ta_approx; atol = 10 * noisefactor) || (T <: Real && t isa DiagonalTensorMap)
+            end
+        end
+
+        @testset "Isometric projections" begin
+            for T in eltypes,
+                    t in (
+                        randn(T, W, W), randn(T, W, W)',
+                        randn(T, W, V1), randn(T, V1, W)',
+                    )
+                t2 = project_isometric(t)
+                @test isisometric(t2)
+                t3 = project_isometric(t2)
+                @test t3 ≈ t2 # stability of the projection
+                @test t2 * (t2' * t) ≈ t
+
+                tc = similar(t)
+                t3 = @constinferred project_isometric!(copy!(tc, t), t2)
+                @test t3 === t2
+                @test isisometric(t2)
+
+                # test that t2 is closer to A then any other isometry
+                for k in 1:10
+                    δt = randn!(similar(t))
+                    t3 = project_isometric(t + δt / 100)
+                    @test norm(t - t3) > norm(t - t2)
+                end
+            end
+        end
     end
 end