From add45859d88058e4f7d2f0c068b1907be13e5b3f Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 27 Oct 2025 10:48:31 -0400 Subject: [PATCH 01/41] bump MatrixAlgebraKit compat to 0.6 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 7d754d817..137aa050c 100644 --- a/Project.toml +++ b/Project.toml @@ -35,7 +35,7 @@ Combinatorics = "1" FiniteDifferences = "0.12" LRUCache = "1.0.2" LinearAlgebra = "1" -MatrixAlgebraKit = "0.5.0" +MatrixAlgebraKit = "0.6.0" OhMyThreads = "0.8.0" PackageExtensionCompat = "1" Printf = "1" From b102198406da9d8912cb1d0f8f85e442a49e30f8 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 27 Oct 2025 10:56:54 -0400 Subject: [PATCH 02/41] update `ishermitian` implementations --- src/factorizations/factorizations.jl | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index 29d2c16a6..1e0a37078 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -61,14 +61,17 @@ LinearAlgebra.svdvals!(t::AbstractTensorMap) = diagview(svd_vals!(t)) #--------------------------------------------------# # Checks for hermiticity and positive definiteness # #--------------------------------------------------# -function LinearAlgebra.ishermitian(t::AbstractTensorMap) - domain(t) == codomain(t) || return false - InnerProductStyle(t) === EuclideanInnerProduct() || return false # hermiticity only defined for euclidean - for (c, b) in blocks(t) - ishermitian(b) || return false - end - return true +function MAK.ishermitian(t::AbstractTensorMap; kwargs...) + return InnerProductStyle(t) === EuclideanInnerProduct() && + domain(t) == codomain(t) && + all(cb -> MatrixAlgebraKit.ishermitian(cb[2]; kwargs...), blocks(t)) +end +function MAK.isantihermitian(t::AbstractTensorMap; kwargs...) + return InnerProductStyle(t) === EuclideanInnerProduct() && + domain(t) == codomain(t) && + all(cb -> MatrixAlgebraKit.isantihermitian(cb[2]; kwargs...), blocks(t)) end +LinearAlgebra.ishermitian(t::AbstractTensorMap) = MAK.ishermitian(t) function LinearAlgebra.isposdef(t::AbstractTensorMap) return isposdef!(copy_oftype(t, factorisation_scalartype(isposdef, t))) From f7fae127d286168fdcae7449ebf0cae1a17d2903 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 27 Oct 2025 10:57:08 -0400 Subject: [PATCH 03/41] update `isisometric` implementations --- src/factorizations/factorizations.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index 1e0a37078..42705f432 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -87,14 +87,14 @@ function LinearAlgebra.isposdef!(t::AbstractTensorMap) end # TODO: tolerances are per-block, not global or weighted - does that matter? -function MatrixAlgebraKit.is_left_isometry(t::AbstractTensorMap; kwargs...) +function MAK.is_left_isometric(t::AbstractTensorMap; kwargs...) domain(t) ≾ codomain(t) || return false - f((c, b)) = MatrixAlgebraKit.is_left_isometry(b; kwargs...) + f((c, b)) = MAK.is_left_isometric(b; kwargs...) return all(f, blocks(t)) end -function MatrixAlgebraKit.is_right_isometry(t::AbstractTensorMap; kwargs...) +function MAK.is_right_isometric(t::AbstractTensorMap; kwargs...) domain(t) ≿ codomain(t) || return false - f((c, b)) = MatrixAlgebraKit.is_right_isometry(b; kwargs...) + f((c, b)) = MAK.is_right_isometric(b; kwargs...) return all(f, blocks(t)) end From f2067864ed5b53809aa783a255c6ad18f6e4b733 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 27 Oct 2025 11:13:26 -0400 Subject: [PATCH 04/41] add projections --- src/factorizations/matrixalgebrakit.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/factorizations/matrixalgebrakit.jl b/src/factorizations/matrixalgebrakit.jl index fc5f84502..5a51d3aff 100644 --- a/src/factorizations/matrixalgebrakit.jl +++ b/src/factorizations/matrixalgebrakit.jl @@ -5,6 +5,7 @@ for f in :svd_compact, :svd_full, :svd_trunc, :svd_vals, :qr_compact, :qr_full, :qr_null, :lq_compact, :lq_full, :lq_null, :eig_full, :eig_trunc, :eig_vals, :eigh_full, :eigh_trunc, :eigh_vals, :left_polar, :right_polar, + :project_hermitian, :project_antihermitian, :project_isometric, ] f! = Symbol(f, :!) @eval function MAK.default_algorithm(::typeof($f!), ::Type{T}; kwargs...) where {T <: AbstractTensorMap} @@ -655,3 +656,17 @@ for (f!, f_svd!) in zip((:left_null!, :right_null!), (:left_null_svd!, :right_nu return $f!(t, N; alg_svd = alg) end end + +# Projections +# ----------- +function MAK.check_input(::typeof(project_hermitian!), tsrc::AbstractTensorMap, tdst::AbstractTensorMap) + 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.initialize_output(::typeof(project_hermitian!), tsrc::AbstractTensorMap) = tsrc +MAK.initialize_output(::typeof(project_antihermitian!), tsrc::AbstractTensorMap) = tsrc From b7db0604c1c7398d224e53ff0f586f3964462dc5 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Fri, 14 Nov 2025 08:06:28 -0500 Subject: [PATCH 05/41] Incremental fixes --- src/TensorKit.jl | 2 +- src/factorizations/adjoint.jl | 8 +-- src/factorizations/factorizations.jl | 7 +-- src/factorizations/matrixalgebrakit.jl | 79 +------------------------- test/tensors/factorizations.jl | 66 ++++++++++----------- 5 files changed, 42 insertions(+), 120 deletions(-) diff --git a/src/TensorKit.jl b/src/TensorKit.jl index 74d23c8ce..85bc9a051 100644 --- a/src/TensorKit.jl +++ b/src/TensorKit.jl @@ -82,7 +82,7 @@ export left_orth, right_orth, left_null, right_null, eigh_full!, eigh_full, eigh_trunc!, eigh_trunc, eig_full!, eig_full, eig_trunc!, eig_trunc, eigh_vals!, eigh_vals, eig_vals!, eig_vals, - isposdef, isposdef!, ishermitian, isisometry, isunitary, sylvester, rank, cond + isposdef, isposdef!, ishermitian, isisometric, isunitary, sylvester, rank, cond export braid, braid!, permute, permute!, transpose, transpose!, twist, twist!, repartition, repartition! diff --git a/src/factorizations/adjoint.jl b/src/factorizations/adjoint.jl index 30a68c138..ef7e5e265 100644 --- a/src/factorizations/adjoint.jl +++ b/src/factorizations/adjoint.jl @@ -29,11 +29,11 @@ function MAK.right_null!(t::AdjointTensorMap, N, alg::AbstractAlgorithm) return N end -function MAK.is_left_isometry(t::AdjointTensorMap; kwargs...) - return is_right_isometry(adjoint(t); kwargs...) +function MAK.is_left_isometric(t::AdjointTensorMap; kwargs...) + return is_right_isometric(adjoint(t); kwargs...) end -function MAK.is_right_isometry(t::AdjointTensorMap; kwargs...) - return is_left_isometry(adjoint(t); kwargs...) +function MAK.is_right_isometric(t::AdjointTensorMap; kwargs...) + return is_left_isometric(adjoint(t); kwargs...) end # 2-arg functions diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index 42705f432..6961d55aa 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -18,8 +18,7 @@ import MatrixAlgebraKit as MAK using MatrixAlgebraKit: AbstractAlgorithm, TruncatedAlgorithm, DiagonalAlgorithm using MatrixAlgebraKit: TruncationStrategy, NoTruncation, TruncationByValue, TruncationByError, TruncationIntersection, TruncationByFilter, TruncationByOrder -using MatrixAlgebraKit: left_orth_polar!, right_orth_polar!, left_orth_svd!, - right_orth_svd!, left_null_svd!, right_null_svd!, diagview +using MatrixAlgebraKit: diagview, isisometric include("utility.jl") include("matrixalgebrakit.jl") @@ -30,9 +29,9 @@ include("pullbacks.jl") TensorKit.one!(A::AbstractMatrix) = MatrixAlgebraKit.one!(A) -function MatrixAlgebraKit.isisometry(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple) +function MatrixAlgebraKit.isisometric(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple) t = permute(t, (p₁, p₂); copy = false) - return isisometry(t) + return isisometric(t) end #------------------------------# diff --git a/src/factorizations/matrixalgebrakit.jl b/src/factorizations/matrixalgebrakit.jl index 5a51d3aff..974b577bd 100644 --- a/src/factorizations/matrixalgebrakit.jl +++ b/src/factorizations/matrixalgebrakit.jl @@ -26,8 +26,7 @@ end for f! in ( :qr_compact!, :qr_full!, :lq_compact!, :lq_full!, :eig_full!, :eigh_full!, :svd_compact!, :svd_full!, - :left_polar!, :left_orth_polar!, :right_polar!, :right_orth_polar!, - :left_orth!, :right_orth!, + :left_polar!, :right_polar!, :left_orth!, :right_orth!, ) @eval function MAK.$f!(t::AbstractTensorMap, F, alg::AbstractAlgorithm) MAK.check_input($f!, t, F, alg) @@ -417,26 +416,6 @@ function MAK.check_input(::typeof(left_polar!), t::AbstractTensorMap, WP, ::Abst return nothing end -function MAK.check_input(::typeof(left_orth_polar!), t::AbstractTensorMap, WP, ::AbstractAlgorithm) - codomain(t) ≿ domain(t) || - throw(ArgumentError("Polar decomposition requires `codomain(t) ≿ domain(t)`")) - - W, P = WP - @assert W isa AbstractTensorMap - @assert P isa AbstractTensorMap - - # scalartype checks - @check_scalar W t - @check_scalar P t - - # space checks - VW = fuse(domain(t)) - @check_space(W, codomain(t) ← VW) - @check_space(P, VW ← domain(t)) - - return nothing -end - function MAK.initialize_output(::typeof(left_polar!), t::AbstractTensorMap, ::AbstractAlgorithm) W = similar(t, space(t)) P = similar(t, domain(t) ← domain(t)) @@ -462,26 +441,6 @@ function MAK.check_input(::typeof(right_polar!), t::AbstractTensorMap, PWᴴ, :: return nothing end -function MAK.check_input(::typeof(right_orth_polar!), t::AbstractTensorMap, PWᴴ, ::AbstractAlgorithm) - codomain(t) ≾ domain(t) || - throw(ArgumentError("Polar decomposition requires `domain(t) ≿ codomain(t)`")) - - P, Wᴴ = PWᴴ - @assert P isa AbstractTensorMap - @assert Wᴴ isa AbstractTensorMap - - # scalartype checks - @check_scalar P t - @check_scalar Wᴴ t - - # space checks - VW = fuse(codomain(t)) - @check_space(P, codomain(t) ← VW) - @check_space(Wᴴ, VW ← domain(t)) - - return nothing -end - function MAK.initialize_output(::typeof(right_polar!), t::AbstractTensorMap, ::AbstractAlgorithm) P = similar(t, codomain(t) ← codomain(t)) Wᴴ = similar(t, space(t)) @@ -582,36 +541,6 @@ function MAK.right_orth!( end end -function MAK.left_orth_polar!(t::AbstractTensorMap; alg = nothing, kwargs...) - alg′ = MAK.select_algorithm(left_polar!, t, alg; kwargs...) - VC = MAK.initialize_output(left_orth!, t) - return left_orth_polar!(t, VC, alg′) -end -function MAK.left_orth_polar!(t::AbstractTensorMap, VC, alg) - alg′ = MAK.select_algorithm(left_polar!, t, alg) - return left_orth_polar!(t, VC, alg′) -end -function MAK.right_orth_polar!(t::AbstractTensorMap; alg = nothing, kwargs...) - alg′ = MAK.select_algorithm(right_polar!, t, alg; kwargs...) - CVᴴ = MAK.initialize_output(right_orth!, t) - return right_orth_polar!(t, CVᴴ, alg′) -end -function MAK.right_orth_polar!(t::AbstractTensorMap, CVᴴ, alg) - alg′ = MAK.select_algorithm(right_polar!, t, alg) - return right_orth_polar!(t, CVᴴ, alg′) -end - -function MAK.left_orth_svd!(t::AbstractTensorMap; trunc = notrunc(), kwargs...) - U, S, Vᴴ = trunc == notrunc() ? svd_compact!(t; kwargs...) : - svd_trunc!(t; trunc, kwargs...) - return U, lmul!(S, Vᴴ) -end -function MAK.right_orth_svd!(t::AbstractTensorMap; trunc = notrunc(), kwargs...) - U, S, Vᴴ = trunc == notrunc() ? svd_compact!(t; kwargs...) : - svd_trunc!(t; trunc, kwargs...) - return rmul!(U, S), Vᴴ -end - # Nullspace # --------- function MAK.check_input(::typeof(left_null!), t::AbstractTensorMap, N, ::AbstractAlgorithm) @@ -651,12 +580,6 @@ function MAK.initialize_output(::typeof(right_null!), t::AbstractTensorMap) return N end -for (f!, f_svd!) in zip((:left_null!, :right_null!), (:left_null_svd!, :right_null_svd!)) - @eval function MAK.$f_svd!(t::AbstractTensorMap, N, alg, ::Nothing = nothing) - return $f!(t, N; alg_svd = alg) - end -end - # Projections # ----------- function MAK.check_input(::typeof(project_hermitian!), tsrc::AbstractTensorMap, tdst::AbstractTensorMap) diff --git a/test/tensors/factorizations.jl b/test/tensors/factorizations.jl index c5b4e8748..8c9bde170 100644 --- a/test/tensors/factorizations.jl +++ b/test/tensors/factorizations.jl @@ -47,18 +47,18 @@ for V in spacelist Q, R = @constinferred qr_compact(t) @test Q * R ≈ t - @test isisometry(Q) + @test isisometric(Q) Q, R = @constinferred left_orth(t; kind = :qr) @test Q * R ≈ t - @test isisometry(Q) + @test isisometric(Q) N = @constinferred qr_null(t) - @test isisometry(N) + @test isisometric(N) @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t)) N = @constinferred left_null(t; kind = :qr) - @test isisometry(N) + @test isisometric(N) @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t)) end @@ -73,12 +73,12 @@ for V in spacelist Q, R = @constinferred qr_compact(t) @test Q * R ≈ t - @test isisometry(Q) + @test isisometric(Q) @test dim(Q) == dim(R) == dim(t) Q, R = @constinferred left_orth(t; kind = :qr) @test Q * R ≈ t - @test isisometry(Q) + @test isisometric(Q) @test dim(Q) == dim(R) == dim(t) N = @constinferred qr_null(t) @@ -100,14 +100,14 @@ for V in spacelist L, Q = @constinferred lq_compact(t) @test L * Q ≈ t - @test isisometry(Q; side = :right) + @test isisometric(Q; side = :right) L, Q = @constinferred right_orth(t; kind = :lq) @test L * Q ≈ t - @test isisometry(Q; side = :right) + @test isisometric(Q; side = :right) Nᴴ = @constinferred lq_null(t) - @test isisometry(Nᴴ; side = :right) + @test isisometric(Nᴴ; side = :right) @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t)) end @@ -122,12 +122,12 @@ for V in spacelist L, Q = @constinferred lq_compact(t) @test L * Q ≈ t - @test isisometry(Q; side = :right) + @test isisometric(Q; side = :right) @test dim(Q) == dim(L) == dim(t) L, Q = @constinferred right_orth(t; kind = :lq) @test L * Q ≈ t - @test isisometry(Q; side = :right) + @test isisometric(Q; side = :right) @test dim(Q) == dim(L) == dim(t) Nᴴ = @constinferred lq_null(t) @@ -146,12 +146,12 @@ for V in spacelist @assert domain(t) ≾ codomain(t) w, p = @constinferred left_polar(t) @test w * p ≈ t - @test isisometry(w) + @test isisometric(w) @test isposdef(p) w, p = @constinferred left_orth(t; kind = :polar) @test w * p ≈ t - @test isisometry(w) + @test isisometric(w) end for T in eltypes, @@ -160,12 +160,12 @@ for V in spacelist @assert codomain(t) ≾ domain(t) p, wᴴ = @constinferred right_polar(t) @test p * wᴴ ≈ t - @test isisometry(wᴴ; side = :right) + @test isisometric(wᴴ; side = :right) @test isposdef(p) p, wᴴ = @constinferred right_orth(t; kind = :polar) @test p * wᴴ ≈ t - @test isisometry(wᴴ; side = :right) + @test isisometric(wᴴ; side = :right) end end @@ -185,9 +185,9 @@ for V in spacelist u, s, vᴴ = @constinferred svd_compact(t) @test u * s * vᴴ ≈ t - @test isisometry(u) + @test isisometric(u) @test isposdef(s) - @test isisometry(vᴴ; side = :right) + @test isisometric(vᴴ; side = :right) s′ = LinearAlgebra.diag(s) for (c, b) in LinearAlgebra.svdvals(t) @@ -196,14 +196,14 @@ for V in spacelist v, c = @constinferred left_orth(t; kind = :svd) @test v * c ≈ t - @test isisometry(v) + @test isisometric(v) N = @constinferred left_null(t; kind = :svd) - @test isisometry(N) + @test isisometric(N) @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t)) Nᴴ = @constinferred right_null(t; kind = :svd) - @test isisometry(Nᴴ; side = :right) + @test isisometric(Nᴴ; side = :right) @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t)) end @@ -233,22 +233,22 @@ for V in spacelist U, S, Vᴴ = @constinferred svd_trunc(t; trunc = notrunc()) @test U * S * Vᴴ ≈ t - @test isisometry(U) - @test isisometry(Vᴴ; side = :right) + @test isisometric(U) + @test isisometric(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 isisometric(U1) + @test isisometric(Vᴴ1; side = :right) @test dim(domain(S1)) <= trunc.howmany λ = minimum(minimum, values(LinearAlgebra.diag(S1))) trunc = trunctol(; atol = λ - 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 isisometric(U2) + @test isisometric(Vᴴ2; side = :right) @test minimum(minimum, values(LinearAlgebra.diag(S1))) >= λ @test U2 ≈ U1 @test S2 ≈ S1 @@ -257,22 +257,22 @@ for V in spacelist 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 isisometric(U3) + @test isisometric(Vᴴ3; side = :right) @test space(S3, 1) ≾ space(S2, 1) trunc = truncerror(; atol = 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 isisometric(U4) + @test isisometric(Vᴴ4; side = :right) @test norm(t - U4 * S4 * Vᴴ4) <= 0.5 trunc = truncrank(dim(domain(S)) ÷ 2) & trunctol(; atol = λ - 10eps(λ)) U5, S5, Vᴴ5 = @constinferred svd_trunc(t; trunc) @test t * Vᴴ5' ≈ U5 * S5 - @test isisometry(U5) - @test isisometry(Vᴴ5; side = :right) + @test isisometric(U5) + @test isisometric(Vᴴ5; side = :right) @test minimum(minimum, values(LinearAlgebra.diag(S5))) >= λ @test dim(domain(S5)) ≤ dim(domain(S)) ÷ 2 end @@ -304,7 +304,7 @@ for V in spacelist t2 = (t + t') D, V = eigen(t2) - @test isisometry(V) + @test isisometric(V) D̃, Ṽ = @constinferred eigh_full(t2) @test D ≈ D̃ @test V ≈ Ṽ From 5890ad4e08734ee3e8f1538093aa16743f81538d Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Fri, 14 Nov 2025 15:30:38 +0100 Subject: [PATCH 06/41] orthnull progress --- src/TensorKit.jl | 3 +- src/factorizations/factorizations.jl | 1 + src/factorizations/matrixalgebrakit.jl | 99 +++++++++++++++++--------- test/autodiff/ad.jl | 6 +- test/tensors/factorizations.jl | 20 +++--- 5 files changed, 84 insertions(+), 45 deletions(-) diff --git a/src/TensorKit.jl b/src/TensorKit.jl index 85bc9a051..7ebdb504d 100644 --- a/src/TensorKit.jl +++ b/src/TensorKit.jl @@ -82,7 +82,8 @@ export left_orth, right_orth, left_null, right_null, eigh_full!, eigh_full, eigh_trunc!, eigh_trunc, eig_full!, eig_full, eig_trunc!, eig_trunc, eigh_vals!, eigh_vals, eig_vals!, eig_vals, - isposdef, isposdef!, ishermitian, isisometric, isunitary, sylvester, rank, cond + isposdef, isposdef!, ishermitian, isisometric, isunitary, sylvester, rank, cond, + LeftOrthAlgorithm, RightOrthAlgorithm export braid, braid!, permute, permute!, transpose, transpose!, twist, twist!, repartition, repartition! diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index 6961d55aa..d618bca09 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -16,6 +16,7 @@ using TensorOperations: Index2Tuple using MatrixAlgebraKit import MatrixAlgebraKit as MAK using MatrixAlgebraKit: AbstractAlgorithm, TruncatedAlgorithm, DiagonalAlgorithm +using MatrixAlgebraKit: LeftOrthAlgorithm, RightOrthAlgorithm using MatrixAlgebraKit: TruncationStrategy, NoTruncation, TruncationByValue, TruncationByError, TruncationIntersection, TruncationByFilter, TruncationByOrder using MatrixAlgebraKit: diagview, isisometric diff --git a/src/factorizations/matrixalgebrakit.jl b/src/factorizations/matrixalgebrakit.jl index 974b577bd..73fd04a6b 100644 --- a/src/factorizations/matrixalgebrakit.jl +++ b/src/factorizations/matrixalgebrakit.jl @@ -26,7 +26,7 @@ end for f! in ( :qr_compact!, :qr_full!, :lq_compact!, :lq_full!, :eig_full!, :eigh_full!, :svd_compact!, :svd_full!, - :left_polar!, :right_polar!, :left_orth!, :right_orth!, + :left_polar!, :right_polar! ) @eval function MAK.$f!(t::AbstractTensorMap, F, alg::AbstractAlgorithm) MAK.check_input($f!, t, F, alg) @@ -45,6 +45,43 @@ for f! in ( end end +for alg in (MAK.LeftOrthAlgorithm{:qr}, MAK.LeftOrthAlgorithm{:svd}, MAK.LeftOrthAlgorithm{:polar}) + @eval begin + function MAK.left_orth!(t::AbstractTensorMap, F, alg::$alg) + MAK.check_input(left_orth!, t, F, alg) + + foreachblock(t, F...) do _, bs + factors = Base.tail(bs) + factors′ = left_orth!(first(bs), factors, alg) + # deal with the case where the output is not in-place + for (f′, f) in zip(factors′, factors) + f′ === f || copy!(f, f′) + end + return nothing + end + return F + end + end +end + +for alg in (MAK.RightOrthAlgorithm{:lq}, MAK.RightOrthAlgorithm{:svd}, MAK.RightOrthAlgorithm{:polar}) + @eval function MAK.right_orth!(t::AbstractTensorMap, F, alg::$alg) + MAK.check_input(right_orth!, t, F, alg) + + foreachblock(t, F...) do _, bs + factors = Base.tail(bs) + factors′ = right_orth!(first(bs), factors, alg) + # deal with the case where the output is not in-place + for (f′, f) in zip(factors′, factors) + f′ === f || copy!(f, f′) + end + return nothing + end + return F + end +end + + # Handle these separately because single output instead of tuple for f! in (:qr_null!, :lq_null!) @eval function MAK.$f!(t::AbstractTensorMap, N, alg::AbstractAlgorithm) @@ -464,6 +501,18 @@ function MAK.check_input(::typeof(left_orth!), t::AbstractTensorMap, VC, ::Abstr return nothing end +function MAK.check_input(::typeof(left_orth!), t::AbstractTensorMap, VC, alg::MAK.LeftOrthAlgorithm{:qr}) + return MAK.check_input(qr_compact!, t, VC, alg.alg) +end + +function MAK.check_input(::typeof(left_orth!), t::AbstractTensorMap, VC, alg::MAK.LeftOrthAlgorithm{:svd}) + return MAK.check_input(svd_compact!, t, VC, alg.alg) +end + +function MAK.check_input(::typeof(left_orth!), t::AbstractTensorMap, VC, alg::MAK.LeftOrthAlgorithm{:polar}) + return MAK.check_input(left_polar!, t, VC, alg.alg) +end + function MAK.check_input(::typeof(right_orth!), t::AbstractTensorMap, CVᴴ, ::AbstractAlgorithm) C, Vᴴ = CVᴴ @@ -479,6 +528,18 @@ function MAK.check_input(::typeof(right_orth!), t::AbstractTensorMap, CVᴴ, ::A return nothing end +function MAK.check_input(::typeof(right_orth!), t::AbstractTensorMap, VC, alg::MAK.RightOrthAlgorithm{:lq}) + return MAK.check_input(lq_compact!, t, VC, alg.alg) +end + +function MAK.check_input(::typeof(right_orth!), t::AbstractTensorMap, VC, alg::MAK.RightOrthAlgorithm{:svd}) + return MAK.check_input(svd_compact!, t, VC, alg.alg) +end + +function MAK.check_input(::typeof(right_orth!), t::AbstractTensorMap, VC, alg::MAK.RightOrthAlgorithm{:polar}) + return MAK.check_input(right_polar!, t, VC, alg.alg) +end + function MAK.initialize_output(::typeof(left_orth!), t::AbstractTensorMap) V_C = infimum(fuse(codomain(t)), fuse(domain(t))) V = similar(t, codomain(t) ← V_C) @@ -501,44 +562,16 @@ end function MAK.left_orth!( t::AbstractTensorMap; trunc::TruncationStrategy = notrunc(), - kind = trunc == notrunc() ? :qr : :svd, - alg_qr = (; positive = true), alg_polar = (;), alg_svd = (;) + alg::AbstractAlgorithm = (trunc == notrunc()) ? MAK.select_algorithm(left_orth!, t, Val(:qr)) : MAK.select_algorithm(left_orth!, t, Val(:svd); trunc) ) - trunc == notrunc() || kind === :svd || - throw(ArgumentError("truncation not supported for left_orth with kind = $kind")) - - return if kind === :qr - alg_qr isa NamedTuple ? qr_compact!(t; alg_qr...) : qr_compact!(t; alg = alg_qr) - elseif kind === :polar - alg_polar isa NamedTuple ? left_orth_polar!(t; alg_polar...) : - left_orth_polar!(t; alg = alg_polar) - elseif kind === :svd - alg_svd isa NamedTuple ? left_orth_svd!(t; trunc, alg_svd...) : - left_orth_svd!(t; trunc, alg = alg_svd) - else - throw(ArgumentError(lazy"`left_orth!` received unknown value `kind = $kind`")) - end + MAK.left_orth!(t, MAK.initialize_output(left_orth!, t, alg), alg) end function MAK.right_orth!( t::AbstractTensorMap; trunc::TruncationStrategy = notrunc(), - kind = trunc == notrunc() ? :lq : :svd, - alg_lq = (; positive = true), alg_polar = (;), alg_svd = (;) + alg::AbstractAlgorithm = (trunc == notrunc()) ? MAK.select_algorithm(right_orth!, t, Val(:lq)) : MAK.select_algorithm(right_orth!, t, Val(:svd); trunc) ) - trunc == notrunc() || kind === :svd || - throw(ArgumentError("truncation not supported for right_orth with kind = $kind")) - - return if kind === :lq - alg_lq isa NamedTuple ? lq_compact!(t; alg_lq...) : lq_compact!(t; alg = alg_lq) - elseif kind === :polar - alg_polar isa NamedTuple ? right_orth_polar!(t; alg_polar...) : - right_orth_polar!(t; alg = alg_polar) - elseif kind === :svd - alg_svd isa NamedTuple ? right_orth_svd!(t; trunc, alg_svd...) : - right_orth_svd!(t; trunc, alg = alg_svd) - else - throw(ArgumentError(lazy"`right_orth!` received unknown value `kind = $kind`")) - end + MAK.right_orth!(t, MAK.initialize_output(right_orth!, t, alg), alg) end # Nullspace diff --git a/test/autodiff/ad.jl b/test/autodiff/ad.jl index 15c8c387e..58b90c61c 100644 --- a/test/autodiff/ad.jl +++ b/test/autodiff/ad.jl @@ -103,7 +103,7 @@ function remove_lqgauge_dependence!(ΔQ, t, Q) return ΔQ end function remove_eiggauge_dependence!( - ΔV, D, V; degeneracy_atol = MatrixAlgebraKit.default_pullback_gaugetol(D) + ΔV, D, V; degeneracy_atol = MatrixAlgebraKit.default_pullback_gauge_atol(D) ) gaugepart = V' * ΔV for (c, b) in blocks(gaugepart) @@ -119,7 +119,7 @@ function remove_eiggauge_dependence!( return ΔV end function remove_eighgauge_dependence!( - ΔV, D, V; degeneracy_atol = MatrixAlgebraKit.default_pullback_gaugetol(D) + ΔV, D, V; degeneracy_atol = MatrixAlgebraKit.default_pullback_gauge_atol(D) ) gaugepart = V' * ΔV gaugepart = (gaugepart - gaugepart') / 2 @@ -136,7 +136,7 @@ function remove_eighgauge_dependence!( return ΔV end function remove_svdgauge_dependence!( - ΔU, ΔVᴴ, U, S, Vᴴ; degeneracy_atol = MatrixAlgebraKit.default_pullback_gaugetol(S) + ΔU, ΔVᴴ, U, S, Vᴴ; degeneracy_atol = MatrixAlgebraKit.default_pullback_gauge_atol(S) ) gaugepart = U' * ΔU + Vᴴ * ΔVᴴ' gaugepart = (gaugepart - gaugepart') / 2 diff --git a/test/tensors/factorizations.jl b/test/tensors/factorizations.jl index 8c9bde170..c7950267b 100644 --- a/test/tensors/factorizations.jl +++ b/test/tensors/factorizations.jl @@ -49,7 +49,7 @@ for V in spacelist @test Q * R ≈ t @test isisometric(Q) - Q, R = @constinferred left_orth(t; kind = :qr) + Q, R = @constinferred left_orth(t) @test Q * R ≈ t @test isisometric(Q) @@ -57,7 +57,7 @@ for V in spacelist @test isisometric(N) @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t)) - N = @constinferred left_null(t; kind = :qr) + N = @constinferred left_null(t) @test isisometric(N) @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t)) end @@ -76,7 +76,7 @@ for V in spacelist @test isisometric(Q) @test dim(Q) == dim(R) == dim(t) - Q, R = @constinferred left_orth(t; kind = :qr) + Q, R = @constinferred left_orth(t) @test Q * R ≈ t @test isisometric(Q) @test dim(Q) == dim(R) == dim(t) @@ -102,7 +102,7 @@ for V in spacelist @test L * Q ≈ t @test isisometric(Q; side = :right) - L, Q = @constinferred right_orth(t; kind = :lq) + L, Q = @constinferred right_orth(t) @test L * Q ≈ t @test isisometric(Q; side = :right) @@ -125,7 +125,7 @@ for V in spacelist @test isisometric(Q; side = :right) @test dim(Q) == dim(L) == dim(t) - L, Q = @constinferred right_orth(t; kind = :lq) + L, Q = @constinferred right_orth(t) @test L * Q ≈ t @test isisometric(Q; side = :right) @test dim(Q) == dim(L) == dim(t) @@ -149,7 +149,7 @@ for V in spacelist @test isisometric(w) @test isposdef(p) - w, p = @constinferred left_orth(t; kind = :polar) + w, p = @constinferred left_orth(t; alg = TensorKit.LeftOrthAlgorithm{:polar}) @test w * p ≈ t @test isisometric(w) end @@ -163,7 +163,7 @@ for V in spacelist @test isisometric(wᴴ; side = :right) @test isposdef(p) - p, wᴴ = @constinferred right_orth(t; kind = :polar) + p, wᴴ = @constinferred right_orth(t; alg = TensorKit.RightOrthAlgorithm{:polar}) @test p * wᴴ ≈ t @test isisometric(wᴴ; side = :right) end @@ -194,9 +194,13 @@ for V in spacelist @test b ≈ s′[c] end - v, c = @constinferred left_orth(t; kind = :svd) + v, c = @constinferred left_orth(t; alg = TensorKit.LeftOrthAlgorithm{:svd}) @test v * c ≈ t @test isisometric(v) + + c, vᴴ = @constinferred right_orth(t; alg = TensorKit.RightOrthAlgorithm{:svd}) + @test c * vᴴ ≈ t + @test isisometric(v; side = :right) N = @constinferred left_null(t; kind = :svd) @test isisometric(N) From f4412fccc8be032709c60eb0d6c93ca62a3ee43e Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Fri, 14 Nov 2025 16:21:06 +0100 Subject: [PATCH 07/41] Adjoint factorizations (AD not working) --- src/factorizations/adjoint.jl | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/factorizations/adjoint.jl b/src/factorizations/adjoint.jl index ef7e5e265..ac15dadb1 100644 --- a/src/factorizations/adjoint.jl +++ b/src/factorizations/adjoint.jl @@ -69,6 +69,24 @@ for (left_f!, right_f!) in zip( end end +for (left_f, right_f) in zip( + (:qr_full, :qr_compact, :left_polar, :left_orth), + (:lq_full, :lq_compact, :right_polar, :right_orth) + ) + @eval function MAK.$left_f(t::AdjointTensorMap; kwargs...) + return reverse(adjoint.($right_f(adjoint(t); kwargs...))) + end + @eval function MAK.$right_f(t::AdjointTensorMap; kwargs...) + return reverse(adjoint.($left_f(adjoint(t); kwargs...))) + end +end + +function MAK.eig_full(t::AdjointTensorMap; kwargs...) + DV = eig_full(adjoint(t); kwargs...) + return (DV[1], adjoint(DV[2])) +end + + # 3-arg functions for f! in (:svd_full!, :svd_compact!, :svd_trunc!) @eval function MAK.copy_input(::typeof($f!), t::AdjointTensorMap) From 12e549f1465c46335a1da3df222309473262a4c3 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 09:10:14 -0500 Subject: [PATCH 08/41] add missing MAK. --- src/factorizations/adjoint.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/factorizations/adjoint.jl b/src/factorizations/adjoint.jl index ac15dadb1..c3e96ec1a 100644 --- a/src/factorizations/adjoint.jl +++ b/src/factorizations/adjoint.jl @@ -30,10 +30,10 @@ function MAK.right_null!(t::AdjointTensorMap, N, alg::AbstractAlgorithm) end function MAK.is_left_isometric(t::AdjointTensorMap; kwargs...) - return is_right_isometric(adjoint(t); kwargs...) + return MAK.is_right_isometric(adjoint(t); kwargs...) end function MAK.is_right_isometric(t::AdjointTensorMap; kwargs...) - return is_left_isometric(adjoint(t); kwargs...) + return MAK.is_left_isometric(adjoint(t); kwargs...) end # 2-arg functions From f84d97154810d22d50626043678b75984263c78f Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 09:10:56 -0500 Subject: [PATCH 09/41] remove `eig(::Adjoint)` This doesn't actually hold for non-hermitian matrices --- src/factorizations/adjoint.jl | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/factorizations/adjoint.jl b/src/factorizations/adjoint.jl index c3e96ec1a..8b209d43d 100644 --- a/src/factorizations/adjoint.jl +++ b/src/factorizations/adjoint.jl @@ -81,12 +81,6 @@ for (left_f, right_f) in zip( end end -function MAK.eig_full(t::AdjointTensorMap; kwargs...) - DV = eig_full(adjoint(t); kwargs...) - return (DV[1], adjoint(DV[2])) -end - - # 3-arg functions for f! in (:svd_full!, :svd_compact!, :svd_trunc!) @eval function MAK.copy_input(::typeof($f!), t::AdjointTensorMap) From 9b37806fbfd8e09bab83a39a6998541df2f95cf5 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 09:11:21 -0500 Subject: [PATCH 10/41] add truncation error implementation --- src/factorizations/truncation.jl | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 19f39759c..478d5ac12 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -242,3 +242,17 @@ function MAK.findtruncated_svd(values::SectorDict, strategy::TruncationIntersect ) for c in intersect(map(keys, inds)...) ) end + +# Truncation error +# ---------------- +MAK.truncation_error(values::SectorDict, ind) = + MAK.truncation_error!(SectorDict(c => copy(v) for (c, v) in values), ind) + +function MAK.truncation_error!(values::SectorDict, ind) + for (c, v) in values + ind_c = get(ind, c, nothing) + isnothing(ind_c) && continue + v[ind_c] .= zero(eltype(v)) + end + return TensorKit._norm(values, 2, zero(real(eltype(valtype(values))))) +end From effb03e5080a7a73bc9849a2d031f8819c293bb4 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 09:11:37 -0500 Subject: [PATCH 11/41] add nullspace truncation specializations --- src/factorizations/truncation.jl | 38 +++++++++++++++++++++++++++++--- 1 file changed, 35 insertions(+), 3 deletions(-) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 478d5ac12..23fc4c985 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -76,9 +76,7 @@ function MAK.truncate( end function MAK.truncate( - ::typeof(left_null!), - (U, S)::Tuple{AbstractTensorMap, AbstractTensorMap}, - strategy::MatrixAlgebraKit.TruncationStrategy + ::typeof(left_null!), (U, S)::NTuple{2, AbstractTensorMap}, strategy::TruncationStrategy ) extended_S = SectorDict( c => vcat(diagview(b), zeros(eltype(b), max(0, size(b, 2) - size(b, 1)))) @@ -90,6 +88,40 @@ function MAK.truncate( truncate_domain!(Ũ, U, ind) return Ũ, ind end +function MAK.truncate( + ::typeof(right_null!), (S, Vᴴ)::NTuple{2, AbstractTensorMap}, strategy::TruncationStrategy + ) + extended_S = SectorDict( + c => vcat(diagview(b), zeros(eltype(b), max(0, size(b, 1) - size(b, 2)))) + for (c, b) in blocks(S) + ) + ind = MAK.findtruncated(extended_S, strategy) + V_truncated = truncate_space(space(Vᴴ, 1), ind) + Ṽᴴ = similar(Vᴴ, V_truncated ← domain(Vᴴ)) + truncate_codomain!(Ṽᴴ, Vᴴ, ind) + return Ṽᴴ, ind +end + +# special case `NoTruncation` for null: should keep exact zeros due to rectangularity +# need to specialize to avoid ambiguity with special case in MatrixAlgebraKit +function MAK.truncate( + ::typeof(left_null!), (U, S)::NTuple{2, AbstractTensorMap}, strategy::NoTruncation + ) + ind = SectorDict(c => (size(b, 2) + 1):size(b, 1) for (c, b) in blocks(S)) + V_truncated = truncate_space(space(S, 1), ind) + Ũ = similar(U, codomain(U) ← V_truncated) + truncate_domain!(Ũ, U, ind) + return Ũ, ind +end +function MAK.truncate( + ::typeof(right_null!), (S, Vᴴ)::NTuple{2, AbstractTensorMap}, strategy::NoTruncation + ) + ind = SectorDict(c => (size(b, 1) + 1):size(b, 2) for (c, b) in blocks(S)) + V_truncated = truncate_space(space(Vᴴ, 1), ind) + Ṽᴴ = similar(Vᴴ, V_truncated ← domain(Vᴴ)) + truncate_codomain!(Ṽᴴ, Vᴴ, ind) + return Ṽᴴ, ind +end for f! in (:eig_trunc!, :eigh_trunc!) @eval function MAK.truncate( From 9b408e6b7a1b17ffa26febdb6bc78d0c76571627 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 09:22:45 -0500 Subject: [PATCH 12/41] update orth interface in tests --- test/tensors/factorizations.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/test/tensors/factorizations.jl b/test/tensors/factorizations.jl index c7950267b..af1e9ef6c 100644 --- a/test/tensors/factorizations.jl +++ b/test/tensors/factorizations.jl @@ -149,7 +149,7 @@ for V in spacelist @test isisometric(w) @test isposdef(p) - w, p = @constinferred left_orth(t; alg = TensorKit.LeftOrthAlgorithm{:polar}) + w, p = @constinferred left_orth(t; alg = :polar) @test w * p ≈ t @test isisometric(w) end @@ -163,7 +163,7 @@ for V in spacelist @test isisometric(wᴴ; side = :right) @test isposdef(p) - p, wᴴ = @constinferred right_orth(t; alg = TensorKit.RightOrthAlgorithm{:polar}) + p, wᴴ = @constinferred right_orth(t; alg = :polar) @test p * wᴴ ≈ t @test isisometric(wᴴ; side = :right) end @@ -194,19 +194,19 @@ for V in spacelist @test b ≈ s′[c] end - v, c = @constinferred left_orth(t; alg = TensorKit.LeftOrthAlgorithm{:svd}) + v, c = @constinferred left_orth(t; alg = :svd) @test v * c ≈ t @test isisometric(v) - - c, vᴴ = @constinferred right_orth(t; alg = TensorKit.RightOrthAlgorithm{:svd}) + + c, vᴴ = @constinferred right_orth(t; alg = :svd) @test c * vᴴ ≈ t @test isisometric(v; side = :right) - N = @constinferred left_null(t; kind = :svd) + N = @constinferred left_null(t; alg = :svd) @test isisometric(N) @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t)) - Nᴴ = @constinferred right_null(t; kind = :svd) + Nᴴ = @constinferred right_null(t; alg = :svd) @test isisometric(Nᴴ; side = :right) @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t)) end From b8d9d13ff0b71874c9aac18ca22740de4a4a5fbf Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 09:22:54 -0500 Subject: [PATCH 13/41] fix wrong variable --- test/tensors/factorizations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/tensors/factorizations.jl b/test/tensors/factorizations.jl index af1e9ef6c..5ca880582 100644 --- a/test/tensors/factorizations.jl +++ b/test/tensors/factorizations.jl @@ -200,7 +200,7 @@ for V in spacelist c, vᴴ = @constinferred right_orth(t; alg = :svd) @test c * vᴴ ≈ t - @test isisometric(v; side = :right) + @test isisometric(vᴴ; side = :right) N = @constinferred left_null(t; alg = :svd) @test isisometric(N) From 6eb246921d24f5303d67870475e54de07df68660 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 09:23:58 -0500 Subject: [PATCH 14/41] improve adjoint support --- src/factorizations/adjoint.jl | 23 +++++++++++++++++++---- src/tensors/adjoint.jl | 2 ++ 2 files changed, 21 insertions(+), 4 deletions(-) diff --git a/src/factorizations/adjoint.jl b/src/factorizations/adjoint.jl index 8b209d43d..406fc8082 100644 --- a/src/factorizations/adjoint.jl +++ b/src/factorizations/adjoint.jl @@ -9,6 +9,21 @@ _adjoint(alg::MAK.LAPACK_HouseholderRQ) = MAK.LAPACK_HouseholderQL(; alg.kwargs. _adjoint(alg::MAK.PolarViaSVD) = MAK.PolarViaSVD(_adjoint(alg.svdalg)) _adjoint(alg::AbstractAlgorithm) = alg +for f in + [ + :svd_compact, :svd_full, :svd_trunc, :svd_vals, :qr_compact, :qr_full, :qr_null, + :lq_compact, :lq_full, :lq_null, :eig_full, :eig_trunc, :eig_vals, :eigh_full, + :eigh_trunc, :eigh_vals, :left_polar, :right_polar, + :project_hermitian, :project_antihermitian, :project_isometric, + ] + f! = Symbol(f, :!) + # just return the algorithm for the parent type since we are mapping this with + # `_adjoint` afterwards anyways. + # TODO: properly handle these cases + @eval MAK.default_algorithm(::typeof($f!), ::Type{T}; kwargs...) where {T <: AdjointTensorMap} = + MAK.default_algorithm($f!, TensorKit.parenttype(T); kwargs...) +end + # 1-arg functions function MAK.initialize_output(::typeof(left_null!), t::AdjointTensorMap, alg::AbstractAlgorithm) return adjoint(MAK.initialize_output(right_null!, adjoint(t), _adjoint(alg))) @@ -38,8 +53,8 @@ end # 2-arg functions for (left_f!, right_f!) in zip( - (:qr_full!, :qr_compact!, :left_polar!, :left_orth!), - (:lq_full!, :lq_compact!, :right_polar!, :right_orth!) + (:qr_full!, :qr_compact!, :left_polar!), + (:lq_full!, :lq_compact!, :right_polar!) ) @eval function MAK.copy_input(::typeof($left_f!), t::AdjointTensorMap) return adjoint(MAK.copy_input($right_f!, adjoint(t))) @@ -70,8 +85,8 @@ for (left_f!, right_f!) in zip( end for (left_f, right_f) in zip( - (:qr_full, :qr_compact, :left_polar, :left_orth), - (:lq_full, :lq_compact, :right_polar, :right_orth) + (:qr_full, :qr_compact, :left_polar), + (:lq_full, :lq_compact, :right_polar) ) @eval function MAK.$left_f(t::AdjointTensorMap; kwargs...) return reverse(adjoint.($right_f(adjoint(t); kwargs...))) diff --git a/src/tensors/adjoint.jl b/src/tensors/adjoint.jl index 2a229f2c5..ca484e77b 100644 --- a/src/tensors/adjoint.jl +++ b/src/tensors/adjoint.jl @@ -11,6 +11,8 @@ struct AdjointTensorMap{T, S, N₁, N₂, TT <: AbstractTensorMap{T, S, N₂, N parent::TT end Base.parent(t::AdjointTensorMap) = t.parent +parenttype(t::AdjointTensorMap) = parenttype(typeof(t)) +parenttype(::Type{AdjointTensorMap{T, S, N₁, N₂, TT}}) where {T, S, N₁, N₂, TT} = TT # Constructor: construct from taking adjoint of a tensor Base.adjoint(t::AdjointTensorMap) = parent(t) From 95ea696dd806c6f24efb5af25efc467d9870168d Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 09:26:16 -0500 Subject: [PATCH 15/41] remove unnecessary left/rightorth functions --- src/TensorKit.jl | 3 +- src/factorizations/factorizations.jl | 1 - src/factorizations/matrixalgebrakit.jl | 168 +------------------------ 3 files changed, 2 insertions(+), 170 deletions(-) diff --git a/src/TensorKit.jl b/src/TensorKit.jl index 7ebdb504d..85bc9a051 100644 --- a/src/TensorKit.jl +++ b/src/TensorKit.jl @@ -82,8 +82,7 @@ export left_orth, right_orth, left_null, right_null, eigh_full!, eigh_full, eigh_trunc!, eigh_trunc, eig_full!, eig_full, eig_trunc!, eig_trunc, eigh_vals!, eigh_vals, eig_vals!, eig_vals, - isposdef, isposdef!, ishermitian, isisometric, isunitary, sylvester, rank, cond, - LeftOrthAlgorithm, RightOrthAlgorithm + isposdef, isposdef!, ishermitian, isisometric, isunitary, sylvester, rank, cond export braid, braid!, permute, permute!, transpose, transpose!, twist, twist!, repartition, repartition! diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index d618bca09..6961d55aa 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -16,7 +16,6 @@ using TensorOperations: Index2Tuple using MatrixAlgebraKit import MatrixAlgebraKit as MAK using MatrixAlgebraKit: AbstractAlgorithm, TruncatedAlgorithm, DiagonalAlgorithm -using MatrixAlgebraKit: LeftOrthAlgorithm, RightOrthAlgorithm using MatrixAlgebraKit: TruncationStrategy, NoTruncation, TruncationByValue, TruncationByError, TruncationIntersection, TruncationByFilter, TruncationByOrder using MatrixAlgebraKit: diagview, isisometric diff --git a/src/factorizations/matrixalgebrakit.jl b/src/factorizations/matrixalgebrakit.jl index 73fd04a6b..3478a6d27 100644 --- a/src/factorizations/matrixalgebrakit.jl +++ b/src/factorizations/matrixalgebrakit.jl @@ -26,7 +26,7 @@ end for f! in ( :qr_compact!, :qr_full!, :lq_compact!, :lq_full!, :eig_full!, :eigh_full!, :svd_compact!, :svd_full!, - :left_polar!, :right_polar! + :left_polar!, :right_polar!, ) @eval function MAK.$f!(t::AbstractTensorMap, F, alg::AbstractAlgorithm) MAK.check_input($f!, t, F, alg) @@ -45,43 +45,6 @@ for f! in ( end end -for alg in (MAK.LeftOrthAlgorithm{:qr}, MAK.LeftOrthAlgorithm{:svd}, MAK.LeftOrthAlgorithm{:polar}) - @eval begin - function MAK.left_orth!(t::AbstractTensorMap, F, alg::$alg) - MAK.check_input(left_orth!, t, F, alg) - - foreachblock(t, F...) do _, bs - factors = Base.tail(bs) - factors′ = left_orth!(first(bs), factors, alg) - # deal with the case where the output is not in-place - for (f′, f) in zip(factors′, factors) - f′ === f || copy!(f, f′) - end - return nothing - end - return F - end - end -end - -for alg in (MAK.RightOrthAlgorithm{:lq}, MAK.RightOrthAlgorithm{:svd}, MAK.RightOrthAlgorithm{:polar}) - @eval function MAK.right_orth!(t::AbstractTensorMap, F, alg::$alg) - MAK.check_input(right_orth!, t, F, alg) - - foreachblock(t, F...) do _, bs - factors = Base.tail(bs) - factors′ = right_orth!(first(bs), factors, alg) - # deal with the case where the output is not in-place - for (f′, f) in zip(factors′, factors) - f′ === f || copy!(f, f′) - end - return nothing - end - return F - end -end - - # Handle these separately because single output instead of tuple for f! in (:qr_null!, :lq_null!) @eval function MAK.$f!(t::AbstractTensorMap, N, alg::AbstractAlgorithm) @@ -484,135 +447,6 @@ function MAK.initialize_output(::typeof(right_polar!), t::AbstractTensorMap, ::A return P, Wᴴ end -# Orthogonalization -# ----------------- -function MAK.check_input(::typeof(left_orth!), t::AbstractTensorMap, VC, ::AbstractAlgorithm) - V, C = VC - - # scalartype checks - @check_scalar V t - isnothing(C) || @check_scalar C t - - # space checks - V_C = infimum(fuse(codomain(t)), fuse(domain(t))) - @check_space(V, codomain(t) ← V_C) - isnothing(C) || @check_space(C, V_C ← domain(t)) - - return nothing -end - -function MAK.check_input(::typeof(left_orth!), t::AbstractTensorMap, VC, alg::MAK.LeftOrthAlgorithm{:qr}) - return MAK.check_input(qr_compact!, t, VC, alg.alg) -end - -function MAK.check_input(::typeof(left_orth!), t::AbstractTensorMap, VC, alg::MAK.LeftOrthAlgorithm{:svd}) - return MAK.check_input(svd_compact!, t, VC, alg.alg) -end - -function MAK.check_input(::typeof(left_orth!), t::AbstractTensorMap, VC, alg::MAK.LeftOrthAlgorithm{:polar}) - return MAK.check_input(left_polar!, t, VC, alg.alg) -end - -function MAK.check_input(::typeof(right_orth!), t::AbstractTensorMap, CVᴴ, ::AbstractAlgorithm) - C, Vᴴ = CVᴴ - - # scalartype checks - isnothing(C) || @check_scalar C t - @check_scalar Vᴴ t - - # space checks - V_C = infimum(fuse(codomain(t)), fuse(domain(t))) - isnothing(C) || @check_space(C, codomain(t) ← V_C) - @check_space(Vᴴ, V_C ← domain(t)) - - return nothing -end - -function MAK.check_input(::typeof(right_orth!), t::AbstractTensorMap, VC, alg::MAK.RightOrthAlgorithm{:lq}) - return MAK.check_input(lq_compact!, t, VC, alg.alg) -end - -function MAK.check_input(::typeof(right_orth!), t::AbstractTensorMap, VC, alg::MAK.RightOrthAlgorithm{:svd}) - return MAK.check_input(svd_compact!, t, VC, alg.alg) -end - -function MAK.check_input(::typeof(right_orth!), t::AbstractTensorMap, VC, alg::MAK.RightOrthAlgorithm{:polar}) - return MAK.check_input(right_polar!, t, VC, alg.alg) -end - -function MAK.initialize_output(::typeof(left_orth!), t::AbstractTensorMap) - V_C = infimum(fuse(codomain(t)), fuse(domain(t))) - V = similar(t, codomain(t) ← V_C) - C = similar(t, V_C ← domain(t)) - return V, C -end - -function MAK.initialize_output(::typeof(right_orth!), t::AbstractTensorMap) - V_C = infimum(fuse(codomain(t)), fuse(domain(t))) - C = similar(t, codomain(t) ← V_C) - Vᴴ = similar(t, V_C ← domain(t)) - return C, Vᴴ -end - -# This is a rework of the dispatch logic in order to avoid having to deal with having to -# allocate the output before knowing the kind of decomposition. In particular, here I disable -# providing output arguments for left_ and right_orth. -# This is mainly because polar decompositions have different shapes, and SVD for Diagonal -# also does -function MAK.left_orth!( - t::AbstractTensorMap; - trunc::TruncationStrategy = notrunc(), - alg::AbstractAlgorithm = (trunc == notrunc()) ? MAK.select_algorithm(left_orth!, t, Val(:qr)) : MAK.select_algorithm(left_orth!, t, Val(:svd); trunc) - ) - MAK.left_orth!(t, MAK.initialize_output(left_orth!, t, alg), alg) -end -function MAK.right_orth!( - t::AbstractTensorMap; - trunc::TruncationStrategy = notrunc(), - alg::AbstractAlgorithm = (trunc == notrunc()) ? MAK.select_algorithm(right_orth!, t, Val(:lq)) : MAK.select_algorithm(right_orth!, t, Val(:svd); trunc) - ) - MAK.right_orth!(t, MAK.initialize_output(right_orth!, t, alg), alg) -end - -# Nullspace -# --------- -function MAK.check_input(::typeof(left_null!), t::AbstractTensorMap, N, ::AbstractAlgorithm) - # scalartype checks - @check_scalar N t - - # space checks - V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) - V_N = ⊖(fuse(codomain(t)), V_Q) - @check_space(N, codomain(t) ← V_N) - - return nothing -end - -function MAK.check_input(::typeof(right_null!), t::AbstractTensorMap, N, ::AbstractAlgorithm) - @check_scalar N t - - # space checks - V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) - V_N = ⊖(fuse(domain(t)), V_Q) - @check_space(N, V_N ← domain(t)) - - return nothing -end - -function MAK.initialize_output(::typeof(left_null!), t::AbstractTensorMap) - V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) - V_N = ⊖(fuse(codomain(t)), V_Q) - N = similar(t, codomain(t) ← V_N) - return N -end - -function MAK.initialize_output(::typeof(right_null!), t::AbstractTensorMap) - V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) - V_N = ⊖(fuse(domain(t)), V_Q) - N = similar(t, V_N ← domain(t)) - return N -end - # Projections # ----------- function MAK.check_input(::typeof(project_hermitian!), tsrc::AbstractTensorMap, tdst::AbstractTensorMap) From fb207b61fff9ae00ab4fa6893ba34ba0203a188d Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 09:26:24 -0500 Subject: [PATCH 16/41] pass kwargs in isisometri --- src/factorizations/factorizations.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index 6961d55aa..f9cef84a1 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -29,9 +29,9 @@ include("pullbacks.jl") TensorKit.one!(A::AbstractMatrix) = MatrixAlgebraKit.one!(A) -function MatrixAlgebraKit.isisometric(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple) +function MatrixAlgebraKit.isisometric(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple; kwargs...) t = permute(t, (p₁, p₂); copy = false) - return isisometric(t) + return isisometric(t; kwargs...) end #------------------------------# From 23fa05924638c72ea4898418e9f6be8eef01c191 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 09:34:07 -0500 Subject: [PATCH 17/41] add tests for truncation error --- test/tensors/factorizations.jl | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/test/tensors/factorizations.jl b/test/tensors/factorizations.jl index 5ca880582..2207b6bae 100644 --- a/test/tensors/factorizations.jl +++ b/test/tensors/factorizations.jl @@ -235,48 +235,55 @@ for V in spacelist @constinferred normalize!(t) - U, S, Vᴴ = @constinferred svd_trunc(t; trunc = notrunc()) + U, S, Vᴴ, ϵ = @constinferred svd_trunc(t; trunc = notrunc()) @test U * S * Vᴴ ≈ t + @test ϵ ≈ 0 @test isisometric(U) @test isisometric(Vᴴ; side = :right) trunc = truncrank(dim(domain(S)) ÷ 2) - U1, S1, Vᴴ1 = @constinferred svd_trunc(t; trunc) + U1, S1, Vᴴ1, ϵ1 = @constinferred svd_trunc(t; trunc) @test t * Vᴴ1' ≈ U1 * S1 @test isisometric(U1) @test isisometric(Vᴴ1; side = :right) + @test norm(t - U1 * S1 * Vᴴ1) ≈ ϵ1 atol = eps(real(T))^(4 / 5) @test dim(domain(S1)) <= trunc.howmany λ = minimum(minimum, values(LinearAlgebra.diag(S1))) trunc = trunctol(; atol = λ - 10eps(λ)) - U2, S2, Vᴴ2 = @constinferred svd_trunc(t; trunc) + U2, S2, Vᴴ2, ϵ2 = @constinferred svd_trunc(t; trunc) @test t * Vᴴ2' ≈ U2 * S2 @test isisometric(U2) @test isisometric(Vᴴ2; side = :right) + @test norm(t - U2 * S2 * Vᴴ2) ≈ ϵ2 atol = eps(real(T))^(4 / 5) @test minimum(minimum, values(LinearAlgebra.diag(S1))) >= λ @test U2 ≈ U1 @test S2 ≈ S1 @test Vᴴ2 ≈ Vᴴ1 + @test ϵ1 ≈ ϵ2 trunc = truncspace(space(S2, 1)) - U3, S3, Vᴴ3 = @constinferred svd_trunc(t; trunc) + U3, S3, Vᴴ3, ϵ3 = @constinferred svd_trunc(t; trunc) @test t * Vᴴ3' ≈ U3 * S3 @test isisometric(U3) @test isisometric(Vᴴ3; side = :right) + @test norm(t - U3 * S3 * Vᴴ3) ≈ ϵ3 atol = eps(real(T))^(4 / 5) @test space(S3, 1) ≾ space(S2, 1) - trunc = truncerror(; atol = 0.5) - U4, S4, Vᴴ4 = @constinferred svd_trunc(t; trunc) + trunc = truncerror(; atol = ϵ2) + U4, S4, Vᴴ4, ϵ4 = @constinferred svd_trunc(t; trunc) @test t * Vᴴ4' ≈ U4 * S4 @test isisometric(U4) @test isisometric(Vᴴ4; side = :right) - @test norm(t - U4 * S4 * Vᴴ4) <= 0.5 + @test norm(t - U4 * S4 * Vᴴ4) ≈ ϵ4 atol = eps(real(T))^(4 / 5) + @test ϵ4 ≤ ϵ2 trunc = truncrank(dim(domain(S)) ÷ 2) & trunctol(; atol = λ - 10eps(λ)) - U5, S5, Vᴴ5 = @constinferred svd_trunc(t; trunc) + U5, S5, Vᴴ5, ϵ5 = @constinferred svd_trunc(t; trunc) @test t * Vᴴ5' ≈ U5 * S5 @test isisometric(U5) @test isisometric(Vᴴ5; side = :right) + @test norm(t - U5 * S5 * Vᴴ5) ≈ ϵ5 atol = eps(real(T))^(4 / 5) @test minimum(minimum, values(LinearAlgebra.diag(S5))) >= λ @test dim(domain(S5)) ≤ dim(domain(S)) ÷ 2 end From fb547a8971cc8eaaa05e712e85afd139e67d9eab Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 09:48:00 -0500 Subject: [PATCH 18/41] more nullspace tests and fixes --- src/factorizations/truncation.jl | 8 ++++---- test/tensors/factorizations.jl | 8 ++++++++ 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 23fc4c985..1940ec6de 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -79,7 +79,7 @@ function MAK.truncate( ::typeof(left_null!), (U, S)::NTuple{2, AbstractTensorMap}, strategy::TruncationStrategy ) extended_S = SectorDict( - c => vcat(diagview(b), zeros(eltype(b), max(0, size(b, 2) - size(b, 1)))) + c => vcat(diagview(b), zeros(eltype(b), max(0, size(b, 1) - size(b, 2)))) for (c, b) in blocks(S) ) ind = MAK.findtruncated(extended_S, strategy) @@ -92,11 +92,11 @@ function MAK.truncate( ::typeof(right_null!), (S, Vᴴ)::NTuple{2, AbstractTensorMap}, strategy::TruncationStrategy ) extended_S = SectorDict( - c => vcat(diagview(b), zeros(eltype(b), max(0, size(b, 1) - size(b, 2)))) + c => vcat(diagview(b), zeros(eltype(b), max(0, size(b, 2) - size(b, 1)))) for (c, b) in blocks(S) ) ind = MAK.findtruncated(extended_S, strategy) - V_truncated = truncate_space(space(Vᴴ, 1), ind) + V_truncated = truncate_space(dual(space(S, 2)), ind) Ṽᴴ = similar(Vᴴ, V_truncated ← domain(Vᴴ)) truncate_codomain!(Ṽᴴ, Vᴴ, ind) return Ṽᴴ, ind @@ -117,7 +117,7 @@ function MAK.truncate( ::typeof(right_null!), (S, Vᴴ)::NTuple{2, AbstractTensorMap}, strategy::NoTruncation ) ind = SectorDict(c => (size(b, 1) + 1):size(b, 2) for (c, b) in blocks(S)) - V_truncated = truncate_space(space(Vᴴ, 1), ind) + V_truncated = truncate_space(dual(space(S, 2)), ind) Ṽᴴ = similar(Vᴴ, V_truncated ← domain(Vᴴ)) truncate_codomain!(Ṽᴴ, Vᴴ, ind) return Ṽᴴ, ind diff --git a/test/tensors/factorizations.jl b/test/tensors/factorizations.jl index 2207b6bae..551c382d1 100644 --- a/test/tensors/factorizations.jl +++ b/test/tensors/factorizations.jl @@ -206,9 +206,17 @@ for V in spacelist @test isisometric(N) @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t)) + N = @constinferred left_null(t; trunc = (; atol = 100 * eps(norm(t)))) + @test isisometric(N) + @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t)) + Nᴴ = @constinferred right_null(t; alg = :svd) @test isisometric(Nᴴ; side = :right) @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t)) + + Nᴴ = @constinferred right_null(t; trunc = (; atol = 100 * eps(norm(t)))) + @test isisometric(Nᴴ; side = :right) + @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t)) end # empty tensor From a667e00961f19405c3da43311f42d19fff3f1fbf Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 09:56:52 -0500 Subject: [PATCH 19/41] remove unnecessary specialization --- src/factorizations/matrixalgebrakit.jl | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/factorizations/matrixalgebrakit.jl b/src/factorizations/matrixalgebrakit.jl index 3478a6d27..8bb9fe536 100644 --- a/src/factorizations/matrixalgebrakit.jl +++ b/src/factorizations/matrixalgebrakit.jl @@ -149,11 +149,6 @@ function MAK.initialize_output(::typeof(svd_compact!), t::AbstractTensorMap, ::A return U, S, Vᴴ end -# TODO: remove this once `AbstractMatrix` specialization is removed in MatrixAlgebraKit -function MAK.initialize_output(::typeof(svd_trunc!), t::AbstractTensorMap, alg::TruncatedAlgorithm) - return MAK.initialize_output(svd_compact!, t, alg.alg) -end - function MAK.initialize_output(::typeof(svd_vals!), t::AbstractTensorMap, alg::AbstractAlgorithm) V_cod = infimum(fuse(codomain(t)), fuse(domain(t))) return DiagonalTensorMap{real(scalartype(t))}(undef, V_cod) From 8da18070947bb915b13580b89e3f1edf05f17e11 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 09:59:44 -0500 Subject: [PATCH 20/41] add missing argument --- src/factorizations/truncation.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 1940ec6de..f08a15569 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -145,7 +145,8 @@ end # Find truncation # --------------- # auxiliary functions -rtol_to_atol(S, p, atol, rtol) = rtol > 0 ? max(atol, TensorKit._norm(S, p) * rtol) : atol +rtol_to_atol(S, p, atol, rtol) = + rtol == 0 ? atol : max(atol, TensorKit._norm(S, p, norm(zero(scalartype(valtype(S))))) * rtol) function _compute_truncerr(Σdata, truncdim, p = 2) I = keytype(Σdata) From d7538236603cb077ab77ba5f772d440ced6fff44 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 10:04:49 -0500 Subject: [PATCH 21/41] remove unnecessary specializations --- src/factorizations/adjoint.jl | 28 +++++++--------------------- 1 file changed, 7 insertions(+), 21 deletions(-) diff --git a/src/factorizations/adjoint.jl b/src/factorizations/adjoint.jl index 406fc8082..f2fc3548f 100644 --- a/src/factorizations/adjoint.jl +++ b/src/factorizations/adjoint.jl @@ -52,10 +52,12 @@ function MAK.is_right_isometric(t::AdjointTensorMap; kwargs...) end # 2-arg functions -for (left_f!, right_f!) in zip( - (:qr_full!, :qr_compact!, :left_polar!), - (:lq_full!, :lq_compact!, :right_polar!) +for (left_f, right_f) in zip( + (:qr_full, :qr_compact, :left_polar), + (:lq_full, :lq_compact, :right_polar) ) + left_f! = Symbol(left_f, :!) + right_f! = Symbol(right_f, :!) @eval function MAK.copy_input(::typeof($left_f!), t::AdjointTensorMap) return adjoint(MAK.copy_input($right_f!, adjoint(t))) end @@ -74,26 +76,10 @@ for (left_f!, right_f!) in zip( return reverse(adjoint.(MAK.initialize_output($left_f!, adjoint(t), _adjoint(alg)))) end - @eval function MAK.$left_f!(t::AdjointTensorMap, F, alg::AbstractAlgorithm) + @eval MAK.$left_f!(t::AdjointTensorMap, F, alg::AbstractAlgorithm) = $right_f!(adjoint(t), reverse(adjoint.(F)), _adjoint(alg)) - return F - end - @eval function MAK.$right_f!(t::AdjointTensorMap, F, alg::AbstractAlgorithm) + @eval MAK.$right_f!(t::AdjointTensorMap, F, alg::AbstractAlgorithm) = $left_f!(adjoint(t), reverse(adjoint.(F)), _adjoint(alg)) - return F - end -end - -for (left_f, right_f) in zip( - (:qr_full, :qr_compact, :left_polar), - (:lq_full, :lq_compact, :right_polar) - ) - @eval function MAK.$left_f(t::AdjointTensorMap; kwargs...) - return reverse(adjoint.($right_f(adjoint(t); kwargs...))) - end - @eval function MAK.$right_f(t::AdjointTensorMap; kwargs...) - return reverse(adjoint.($left_f(adjoint(t); kwargs...))) - end end # 3-arg functions From 812f386011f8dee40932750f61263d5bbeeea037 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 10:09:11 -0500 Subject: [PATCH 22/41] slight reorganization --- src/factorizations/factorizations.jl | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index f9cef84a1..33f267e27 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -18,7 +18,7 @@ import MatrixAlgebraKit as MAK using MatrixAlgebraKit: AbstractAlgorithm, TruncatedAlgorithm, DiagonalAlgorithm using MatrixAlgebraKit: TruncationStrategy, NoTruncation, TruncationByValue, TruncationByError, TruncationIntersection, TruncationByFilter, TruncationByOrder -using MatrixAlgebraKit: diagview, isisometric +using MatrixAlgebraKit: diagview include("utility.jl") include("matrixalgebrakit.jl") @@ -29,11 +29,6 @@ include("pullbacks.jl") TensorKit.one!(A::AbstractMatrix) = MatrixAlgebraKit.one!(A) -function MatrixAlgebraKit.isisometric(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple; kwargs...) - t = permute(t, (p₁, p₂); copy = false) - return isisometric(t; kwargs...) -end - #------------------------------# # LinearAlgebra overloads #------------------------------# @@ -96,5 +91,9 @@ function MAK.is_right_isometric(t::AbstractTensorMap; kwargs...) f((c, b)) = MAK.is_right_isometric(b; kwargs...) return all(f, blocks(t)) end +function MAK.isisometric(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple; kwargs...) + t = permute(t, (p₁, p₂); copy = false) + return MAK.isisometric(t; kwargs...) +end end From aea45652031edeb4017d87041e9d4f9b5dfe7b19 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 10:09:16 -0500 Subject: [PATCH 23/41] export new functionality --- src/TensorKit.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/TensorKit.jl b/src/TensorKit.jl index 85bc9a051..ded08e50e 100644 --- a/src/TensorKit.jl +++ b/src/TensorKit.jl @@ -79,10 +79,12 @@ export left_orth, right_orth, left_null, right_null, qr_full!, qr_compact!, qr_null!, lq_full!, lq_compact!, lq_null!, svd_compact!, svd_full!, svd_trunc!, svd_compact, svd_full, svd_trunc, exp, exp!, - eigh_full!, eigh_full, eigh_trunc!, eigh_trunc, eig_full!, eig_full, eig_trunc!, - eig_trunc, - eigh_vals!, eigh_vals, eig_vals!, eig_vals, - isposdef, isposdef!, ishermitian, isisometric, isunitary, sylvester, rank, cond + eigh_full!, eigh_full, eigh_trunc!, eigh_trunc, eigh_vals!, eigh_vals, + eig_full!, eig_full, eig_trunc!, eig_trunc, eig_vals!, eig_vals, + ishermitian, project_hermitian, project_hermitian!, + isantihermitian, project_antihermitian, project_antihermitian!, + isisometric, isunitary, project_isometric, project_isometric!, + isposdef, isposdef!, sylvester, rank, cond export braid, braid!, permute, permute!, transpose, transpose!, twist, twist!, repartition, repartition! From 3a1a2597503f15f61ffa4f1b42712439ec84c618 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 20:25:04 -0500 Subject: [PATCH 24/41] simplify truncationerror implementation --- src/factorizations/truncation.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index f08a15569..0e9e2dda8 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -282,9 +282,8 @@ MAK.truncation_error(values::SectorDict, ind) = MAK.truncation_error!(SectorDict(c => copy(v) for (c, v) in values), ind) function MAK.truncation_error!(values::SectorDict, ind) - for (c, v) in values - ind_c = get(ind, c, nothing) - isnothing(ind_c) && continue + for (c, ind_c) in ind + v = values[c] v[ind_c] .= zero(eltype(v)) end return TensorKit._norm(values, 2, zero(real(eltype(valtype(values))))) From e7a472fd67c6667781e02020de01b29c7ecf158d Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 20:25:11 -0500 Subject: [PATCH 25/41] fix testcase --- test/tensors/tensors.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/tensors/tensors.jl b/test/tensors/tensors.jl index 51c326cf9..34934eb34 100644 --- a/test/tensors/tensors.jl +++ b/test/tensors/tensors.jl @@ -364,7 +364,7 @@ for V in spacelist for T in (Float64, ComplexF64) t1 = randisometry(T, W1, W2) t2 = randisometry(T, W2 ← W2) - @test isisometry(t1) + @test isisometric(t1) @test isunitary(t2) P = t1 * t1' @test P * P ≈ P From 8a6365ac538f2972efbd86c90c229cdb288e4245 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 20:32:23 -0500 Subject: [PATCH 26/41] refrain from defining left_null functions to avoid ambiguities --- src/factorizations/adjoint.jl | 35 ++++++++++++----------------------- 1 file changed, 12 insertions(+), 23 deletions(-) diff --git a/src/factorizations/adjoint.jl b/src/factorizations/adjoint.jl index f2fc3548f..41a0b9513 100644 --- a/src/factorizations/adjoint.jl +++ b/src/factorizations/adjoint.jl @@ -25,31 +25,20 @@ for f in end # 1-arg functions -function MAK.initialize_output(::typeof(left_null!), t::AdjointTensorMap, alg::AbstractAlgorithm) - return adjoint(MAK.initialize_output(right_null!, adjoint(t), _adjoint(alg))) -end -function MAK.initialize_output( - ::typeof(right_null!), t::AdjointTensorMap, - alg::AbstractAlgorithm - ) - return adjoint(MAK.initialize_output(left_null!, adjoint(t), _adjoint(alg))) -end +MAK.initialize_output(::typeof(qr_null!), t::AdjointTensorMap, alg::AbstractAlgorithm) = + adjoint(MAK.initialize_output(lq_null!, adjoint(t), _adjoint(alg))) +MAK.initialize_output(::typeof(lq_null!), t::AdjointTensorMap, alg::AbstractAlgorithm) = + adjoint(MAK.initialize_output(qr_null!, adjoint(t), _adjoint(alg))) -function MAK.left_null!(t::AdjointTensorMap, N, alg::AbstractAlgorithm) - right_null!(adjoint(t), adjoint(N), _adjoint(alg)) - return N -end -function MAK.right_null!(t::AdjointTensorMap, N, alg::AbstractAlgorithm) - left_null!(adjoint(t), adjoint(N), _adjoint(alg)) - return N -end +MAK.qr_null!(t::AdjointTensorMap, N, alg::AbstractAlgorithm) = + lq_null!(adjoint(t), adjoint(N), _adjoint(alg)) +MAK.lq_null!(t::AdjointTensorMap, N, alg::AbstractAlgorithm) = + qr_null!(adjoint(t), adjoint(N), _adjoint(alg)) -function MAK.is_left_isometric(t::AdjointTensorMap; kwargs...) - return MAK.is_right_isometric(adjoint(t); kwargs...) -end -function MAK.is_right_isometric(t::AdjointTensorMap; kwargs...) - return MAK.is_left_isometric(adjoint(t); kwargs...) -end +MAK.is_left_isometric(t::AdjointTensorMap; kwargs...) = + MAK.is_right_isometric(adjoint(t); kwargs...) +MAK.is_right_isometric(t::AdjointTensorMap; kwargs...) = + MAK.is_left_isometric(adjoint(t); kwargs...) # 2-arg functions for (left_f, right_f) in zip( From ea30e60b1a8c0c0622d6c7f1f2e9ff4a8f89c0c1 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 20:48:07 -0500 Subject: [PATCH 27/41] fix ad tests --- test/autodiff/ad.jl | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/test/autodiff/ad.jl b/test/autodiff/ad.jl index 58b90c61c..e72b74873 100644 --- a/test/autodiff/ad.jl +++ b/test/autodiff/ad.jl @@ -90,7 +90,6 @@ function remove_qrgauge_dependence!(ΔQ, t, Q) end return ΔQ end - function remove_lqgauge_dependence!(ΔQ, t, Q) for (c, b) in blocks(ΔQ) m, n = size(block(t, c)) @@ -103,7 +102,7 @@ function remove_lqgauge_dependence!(ΔQ, t, Q) return ΔQ end function remove_eiggauge_dependence!( - ΔV, D, V; degeneracy_atol = MatrixAlgebraKit.default_pullback_gauge_atol(D) + ΔV, D, V; degeneracy_atol = MatrixAlgebraKit.default_pullback_degeneracy_atol(D) ) gaugepart = V' * ΔV for (c, b) in blocks(gaugepart) @@ -119,9 +118,9 @@ function remove_eiggauge_dependence!( return ΔV end function remove_eighgauge_dependence!( - ΔV, D, V; degeneracy_atol = MatrixAlgebraKit.default_pullback_gauge_atol(D) + ΔV, D, V; degeneracy_atol = MatrixAlgebraKit.default_pullback_degeneracy_atol(D) ) - gaugepart = V' * ΔV + gaugepart = project_antihermitian!(V' * ΔV) gaugepart = (gaugepart - gaugepart') / 2 for (c, b) in blocks(gaugepart) Dc = diagview(block(D, c)) @@ -136,10 +135,9 @@ function remove_eighgauge_dependence!( return ΔV end function remove_svdgauge_dependence!( - ΔU, ΔVᴴ, U, S, Vᴴ; degeneracy_atol = MatrixAlgebraKit.default_pullback_gauge_atol(S) + ΔU, ΔVᴴ, U, S, Vᴴ; degeneracy_atol = MatrixAlgebraKit.default_pullback_degeneracy_atol(S) ) - gaugepart = U' * ΔU + Vᴴ * ΔVᴴ' - gaugepart = (gaugepart - gaugepart') / 2 + gaugepart = project_antihermitian!(U' * ΔU + Vᴴ * ΔVᴴ') for (c, b) in blocks(gaugepart) Sd = diagview(block(S, c)) # for some reason this fails only on tests, and I cannot reproduce it in an @@ -153,8 +151,6 @@ function remove_svdgauge_dependence!( return ΔU, ΔVᴴ end -project_hermitian(A) = (A + A') / 2 - # Tests # ----- @@ -602,7 +598,7 @@ for V in spacelist USVᴴ_trunc = svd_trunc(t; trunc) ΔUSVᴴ_trunc = rand_tangent.(USVᴴ_trunc) remove_svdgauge_dependence!( - ΔUSVᴴ_trunc[1], ΔUSVᴴ_trunc[3], USVᴴ_trunc...; degeneracy_atol + ΔUSVᴴ_trunc[1], ΔUSVᴴ_trunc[3], Base.front(USVᴴ_trunc)...; degeneracy_atol ) # test_ad_rrule(svd_trunc, t; # fkwargs=(; trunc), output_tangent=ΔUSVᴴ_trunc, atol, rtol) @@ -611,7 +607,7 @@ for V in spacelist USVᴴ_trunc = svd_trunc(t; trunc) ΔUSVᴴ_trunc = rand_tangent.(USVᴴ_trunc) remove_svdgauge_dependence!( - ΔUSVᴴ_trunc[1], ΔUSVᴴ_trunc[3], USVᴴ_trunc...; degeneracy_atol + ΔUSVᴴ_trunc[1], ΔUSVᴴ_trunc[3], Base.front(USVᴴ_trunc)...; degeneracy_atol ) test_ad_rrule( svd_trunc, t; @@ -632,7 +628,7 @@ for V in spacelist USVᴴ_trunc = svd_trunc(t; trunc) ΔUSVᴴ_trunc = rand_tangent.(USVᴴ_trunc) remove_svdgauge_dependence!( - ΔUSVᴴ_trunc[1], ΔUSVᴴ_trunc[3], USVᴴ_trunc...; degeneracy_atol + ΔUSVᴴ_trunc[1], ΔUSVᴴ_trunc[3], Base.front(USVᴴ_trunc)...; degeneracy_atol ) test_ad_rrule( svd_trunc, t; From 0c056b505954d5d23216b5c0c24c5c9caf98af71 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 21:49:59 -0500 Subject: [PATCH 28/41] add projection tests --- src/TensorKit.jl | 2 +- src/factorizations/matrixalgebrakit.jl | 22 +++++++--- test/tensors/factorizations.jl | 57 ++++++++++++++++++++++++++ 3 files changed, 74 insertions(+), 7 deletions(-) diff --git a/src/TensorKit.jl b/src/TensorKit.jl index ded08e50e..c6e003450 100644 --- a/src/TensorKit.jl +++ b/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 diff --git a/src/factorizations/matrixalgebrakit.jl b/src/factorizations/matrixalgebrakit.jl index 8bb9fe536..72d7015c0 100644 --- a/src/factorizations/matrixalgebrakit.jl +++ b/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) diff --git a/test/tensors/factorizations.jl b/test/tensors/factorizations.jl index 551c382d1..2c1638a30 100644 --- a/test/tensors/factorizations.jl +++ b/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 From edddfc2fa8a3fa8bf6a5e6ff52a11219f95ab0e8 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 15 Nov 2025 22:12:16 -0500 Subject: [PATCH 29/41] temporary AD fixes --- test/autodiff/ad.jl | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/test/autodiff/ad.jl b/test/autodiff/ad.jl index e72b74873..c6824503e 100644 --- a/test/autodiff/ad.jl +++ b/test/autodiff/ad.jl @@ -77,6 +77,9 @@ function test_ad_rrule(f, args...; check_inferred = false, kwargs...) return nothing end +# project_hermitian is non-differentiable for now +_project_hermitian(x) = (x + x') / 2 + # Gauge fixing tangents # --------------------- function remove_qrgauge_dependence!(ΔQ, t, Q) @@ -565,7 +568,7 @@ for V in spacelist remove_eighgauge_dependence!(Δv, d, v) # necessary for FiniteDifferences to not complain - eigh_full′ = eigh_full ∘ project_hermitian + eigh_full′ = eigh_full ∘ _project_hermitian test_ad_rrule(eigh_full′, t; output_tangent = (Δd, Δv), atol, rtol) test_ad_rrule(first ∘ eigh_full′, t; output_tangent = Δd, atol, rtol) @@ -596,7 +599,7 @@ for V in spacelist trunc = truncrank(max(2, round(Int, min(dim(domain(t)), dim(codomain(t))) * (3 / 4)))) USVᴴ_trunc = svd_trunc(t; trunc) - ΔUSVᴴ_trunc = rand_tangent.(USVᴴ_trunc) + ΔUSVᴴ_trunc = (rand_tangent.(Base.front(USVᴴ_trunc))..., zero(last(USVᴴ_trunc))) remove_svdgauge_dependence!( ΔUSVᴴ_trunc[1], ΔUSVᴴ_trunc[3], Base.front(USVᴴ_trunc)...; degeneracy_atol ) @@ -605,7 +608,7 @@ for V in spacelist trunc = truncspace(space(USVᴴ_trunc[2], 1)) USVᴴ_trunc = svd_trunc(t; trunc) - ΔUSVᴴ_trunc = rand_tangent.(USVᴴ_trunc) + ΔUSVᴴ_trunc = (rand_tangent.(Base.front(USVᴴ_trunc))..., zero(last(USVᴴ_trunc))) remove_svdgauge_dependence!( ΔUSVᴴ_trunc[1], ΔUSVᴴ_trunc[3], Base.front(USVᴴ_trunc)...; degeneracy_atol ) @@ -626,14 +629,14 @@ for V in spacelist tol = minimum(((c, b),) -> minimum(diagview(b)), blocks(USVᴴ_trunc[2])) trunc = trunctol(; atol = 10 * tol) USVᴴ_trunc = svd_trunc(t; trunc) - ΔUSVᴴ_trunc = rand_tangent.(USVᴴ_trunc) + ΔUSVᴴ_trunc = (rand_tangent.(Base.front(USVᴴ_trunc))..., zero(last(USVᴴ_trunc))) remove_svdgauge_dependence!( ΔUSVᴴ_trunc[1], ΔUSVᴴ_trunc[3], Base.front(USVᴴ_trunc)...; degeneracy_atol ) - test_ad_rrule( - svd_trunc, t; - fkwargs = (; trunc), output_tangent = ΔUSVᴴ_trunc, atol, rtol - ) + # test_ad_rrule( + # svd_trunc, t; + # fkwargs = (; trunc), output_tangent = ΔUSVᴴ_trunc, atol, rtol + # ) end end From dbb132d03e8bab638399a5457ad3fb4664710ff8 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 16 Nov 2025 08:55:29 -0500 Subject: [PATCH 30/41] try and add back some AD tests --- test/autodiff/ad.jl | 55 +++++++++++++++++++++------------------------ 1 file changed, 26 insertions(+), 29 deletions(-) diff --git a/test/autodiff/ad.jl b/test/autodiff/ad.jl index c6824503e..c3738ada7 100644 --- a/test/autodiff/ad.jl +++ b/test/autodiff/ad.jl @@ -594,19 +594,12 @@ for V in spacelist test_ad_rrule(svd_compact, t; output_tangent = (ΔU, ΔS, ΔVᴴ), atol, rtol) test_ad_rrule(svd_compact, t; output_tangent = (ΔU, ΔS2, ΔVᴴ), atol, rtol) - # TODO: I'm not sure how to properly test with spaces that might change - # with the finite-difference methods, as then the jacobian is ill-defined. - - trunc = truncrank(max(2, round(Int, min(dim(domain(t)), dim(codomain(t))) * (3 / 4)))) - USVᴴ_trunc = svd_trunc(t; trunc) - ΔUSVᴴ_trunc = (rand_tangent.(Base.front(USVᴴ_trunc))..., zero(last(USVᴴ_trunc))) - remove_svdgauge_dependence!( - ΔUSVᴴ_trunc[1], ΔUSVᴴ_trunc[3], Base.front(USVᴴ_trunc)...; degeneracy_atol - ) - # test_ad_rrule(svd_trunc, t; - # fkwargs=(; trunc), output_tangent=ΔUSVᴴ_trunc, atol, rtol) - - trunc = truncspace(space(USVᴴ_trunc[2], 1)) + # Testing truncation with finitedifferences is RNG-prone since the + # Jacobian changes size if the truncation space changes, causing errors. + # So, first test the fixed space case, then do more limited testing on + # some gradients and compare to the fixed space case + V_trunc = spacetype(t)(c => ceil(Int, min(size(b)...) / 2) for (c, b) in blocks(t)) + trunc = truncspace(V_trunc) USVᴴ_trunc = svd_trunc(t; trunc) ΔUSVᴴ_trunc = (rand_tangent.(Base.front(USVᴴ_trunc))..., zero(last(USVᴴ_trunc))) remove_svdgauge_dependence!( @@ -617,26 +610,30 @@ for V in spacelist fkwargs = (; trunc), output_tangent = ΔUSVᴴ_trunc, atol, rtol ) - # ϵ = norm(*(USVᴴ_trunc...) - t) - # trunc = truncerror(; atol=ϵ) - # USVᴴ_trunc = svd_trunc(t; trunc) - # ΔUSVᴴ_trunc = rand_tangent.(USVᴴ_trunc) - # remove_svdgauge_dependence!(ΔUSVᴴ_trunc[1], ΔUSVᴴ_trunc[3], USVᴴ_trunc...; - # degeneracy_atol) - # test_ad_rrule(svd_trunc, t; - # fkwargs=(; trunc), output_tangent=ΔUSVᴴ_trunc, atol, rtol) + # attempt to construct a loss function that doesn't depend on the gauges + function f(t; trunc) + Utr, Str, Vᴴtr, ϵ = svd_trunc(t; trunc) + return LinearAlgebra.tr(Str) + LinearAlgebra.norm(Utr * Vᴴtr) + end + + trunc = truncrank(round(Int, dim(V_trunc))) + USVᴴ_trunc = svd_trunc(t; trunc) + g1, = Zygote.gradient(x -> f(x; trunc), t) + g2, = Zygote.gradient(x -> f(x; trunc = truncspace(space(USVᴴ_trunc[2], 1))), t) + @test g1 ≈ g2 + + trunc = truncerror(; atol = last(USVᴴ_trunc)) + USVᴴ_trunc = svd_trunc(t; trunc) + g1, = Zygote.gradient(x -> f(x; trunc), t) + g2, = Zygote.gradient(x -> f(x; trunc = truncspace(space(USVᴴ_trunc[2], 1))), t) + @test g1 ≈ g2 tol = minimum(((c, b),) -> minimum(diagview(b)), blocks(USVᴴ_trunc[2])) trunc = trunctol(; atol = 10 * tol) USVᴴ_trunc = svd_trunc(t; trunc) - ΔUSVᴴ_trunc = (rand_tangent.(Base.front(USVᴴ_trunc))..., zero(last(USVᴴ_trunc))) - remove_svdgauge_dependence!( - ΔUSVᴴ_trunc[1], ΔUSVᴴ_trunc[3], Base.front(USVᴴ_trunc)...; degeneracy_atol - ) - # test_ad_rrule( - # svd_trunc, t; - # fkwargs = (; trunc), output_tangent = ΔUSVᴴ_trunc, atol, rtol - # ) + g1, = Zygote.gradient(x -> f(x; trunc), t) + g2, = Zygote.gradient(x -> f(x; trunc = truncspace(space(USVᴴ_trunc[2], 1))), t) + @test g1 ≈ g2 end end From fe4efec48d07f5348c036ad3a8c6cf15f9555f42 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 17 Nov 2025 14:14:09 -0500 Subject: [PATCH 31/41] apply code suggestions --- src/factorizations/factorizations.jl | 8 ++------ test/autodiff/ad.jl | 3 +-- 2 files changed, 3 insertions(+), 8 deletions(-) diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index 33f267e27..ad775702e 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -58,12 +58,12 @@ LinearAlgebra.svdvals!(t::AbstractTensorMap) = diagview(svd_vals!(t)) function MAK.ishermitian(t::AbstractTensorMap; kwargs...) return InnerProductStyle(t) === EuclideanInnerProduct() && domain(t) == codomain(t) && - all(cb -> MatrixAlgebraKit.ishermitian(cb[2]; kwargs...), blocks(t)) + all(cb -> MAK.ishermitian(cb[2]; kwargs...), blocks(t)) end function MAK.isantihermitian(t::AbstractTensorMap; kwargs...) return InnerProductStyle(t) === EuclideanInnerProduct() && domain(t) == codomain(t) && - all(cb -> MatrixAlgebraKit.isantihermitian(cb[2]; kwargs...), blocks(t)) + all(cb -> MAK.isantihermitian(cb[2]; kwargs...), blocks(t)) end LinearAlgebra.ishermitian(t::AbstractTensorMap) = MAK.ishermitian(t) @@ -91,9 +91,5 @@ function MAK.is_right_isometric(t::AbstractTensorMap; kwargs...) f((c, b)) = MAK.is_right_isometric(b; kwargs...) return all(f, blocks(t)) end -function MAK.isisometric(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple; kwargs...) - t = permute(t, (p₁, p₂); copy = false) - return MAK.isisometric(t; kwargs...) -end end diff --git a/test/autodiff/ad.jl b/test/autodiff/ad.jl index c3738ada7..66a2efb20 100644 --- a/test/autodiff/ad.jl +++ b/test/autodiff/ad.jl @@ -124,7 +124,6 @@ function remove_eighgauge_dependence!( ΔV, D, V; degeneracy_atol = MatrixAlgebraKit.default_pullback_degeneracy_atol(D) ) gaugepart = project_antihermitian!(V' * ΔV) - gaugepart = (gaugepart - gaugepart') / 2 for (c, b) in blocks(gaugepart) Dc = diagview(block(D, c)) # for some reason this fails only on tests, and I cannot reproduce it in an @@ -598,7 +597,7 @@ for V in spacelist # Jacobian changes size if the truncation space changes, causing errors. # So, first test the fixed space case, then do more limited testing on # some gradients and compare to the fixed space case - V_trunc = spacetype(t)(c => ceil(Int, min(size(b)...) / 2) for (c, b) in blocks(t)) + V_trunc = spacetype(t)(c => div(min(size(b)...), 2) for (c, b) in blocks(t)) trunc = truncspace(V_trunc) USVᴴ_trunc = svd_trunc(t; trunc) ΔUSVᴴ_trunc = (rand_tangent.(Base.front(USVᴴ_trunc))..., zero(last(USVᴴ_trunc))) From 4d58769f0bf9bb98270fd7b7e2e8467a198a8090 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 18 Nov 2025 09:21:33 -0500 Subject: [PATCH 32/41] stabilize AD test --- test/autodiff/ad.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/test/autodiff/ad.jl b/test/autodiff/ad.jl index 66a2efb20..53f007b21 100644 --- a/test/autodiff/ad.jl +++ b/test/autodiff/ad.jl @@ -615,23 +615,23 @@ for V in spacelist return LinearAlgebra.tr(Str) + LinearAlgebra.norm(Utr * Vᴴtr) end - trunc = truncrank(round(Int, dim(V_trunc))) - USVᴴ_trunc = svd_trunc(t; trunc) + trunc = truncrank(ceil(Int, dim(V_trunc))) + USVᴴ_trunc′ = svd_trunc(t; trunc) g1, = Zygote.gradient(x -> f(x; trunc), t) - g2, = Zygote.gradient(x -> f(x; trunc = truncspace(space(USVᴴ_trunc[2], 1))), t) + g2, = Zygote.gradient(x -> f(x; trunc = truncspace(space(USVᴴ_trunc′[2], 1))), t) @test g1 ≈ g2 trunc = truncerror(; atol = last(USVᴴ_trunc)) - USVᴴ_trunc = svd_trunc(t; trunc) + USVᴴ_trunc′ = svd_trunc(t; trunc) g1, = Zygote.gradient(x -> f(x; trunc), t) - g2, = Zygote.gradient(x -> f(x; trunc = truncspace(space(USVᴴ_trunc[2], 1))), t) + g2, = Zygote.gradient(x -> f(x; trunc = truncspace(space(USVᴴ_trunc′[2], 1))), t) @test g1 ≈ g2 - tol = minimum(((c, b),) -> minimum(diagview(b)), blocks(USVᴴ_trunc[2])) + tol = minimum(((c, b),) -> minimum(diagview(b)), blocks(USVᴴ_trunc[2]); init = zero(scalartype(USVᴴ_trunc[2]))) trunc = trunctol(; atol = 10 * tol) - USVᴴ_trunc = svd_trunc(t; trunc) + USVᴴ_trunc′ = svd_trunc(t; trunc) g1, = Zygote.gradient(x -> f(x; trunc), t) - g2, = Zygote.gradient(x -> f(x; trunc = truncspace(space(USVᴴ_trunc[2], 1))), t) + g2, = Zygote.gradient(x -> f(x; trunc = truncspace(space(USVᴴ_trunc′[2], 1))), t) @test g1 ≈ g2 end end From 74f57ee31f7970bdaaa07894d73919302cac168b Mon Sep 17 00:00:00 2001 From: Jutho Haegeman Date: Sat, 22 Nov 2025 02:05:19 +0100 Subject: [PATCH 33/41] some changes --- src/factorizations/factorizations.jl | 25 +++++++++++++------------ src/factorizations/matrixalgebrakit.jl | 23 +++++++++++------------ 2 files changed, 24 insertions(+), 24 deletions(-) diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index ad775702e..4ed891da4 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -55,15 +55,21 @@ LinearAlgebra.svdvals!(t::AbstractTensorMap) = diagview(svd_vals!(t)) #--------------------------------------------------# # Checks for hermiticity and positive definiteness # #--------------------------------------------------# +function _blockmap(f; kwargs...) + return function ((c, b)) + return f(b; kwargs...) + end +end + function MAK.ishermitian(t::AbstractTensorMap; kwargs...) return InnerProductStyle(t) === EuclideanInnerProduct() && - domain(t) == codomain(t) && - all(cb -> MAK.ishermitian(cb[2]; kwargs...), blocks(t)) + domain(t) == codomain(t) && + all(_blockmap(MAK.ishermitian; kwargs...), blocks(t)) end function MAK.isantihermitian(t::AbstractTensorMap; kwargs...) return InnerProductStyle(t) === EuclideanInnerProduct() && - domain(t) == codomain(t) && - all(cb -> MAK.isantihermitian(cb[2]; kwargs...), blocks(t)) + domain(t) == codomain(t) && + all(_blockmap(MAK.isantihermitian; kwargs...), blocks(t)) end LinearAlgebra.ishermitian(t::AbstractTensorMap) = MAK.ishermitian(t) @@ -74,22 +80,17 @@ function LinearAlgebra.isposdef!(t::AbstractTensorMap) domain(t) == codomain(t) || throw(SpaceMismatch("`isposdef` requires domain and codomain to be the same")) InnerProductStyle(spacetype(t)) === EuclideanInnerProduct() || return false - for (c, b) in blocks(t) - isposdef!(b) || return false - end - return true + return all(_blockmap(isposdef!), blocks(t)) end # TODO: tolerances are per-block, not global or weighted - does that matter? function MAK.is_left_isometric(t::AbstractTensorMap; kwargs...) domain(t) ≾ codomain(t) || return false - f((c, b)) = MAK.is_left_isometric(b; kwargs...) - return all(f, blocks(t)) + return all(_blockmap(MAK.is_left_isometric; kwargs...), blocks(t)) end function MAK.is_right_isometric(t::AbstractTensorMap; kwargs...) domain(t) ≿ codomain(t) || return false - f((c, b)) = MAK.is_right_isometric(b; kwargs...) - return all(f, blocks(t)) + return all(_blockmap(MAK.is_right_isometric; kwargs...), blocks(t)) end end diff --git a/src/factorizations/matrixalgebrakit.jl b/src/factorizations/matrixalgebrakit.jl index 72d7015c0..90c8e7645 100644 --- a/src/factorizations/matrixalgebrakit.jl +++ b/src/factorizations/matrixalgebrakit.jl @@ -31,12 +31,11 @@ for f! in ( @eval function MAK.$f!(t::AbstractTensorMap, F, alg::AbstractAlgorithm) MAK.check_input($f!, t, F, alg) - foreachblock(t, F...) do _, bs - factors = Base.tail(bs) - factors′ = $f!(first(bs), factors, alg) + foreachblock(t, F...) do _, (tblock, Fblocks...) + Fblocks′ = $f!(tblock, Fblocks, alg) # deal with the case where the output is not in-place - for (f′, f) in zip(factors′, factors) - f′ === f || copy!(f, f′) + for (b′, b) in zip(Fblocks′, Fblocks) + b === b′ || copy!(b, b′) end return nothing end @@ -50,10 +49,10 @@ for f! in (:qr_null!, :lq_null!, :project_hermitian!, :project_antihermitian!, : @eval function MAK.$f!(t::AbstractTensorMap, N, alg::AbstractAlgorithm) MAK.check_input($f!, t, N, alg) - foreachblock(t, N) do _, (b, n) - n′ = $f!(b, n, alg) + foreachblock(t, N) do _, (tblock, Nblock) + Nblock′ = $f!(tblock, Nblock, alg) # deal with the case where the output is not the same as the input - n === n′ || copy!(n, n′) + Nblock === Nblock′ || copy!(Nblock, Nblock′) return nothing end @@ -66,10 +65,10 @@ for f! in (:svd_vals!, :eig_vals!, :eigh_vals!) @eval function MAK.$f!(t::AbstractTensorMap, N, alg::AbstractAlgorithm) MAK.check_input($f!, t, N, alg) - foreachblock(t, N) do _, (b, n) - n′ = $f!(b, diagview(n), alg) + foreachblock(t, N) do _, (tblock, Nblock) + Nblock′ = $f!(tblock, diagview(Nblock), alg) # deal with the case where the output is not the same as the input - diagview(n) === n′ || copy!(diagview(n), n′) + diagview(Nblock) === Nblock′ || copy!(diagview(Nblock), Nblock′) return nothing end @@ -445,7 +444,7 @@ end # Projections # ----------- function MAK.check_input(::typeof(project_hermitian!), tsrc::AbstractTensorMap, tdst::AbstractTensorMap, ::AbstractAlgorithm) - domain(tsrc) == codomain(tsrc) || throw(ArgumentError("Hermitian projection requires square input tensor")) + domain(tsrc) == codomain(tsrc) || throw(ArgumentError("(Anti-)Hermitian projection requires square input tensor")) tsrc === tdst || @check_space(tdst, space(tsrc)) return nothing end From 260cc137def561958b2df93bbdae85bb128b89ee Mon Sep 17 00:00:00 2001 From: Jutho Haegeman Date: Sat, 22 Nov 2025 02:11:06 +0100 Subject: [PATCH 34/41] more updates: formatting and deprecations --- src/auxiliary/deprecate.jl | 2 +- src/factorizations/factorizations.jl | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/auxiliary/deprecate.jl b/src/auxiliary/deprecate.jl index a37c40531..81fbe9a50 100644 --- a/src/auxiliary/deprecate.jl +++ b/src/auxiliary/deprecate.jl @@ -186,7 +186,7 @@ function tsvd(t::AbstractTensorMap; kwargs...) Base.depwarn("p is a deprecated kwarg, and should be specified through the truncation strategy", :tsvd) kwargs = _drop_p(; kwargs...) end - return haskey(kwargs, :trunc) ? svd_trunc(t; kwargs...) : svd_compact(t; kwargs...) + return haskey(kwargs, :trunc) ? svd_trunc(t; kwargs...) : (svd_compact(t; kwargs...)..., abs(zero(scalartype(t)))) end function tsvd!(t::AbstractTensorMap; kwargs...) Base.depwarn("`tsvd!` is deprecated, use `svd_compact!`, `svd_full!` or `svd_trunc!` instead", :tsvd!) diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index 4ed891da4..44d7e2315 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -63,13 +63,13 @@ end function MAK.ishermitian(t::AbstractTensorMap; kwargs...) return InnerProductStyle(t) === EuclideanInnerProduct() && - domain(t) == codomain(t) && - all(_blockmap(MAK.ishermitian; kwargs...), blocks(t)) + domain(t) == codomain(t) && + all(_blockmap(MAK.ishermitian; kwargs...), blocks(t)) end function MAK.isantihermitian(t::AbstractTensorMap; kwargs...) return InnerProductStyle(t) === EuclideanInnerProduct() && - domain(t) == codomain(t) && - all(_blockmap(MAK.isantihermitian; kwargs...), blocks(t)) + domain(t) == codomain(t) && + all(_blockmap(MAK.isantihermitian; kwargs...), blocks(t)) end LinearAlgebra.ishermitian(t::AbstractTensorMap) = MAK.ishermitian(t) From cb390edfd96b300a432502fea41ef0ca3f1804f9 Mon Sep 17 00:00:00 2001 From: Jutho Haegeman Date: Mon, 24 Nov 2025 00:36:52 +0100 Subject: [PATCH 35/41] some more review updates --- src/factorizations/adjoint.jl | 16 ++++++++++------ src/factorizations/pullbacks.jl | 12 +++++++++--- src/factorizations/truncation.jl | 6 +++--- 3 files changed, 22 insertions(+), 12 deletions(-) diff --git a/src/factorizations/adjoint.jl b/src/factorizations/adjoint.jl index 41a0b9513..059116686 100644 --- a/src/factorizations/adjoint.jl +++ b/src/factorizations/adjoint.jl @@ -65,10 +65,14 @@ for (left_f, right_f) in zip( return reverse(adjoint.(MAK.initialize_output($left_f!, adjoint(t), _adjoint(alg)))) end - @eval MAK.$left_f!(t::AdjointTensorMap, F, alg::AbstractAlgorithm) = - $right_f!(adjoint(t), reverse(adjoint.(F)), _adjoint(alg)) - @eval MAK.$right_f!(t::AdjointTensorMap, F, alg::AbstractAlgorithm) = - $left_f!(adjoint(t), reverse(adjoint.(F)), _adjoint(alg)) + @eval function MAK.$left_f!(t::AdjointTensorMap, F, alg::AbstractAlgorithm) + F′ = $right_f!(adjoint(t), reverse(adjoint.(F)), _adjoint(alg)) + return reverse(adjoint.(F′)) + end + @eval function MAK.$right_f!(t::AdjointTensorMap, F, alg::AbstractAlgorithm) + F′ = $left_f!(adjoint(t), reverse(adjoint.(F)), _adjoint(alg)) + return reverse(adjoint.(F′)) + end end # 3-arg functions @@ -83,8 +87,8 @@ for f! in (:svd_full!, :svd_compact!, :svd_trunc!) return reverse(adjoint.(MAK.initialize_output($f!, adjoint(t), _adjoint(alg)))) end @eval function MAK.$f!(t::AdjointTensorMap, F, alg::AbstractAlgorithm) - $f!(adjoint(t), reverse(adjoint.(F)), _adjoint(alg)) - return F + F′ = $f!(adjoint(t), reverse(adjoint.(F)), _adjoint(alg)) + return reverse(adjoint.(F′)) end # disambiguate by prohibition diff --git a/src/factorizations/pullbacks.jl b/src/factorizations/pullbacks.jl index eeeec597f..4e48111a0 100644 --- a/src/factorizations/pullbacks.jl +++ b/src/factorizations/pullbacks.jl @@ -32,13 +32,19 @@ for pullback! in (:svd_pullback!, :eig_pullback!, :eigh_pullback!) Δt::AbstractTensorMap, t::AbstractTensorMap, F, ΔF, inds = _notrunc_ind(t); kwargs... ) - for (c, ind) in inds - Δb = block(Δt, c) - b = block(t, c) + foreachblock(Δt, t) do c, (Δb, b) Fc = block.(F, Ref(c)) ΔFc = block.(ΔF, Ref(c)) + ind = inds[c] MAK.$pullback!(Δb, b, Fc, ΔFc, ind; kwargs...) end + # for (c, ind) in inds + # Δb = block(Δt, c) + # b = block(t, c) + # Fc = block.(F, Ref(c)) + # ΔFc = block.(ΔF, Ref(c)) + # MAK.$pullback!(Δb, b, Fc, ΔFc, ind; kwargs...) + # end return Δt end end diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 0e9e2dda8..84f2569e0 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -21,7 +21,7 @@ Truncation strategy to keep the first values for each sector when sorted accordi such that the resulting vector space is no greater than `V`. """ function truncspace(space::ElementarySpace; by = abs, rev::Bool = true) - isdual(space) && throw(ArgumentError("resulting vector space is never dual")) + isdual(space) && throw(ArgumentError("truncation space should not be dual")) return TruncationSpace(space, by, rev) end @@ -37,7 +37,7 @@ function truncate_domain!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap, inds for (c, b) in blocks(tdst) I = get(inds, c, nothing) @assert !isnothing(I) - copy!(b, @view(block(tsrc, c)[:, I])) + copy!(b, view(block(tsrc, c), :, I)) end return tdst end @@ -45,7 +45,7 @@ function truncate_codomain!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap, in for (c, b) in blocks(tdst) I = get(inds, c, nothing) @assert !isnothing(I) - copy!(b, @view(block(tsrc, c)[I, :])) + copy!(b, view(block(tsrc, c), I, :)) end return tdst end From 0a3624c846fc077b393dd09f0fbe823652b71c02 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Fri, 28 Nov 2025 15:09:34 -0500 Subject: [PATCH 36/41] handle fully-truncated blocks --- src/factorizations/pullbacks.jl | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/src/factorizations/pullbacks.jl b/src/factorizations/pullbacks.jl index 4e48111a0..488089809 100644 --- a/src/factorizations/pullbacks.jl +++ b/src/factorizations/pullbacks.jl @@ -33,18 +33,13 @@ for pullback! in (:svd_pullback!, :eig_pullback!, :eigh_pullback!) kwargs... ) foreachblock(Δt, t) do c, (Δb, b) + haskey(inds, c) || return nothing + ind = inds[c] Fc = block.(F, Ref(c)) ΔFc = block.(ΔF, Ref(c)) - ind = inds[c] MAK.$pullback!(Δb, b, Fc, ΔFc, ind; kwargs...) + return nothing end - # for (c, ind) in inds - # Δb = block(Δt, c) - # b = block(t, c) - # Fc = block.(F, Ref(c)) - # ΔFc = block.(ΔF, Ref(c)) - # MAK.$pullback!(Δb, b, Fc, ΔFc, ind; kwargs...) - # end return Δt end end From 2e70b97277f43aa3d29d37ab8846da19975b8233 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Fri, 28 Nov 2025 15:12:27 -0500 Subject: [PATCH 37/41] fix diagonal adjoint implementation --- src/factorizations/adjoint.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/factorizations/adjoint.jl b/src/factorizations/adjoint.jl index 059116686..283bf6797 100644 --- a/src/factorizations/adjoint.jl +++ b/src/factorizations/adjoint.jl @@ -109,6 +109,7 @@ function MAK.svd_trunc!(t::AdjointTensorMap, USVᴴ, alg::TruncatedAlgorithm) USVᴴ′ = svd_compact!(t, USVᴴ, alg.alg) return MAK.truncate(svd_trunc!, USVᴴ′, alg.trunc) end -function MAK.svd_compact!(t::AdjointTensorMap, USVᴴ, alg::DiagonalAlgorithm) - return MAK.svd_compact!(t, USVᴴ, alg.alg) +function MAK.svd_compact!(t::AdjointTensorMap, F, alg::DiagonalAlgorithm) + F′ = svd_compact!(adjoint(t), reverse(adjoint.(F)), _adjoint(alg)) + return reverse(adjoint.(F′)) end From 9442d3f387a6880e70d53b0995aa5276062d7739 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 29 Nov 2025 17:16:49 -0500 Subject: [PATCH 38/41] fix property name for PolarViaSVD --- src/factorizations/adjoint.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/factorizations/adjoint.jl b/src/factorizations/adjoint.jl index 283bf6797..2f8570c1b 100644 --- a/src/factorizations/adjoint.jl +++ b/src/factorizations/adjoint.jl @@ -6,7 +6,7 @@ _adjoint(alg::MAK.LAPACK_HouseholderQR) = MAK.LAPACK_HouseholderLQ(; alg.kwargs. _adjoint(alg::MAK.LAPACK_HouseholderLQ) = MAK.LAPACK_HouseholderQR(; alg.kwargs...) _adjoint(alg::MAK.LAPACK_HouseholderQL) = MAK.LAPACK_HouseholderRQ(; alg.kwargs...) _adjoint(alg::MAK.LAPACK_HouseholderRQ) = MAK.LAPACK_HouseholderQL(; alg.kwargs...) -_adjoint(alg::MAK.PolarViaSVD) = MAK.PolarViaSVD(_adjoint(alg.svdalg)) +_adjoint(alg::MAK.PolarViaSVD) = MAK.PolarViaSVD(_adjoint(alg.svd_alg)) _adjoint(alg::AbstractAlgorithm) = alg for f in From e8de3f5e74a4e71020bbadd27fe304b1e26b78a8 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 29 Nov 2025 17:17:18 -0500 Subject: [PATCH 39/41] copy_input uses `f` instead of `f!` --- src/factorizations/adjoint.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/factorizations/adjoint.jl b/src/factorizations/adjoint.jl index 2f8570c1b..efacf73bd 100644 --- a/src/factorizations/adjoint.jl +++ b/src/factorizations/adjoint.jl @@ -47,11 +47,11 @@ for (left_f, right_f) in zip( ) left_f! = Symbol(left_f, :!) right_f! = Symbol(right_f, :!) - @eval function MAK.copy_input(::typeof($left_f!), t::AdjointTensorMap) - return adjoint(MAK.copy_input($right_f!, adjoint(t))) + @eval function MAK.copy_input(::typeof($left_f), t::AdjointTensorMap) + return adjoint(MAK.copy_input($right_f, adjoint(t))) end - @eval function MAK.copy_input(::typeof($right_f!), t::AdjointTensorMap) - return adjoint(MAK.copy_input($left_f!, adjoint(t))) + @eval function MAK.copy_input(::typeof($right_f), t::AdjointTensorMap) + return adjoint(MAK.copy_input($left_f, adjoint(t))) end @eval function MAK.initialize_output( From 68beac78abfd61517498a98417218161a2eb639c Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 29 Nov 2025 17:20:24 -0500 Subject: [PATCH 40/41] remove more implementations to fix code --- src/factorizations/adjoint.jl | 27 +++++++++++--------------- src/factorizations/matrixalgebrakit.jl | 8 +++++--- 2 files changed, 16 insertions(+), 19 deletions(-) diff --git a/src/factorizations/adjoint.jl b/src/factorizations/adjoint.jl index efacf73bd..29f27203c 100644 --- a/src/factorizations/adjoint.jl +++ b/src/factorizations/adjoint.jl @@ -11,9 +11,11 @@ _adjoint(alg::AbstractAlgorithm) = alg for f in [ - :svd_compact, :svd_full, :svd_trunc, :svd_vals, :qr_compact, :qr_full, :qr_null, - :lq_compact, :lq_full, :lq_null, :eig_full, :eig_trunc, :eig_vals, :eigh_full, - :eigh_trunc, :eigh_vals, :left_polar, :right_polar, + :svd_compact, :svd_full, :svd_vals, + :qr_compact, :qr_full, :qr_null, + :lq_compact, :lq_full, :lq_null, + :eig_full, :eig_vals, :eigh_full, :eigh_trunc, :eigh_vals, + :left_polar, :right_polar, :project_hermitian, :project_antihermitian, :project_isometric, ] f! = Symbol(f, :!) @@ -76,9 +78,10 @@ for (left_f, right_f) in zip( end # 3-arg functions -for f! in (:svd_full!, :svd_compact!, :svd_trunc!) - @eval function MAK.copy_input(::typeof($f!), t::AdjointTensorMap) - return adjoint(MAK.copy_input($f!, adjoint(t))) +for f in (:svd_full, :svd_compact) + f! = Symbol(f, :!) + @eval function MAK.copy_input(::typeof($f), t::AdjointTensorMap) + return adjoint(MAK.copy_input($f, adjoint(t))) end @eval function MAK.initialize_output( @@ -86,6 +89,7 @@ for f! in (:svd_full!, :svd_compact!, :svd_trunc!) ) return reverse(adjoint.(MAK.initialize_output($f!, adjoint(t), _adjoint(alg)))) end + @eval function MAK.$f!(t::AdjointTensorMap, F, alg::AbstractAlgorithm) F′ = $f!(adjoint(t), reverse(adjoint.(F)), _adjoint(alg)) return reverse(adjoint.(F′)) @@ -98,17 +102,8 @@ for f! in (:svd_full!, :svd_compact!, :svd_trunc!) throw(MethodError($f!, (t, alg))) end end + # avoid amgiguity -function MAK.initialize_output( - ::typeof(svd_trunc!), t::AdjointTensorMap, alg::TruncatedAlgorithm - ) - return MAK.initialize_output(svd_compact!, t, alg.alg) -end -# to fix ambiguity -function MAK.svd_trunc!(t::AdjointTensorMap, USVᴴ, alg::TruncatedAlgorithm) - USVᴴ′ = svd_compact!(t, USVᴴ, alg.alg) - return MAK.truncate(svd_trunc!, USVᴴ′, alg.trunc) -end function MAK.svd_compact!(t::AdjointTensorMap, F, alg::DiagonalAlgorithm) F′ = svd_compact!(adjoint(t), reverse(adjoint.(F)), _adjoint(alg)) return reverse(adjoint.(F′)) diff --git a/src/factorizations/matrixalgebrakit.jl b/src/factorizations/matrixalgebrakit.jl index 90c8e7645..a1da88466 100644 --- a/src/factorizations/matrixalgebrakit.jl +++ b/src/factorizations/matrixalgebrakit.jl @@ -2,9 +2,11 @@ # ------------------- for f in [ - :svd_compact, :svd_full, :svd_trunc, :svd_vals, :qr_compact, :qr_full, :qr_null, - :lq_compact, :lq_full, :lq_null, :eig_full, :eig_trunc, :eig_vals, :eigh_full, - :eigh_trunc, :eigh_vals, :left_polar, :right_polar, + :svd_compact, :svd_full, :svd_vals, + :qr_compact, :qr_full, :qr_null, + :lq_compact, :lq_full, :lq_null, + :eig_full, :eig_vals, :eigh_full, :eigh_vals, + :left_polar, :right_polar, :project_hermitian, :project_antihermitian, :project_isometric, ] f! = Symbol(f, :!) From ea7cbc98cbdf2833e70a7c9b4ac90e10c038b88b Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 29 Nov 2025 17:50:40 -0500 Subject: [PATCH 41/41] remove `check_input` --- src/factorizations/diagonal.jl | 96 --------- src/factorizations/matrixalgebrakit.jl | 287 +------------------------ src/factorizations/utility.jl | 8 - 3 files changed, 6 insertions(+), 385 deletions(-) diff --git a/src/factorizations/diagonal.jl b/src/factorizations/diagonal.jl index 4de99f973..cea915995 100644 --- a/src/factorizations/diagonal.jl +++ b/src/factorizations/diagonal.jl @@ -81,44 +81,11 @@ for f! in :eigh_trunc!, :right_orth!, :left_orth!, ) @eval function MAK.$f!(d::DiagonalTensorMap, F, alg::DiagonalAlgorithm) - MAK.check_input($f!, d, F, alg) $f!(_repack_diagonal(d), _repack_diagonal.(F), alg) return F end end -for f! in (:qr_full!, :qr_compact!) - @eval function MAK.check_input( - ::typeof($f!), d::AbstractTensorMap, QR, ::DiagonalAlgorithm - ) - Q, R = QR - @assert d isa DiagonalTensorMap - @assert Q isa DiagonalTensorMap && R isa DiagonalTensorMap - @check_scalar Q d - @check_scalar R d - @check_space(Q, space(d)) - @check_space(R, space(d)) - - return nothing - end -end - -for f! in (:lq_full!, :lq_compact!) - @eval function MAK.check_input( - ::typeof($f!), d::AbstractTensorMap, LQ, ::DiagonalAlgorithm - ) - L, Q = LQ - @assert d isa DiagonalTensorMap - @assert Q isa DiagonalTensorMap && L isa DiagonalTensorMap - @check_scalar Q d - @check_scalar L d - @check_space(Q, space(d)) - @check_space(L, space(d)) - - return nothing - end -end - # disambiguate function MAK.svd_compact!(t::AbstractTensorMap, USVᴴ, alg::DiagonalAlgorithm) return svd_full!(t, USVᴴ, alg) @@ -126,10 +93,8 @@ end # f_vals # ------ - for f! in (:eig_vals!, :eigh_vals!, :svd_vals!) @eval function MAK.$f!(d::AbstractTensorMap, V, alg::DiagonalAlgorithm) - MAK.check_input($f!, d, V, alg) $f!(_repack_diagonal(d), diagview(_repack_diagonal(V)), alg) return V end @@ -140,64 +105,3 @@ for f! in (:eig_vals!, :eigh_vals!, :svd_vals!) return DiagonalTensorMap(data, d.domain) end end - -function MAK.check_input(::typeof(eig_full!), t::AbstractTensorMap, DV, ::DiagonalAlgorithm) - domain(t) == codomain(t) || - throw(ArgumentError("Eigenvalue decomposition requires square input tensor")) - - D, V = DV - - @assert D isa DiagonalTensorMap - @assert V isa AbstractTensorMap - - # scalartype checks - @check_scalar D t - @check_scalar V t - - # space checks - @check_space D space(t) - @check_space V space(t) - - return nothing -end - -function MAK.check_input(::typeof(eigh_full!), t::AbstractTensorMap, DV, ::DiagonalAlgorithm) - domain(t) == codomain(t) || - throw(ArgumentError("Eigenvalue decomposition requires square input tensor")) - - D, V = DV - - @assert D isa DiagonalTensorMap - @assert V isa AbstractTensorMap - - # scalartype checks - @check_scalar D t real - @check_scalar V t - - # space checks - @check_space D space(t) - @check_space V space(t) - - return nothing -end - -function MAK.check_input(::typeof(eig_vals!), t::AbstractTensorMap, D, ::DiagonalAlgorithm) - @assert D isa DiagonalTensorMap - @check_scalar D t - @check_space D space(t) - return nothing -end - -function MAK.check_input(::typeof(eigh_vals!), t::AbstractTensorMap, D, ::DiagonalAlgorithm) - @assert D isa DiagonalTensorMap - @check_scalar D t real - @check_space D space(t) - return nothing -end - -function MAK.check_input(::typeof(svd_vals!), t::AbstractTensorMap, D, ::DiagonalAlgorithm) - @assert D isa DiagonalTensorMap - @check_scalar D t real - @check_space D space(t) - return nothing -end diff --git a/src/factorizations/matrixalgebrakit.jl b/src/factorizations/matrixalgebrakit.jl index a1da88466..412fd8761 100644 --- a/src/factorizations/matrixalgebrakit.jl +++ b/src/factorizations/matrixalgebrakit.jl @@ -31,8 +31,6 @@ for f! in ( :left_polar!, :right_polar!, ) @eval function MAK.$f!(t::AbstractTensorMap, F, alg::AbstractAlgorithm) - MAK.check_input($f!, t, F, alg) - foreachblock(t, F...) do _, (tblock, Fblocks...) Fblocks′ = $f!(tblock, Fblocks, alg) # deal with the case where the output is not in-place @@ -41,7 +39,6 @@ for f! in ( end return nothing end - return F end end @@ -49,15 +46,12 @@ end # Handle these separately because single output instead of tuple 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) - foreachblock(t, N) do _, (tblock, Nblock) Nblock′ = $f!(tblock, Nblock, alg) # deal with the case where the output is not the same as the input Nblock === Nblock′ || copy!(Nblock, Nblock′) return nothing end - return N end end @@ -65,74 +59,18 @@ end # Handle these separately because single output instead of tuple for f! in (:svd_vals!, :eig_vals!, :eigh_vals!) @eval function MAK.$f!(t::AbstractTensorMap, N, alg::AbstractAlgorithm) - MAK.check_input($f!, t, N, alg) - foreachblock(t, N) do _, (tblock, Nblock) Nblock′ = $f!(tblock, diagview(Nblock), alg) # deal with the case where the output is not the same as the input diagview(Nblock) === Nblock′ || copy!(diagview(Nblock), Nblock′) return nothing end - return N end end # Singular value decomposition # ---------------------------- -function MAK.check_input(::typeof(svd_full!), t::AbstractTensorMap, USVᴴ, ::AbstractAlgorithm) - U, S, Vᴴ = USVᴴ - - # type checks - @assert U isa AbstractTensorMap - @assert S isa AbstractTensorMap - @assert Vᴴ isa AbstractTensorMap - - # scalartype checks - @check_scalar U t - @check_scalar S t real - @check_scalar Vᴴ t - - # space checks - V_cod = fuse(codomain(t)) - V_dom = fuse(domain(t)) - @check_space(U, codomain(t) ← V_cod) - @check_space(S, V_cod ← V_dom) - @check_space(Vᴴ, V_dom ← domain(t)) - - return nothing -end - -function MAK.check_input(::typeof(svd_compact!), t::AbstractTensorMap, USVᴴ, ::AbstractAlgorithm) - U, S, Vᴴ = USVᴴ - - # type checks - @assert U isa AbstractTensorMap - @assert S isa DiagonalTensorMap - @assert Vᴴ isa AbstractTensorMap - - # scalartype checks - @check_scalar U t - @check_scalar S t real - @check_scalar Vᴴ t - - # space checks - V_cod = V_dom = infimum(fuse(codomain(t)), fuse(domain(t))) - @check_space(U, codomain(t) ← V_cod) - @check_space(S, V_cod ← V_dom) - @check_space(Vᴴ, V_dom ← domain(t)) - - return nothing -end - -function MAK.check_input(::typeof(svd_vals!), t::AbstractTensorMap, D, ::AbstractAlgorithm) - @check_scalar D t real - @assert D isa DiagonalTensorMap - V_cod = V_dom = infimum(fuse(codomain(t)), fuse(domain(t))) - @check_space(D, V_cod ← V_dom) - return nothing -end - function MAK.initialize_output(::typeof(svd_full!), t::AbstractTensorMap, ::AbstractAlgorithm) V_cod = fuse(codomain(t)) V_dom = fuse(domain(t)) @@ -157,66 +95,6 @@ end # Eigenvalue decomposition # ------------------------ -function MAK.check_input(::typeof(eigh_full!), t::AbstractTensorMap, DV, ::AbstractAlgorithm) - domain(t) == codomain(t) || - throw(ArgumentError("Eigenvalue decomposition requires square input tensor")) - - D, V = DV - - # type checks - @assert D isa DiagonalTensorMap - @assert V isa AbstractTensorMap - - # scalartype checks - @check_scalar D t real - @check_scalar V t - - # space checks - V_D = fuse(domain(t)) - @check_space(D, V_D ← V_D) - @check_space(V, codomain(t) ← V_D) - - return nothing -end - -function MAK.check_input(::typeof(eig_full!), t::AbstractTensorMap, DV, ::AbstractAlgorithm) - domain(t) == codomain(t) || - throw(ArgumentError("Eigenvalue decomposition requires square input tensor")) - - D, V = DV - - # type checks - @assert D isa DiagonalTensorMap - @assert V isa AbstractTensorMap - - # scalartype checks - @check_scalar D t complex - @check_scalar V t complex - - # space checks - V_D = fuse(domain(t)) - @check_space(D, V_D ← V_D) - @check_space(V, codomain(t) ← V_D) - - return nothing -end - -function MAK.check_input(::typeof(eigh_vals!), t::AbstractTensorMap, D, ::AbstractAlgorithm) - @check_scalar D t real - @assert D isa DiagonalTensorMap - V_D = fuse(domain(t)) - @check_space(D, V_D ← V_D) - return nothing -end - -function MAK.check_input(::typeof(eig_vals!), t::AbstractTensorMap, D, ::AbstractAlgorithm) - @check_scalar D t complex - @assert D isa DiagonalTensorMap - V_D = fuse(domain(t)) - @check_space(D, V_D ← V_D) - return nothing -end - function MAK.initialize_output(::typeof(eigh_full!), t::AbstractTensorMap, ::AbstractAlgorithm) V_D = fuse(domain(t)) T = real(scalartype(t)) @@ -247,56 +125,6 @@ end # QR decomposition # ---------------- -function MAK.check_input(::typeof(qr_full!), t::AbstractTensorMap, QR, ::AbstractAlgorithm) - Q, R = QR - - # type checks - @assert Q isa AbstractTensorMap - @assert R isa AbstractTensorMap - - # scalartype checks - @check_scalar Q t - @check_scalar R t - - # space checks - V_Q = fuse(codomain(t)) - @check_space(Q, codomain(t) ← V_Q) - @check_space(R, V_Q ← domain(t)) - - return nothing -end - -function MAK.check_input(::typeof(qr_compact!), t::AbstractTensorMap, QR, ::AbstractAlgorithm) - Q, R = QR - - # type checks - @assert Q isa AbstractTensorMap - @assert R isa AbstractTensorMap - - # scalartype checks - @check_scalar Q t - @check_scalar R t - - # space checks - V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) - @check_space(Q, codomain(t) ← V_Q) - @check_space(R, V_Q ← domain(t)) - - return nothing -end - -function MAK.check_input(::typeof(qr_null!), t::AbstractTensorMap, N, ::AbstractAlgorithm) - # scalartype checks - @check_scalar N t - - # space checks - V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) - V_N = ⊖(fuse(codomain(t)), V_Q) - @check_space(N, codomain(t) ← V_N) - - return nothing -end - function MAK.initialize_output(::typeof(qr_full!), t::AbstractTensorMap, ::AbstractAlgorithm) V_Q = fuse(codomain(t)) Q = similar(t, codomain(t) ← V_Q) @@ -320,56 +148,6 @@ end # LQ decomposition # ---------------- -function MAK.check_input(::typeof(lq_full!), t::AbstractTensorMap, LQ, ::AbstractAlgorithm) - L, Q = LQ - - # type checks - @assert L isa AbstractTensorMap - @assert Q isa AbstractTensorMap - - # scalartype checks - @check_scalar L t - @check_scalar Q t - - # space checks - V_Q = fuse(domain(t)) - @check_space(L, codomain(t) ← V_Q) - @check_space(Q, V_Q ← domain(t)) - - return nothing -end - -function MAK.check_input(::typeof(lq_compact!), t::AbstractTensorMap, LQ, ::AbstractAlgorithm) - L, Q = LQ - - # type checks - @assert L isa AbstractTensorMap - @assert Q isa AbstractTensorMap - - # scalartype checks - @check_scalar L t - @check_scalar Q t - - # space checks - V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) - @check_space(L, codomain(t) ← V_Q) - @check_space(Q, V_Q ← domain(t)) - - return nothing -end - -function MAK.check_input(::typeof(lq_null!), t::AbstractTensorMap, N, ::AbstractAlgorithm) - # scalartype checks - @check_scalar N t - - # space checks - V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) - V_N = ⊖(fuse(domain(t)), V_Q) - @check_space(N, V_N ← domain(t)) - - return nothing -end - function MAK.initialize_output(::typeof(lq_full!), t::AbstractTensorMap, ::AbstractAlgorithm) V_Q = fuse(domain(t)) L = similar(t, codomain(t) ← V_Q) @@ -393,50 +171,12 @@ end # Polar decomposition # ------------------- -function MAK.check_input(::typeof(left_polar!), t::AbstractTensorMap, WP, ::AbstractAlgorithm) - codomain(t) ≿ domain(t) || - throw(ArgumentError("Polar decomposition requires `codomain(t) ≿ domain(t)`")) - - W, P = WP - @assert W isa AbstractTensorMap - @assert P isa AbstractTensorMap - - # scalartype checks - @check_scalar W t - @check_scalar P t - - # space checks - @check_space(W, space(t)) - @check_space(P, domain(t) ← domain(t)) - - return nothing -end - function MAK.initialize_output(::typeof(left_polar!), t::AbstractTensorMap, ::AbstractAlgorithm) W = similar(t, space(t)) P = similar(t, domain(t) ← domain(t)) return W, P end -function MAK.check_input(::typeof(right_polar!), t::AbstractTensorMap, PWᴴ, ::AbstractAlgorithm) - codomain(t) ≾ domain(t) || - throw(ArgumentError("Polar decomposition requires `domain(t) ≿ codomain(t)`")) - - P, Wᴴ = PWᴴ - @assert P isa AbstractTensorMap - @assert Wᴴ isa AbstractTensorMap - - # scalartype checks - @check_scalar P t - @check_scalar Wᴴ t - - # space checks - @check_space(P, codomain(t) ← codomain(t)) - @check_space(Wᴴ, space(t)) - - return nothing -end - function MAK.initialize_output(::typeof(right_polar!), t::AbstractTensorMap, ::AbstractAlgorithm) P = similar(t, codomain(t) ← codomain(t)) Wᴴ = similar(t, space(t)) @@ -445,24 +185,9 @@ end # Projections # ----------- -function MAK.check_input(::typeof(project_hermitian!), tsrc::AbstractTensorMap, tdst::AbstractTensorMap, ::AbstractAlgorithm) - domain(tsrc) == codomain(tsrc) || throw(ArgumentError("(Anti-)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, alg::AbstractAlgorithm) = - MAK.check_input(project_hermitian!, tsrc, tdst, alg) - -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) +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) diff --git a/src/factorizations/utility.jl b/src/factorizations/utility.jl index 0d2c3575b..4d5e1bb00 100644 --- a/src/factorizations/utility.jl +++ b/src/factorizations/utility.jl @@ -1,11 +1,3 @@ -# convenience to set default -macro check_space(x, V) - return esc(:($MatrixAlgebraKit.@check_size($x, $V, $space))) -end -macro check_scalar(x, y, op = :identity, eltype = :scalartype) - return esc(:($MatrixAlgebraKit.@check_scalar($x, $y, $op, $eltype))) -end - function factorisation_scalartype(t::AbstractTensorMap) T = scalartype(t) return promote_type(Float32, typeof(zero(T) / sqrt(abs2(one(T)))))