From 9c9b9dd4104d661968ce28b9fbf08107a75d5a2a Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Fri, 3 Oct 2025 16:11:13 -0400 Subject: [PATCH 1/3] Adapt to v0.15 --- Project.toml | 6 ++-- src/TensorKitManifolds.jl | 61 ++++++++++----------------------------- src/auxiliary.jl | 51 ++++++++++++++++---------------- src/grassmann.jl | 4 +-- src/unitary.jl | 7 +++-- test/runtests.jl | 4 +-- 6 files changed, 52 insertions(+), 81 deletions(-) diff --git a/Project.toml b/Project.toml index 0331ca1..a363a81 100644 --- a/Project.toml +++ b/Project.toml @@ -1,14 +1,16 @@ name = "TensorKitManifolds" uuid = "11fa318c-39cb-4a83-b1ed-cdc7ba1e3684" authors = ["Jutho Haegeman ", "Markus Hauru "] -version = "0.7.2" +version = "0.7.3" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4" TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" [compat] -TensorKit = "0.13,0.14" +MatrixAlgebraKit = "0.5.0" +TensorKit = "0.15" julia = "1.10" [extras] diff --git a/src/TensorKitManifolds.jl b/src/TensorKitManifolds.jl index 900b784..b2f31ea 100644 --- a/src/TensorKitManifolds.jl +++ b/src/TensorKitManifolds.jl @@ -7,6 +7,9 @@ export Grassmann, Stiefel, Unitary export inner, retract, transport, transport! using TensorKit +using MatrixAlgebraKit: MatrixAlgebraKit, AbstractAlgorithm, Algorithm, PolarViaSVD, + LAPACK_DivideAndConquer, diagview +import MatrixAlgebraKit as MAK # Every submodule -- Grassmann, Stiefel, and Unitary -- implements their own methods for # these. The signatures should be @@ -28,23 +31,8 @@ checkbase(x, y, z, args...) = checkbase(checkbase(x, y), z, args...) # the machine epsilon for the elements of an object X, name inspired from eltype scalareps(X) = eps(real(scalartype(X))) -# default SVD algorithm used in the algorithms -default_svd_alg(::AbstractTensorMap) = TensorKit.SVD() - -function isisometry(W::AbstractTensorMap; tol=10 * scalareps(W)) - WdW = W' * W - s = zero(float(real(scalartype(W)))) - for (c, b) in blocks(WdW) - _subtractone!(b) - s += dim(c) * length(b) - end - return norm(WdW) <= tol * sqrt(s) -end - -function isunitary(W::AbstractTensorMap; tol=10 * scalareps(W)) - return isisometry(W; tol=tol) && isisometry(W'; tol=tol) -end - +# TODO: these functions should be replaced by MAK functions +projecthermitian(W::AbstractTensorMap) = projecthermitian!(copy(W)) function projecthermitian!(W::AbstractTensorMap) codomain(W) == domain(W) || throw(DomainError("Tensor with distinct domain and codomain cannot be hermitian.")) @@ -53,6 +41,8 @@ function projecthermitian!(W::AbstractTensorMap) end return W end + +projectantihermitian(W::AbstractTensorMap) = projectantihermitian!(copy(W)) function projectantihermitian!(W::AbstractTensorMap) codomain(W) == domain(W) || throw(DomainError("Tensor with distinct domain and codomain cannot be anithermitian.")) @@ -62,27 +52,18 @@ function projectantihermitian!(W::AbstractTensorMap) return W end -struct PolarNewton <: TensorKit.OrthogonalFactorizationAlgorithm -end -function projectisometric!(W::AbstractTensorMap; alg=default_svd_alg(W)) - if alg isa TensorKit.Polar || alg isa TensorKit.SDD - foreach(blocks(W)) do (c, b) - return _polarsdd!(b) - end - elseif alg isa TensorKit.SVD - foreach(blocks(W)) do (c, b) - return _polarsvd!(b) - end - elseif alg isa PolarNewton - foreach(blocks(W)) do (c, b) - return _polarnewton!(b) - end - else - throw(ArgumentError("unkown algorithm for projectisometric!: alg = $alg")) +projectisometric(W::AbstractTensorMap; kwargs...) = projectisometric!(copy(W); kwargs...) +function projectisometric!(W::AbstractTensorMap; + alg::AbstractAlgorithm=MAK.select_algorithm(left_polar!, W)) + TensorKit.foreachblock(W) do c, (b,) + return _left_polar!(b, alg) end return W end +function projectcomplement(X::AbstractTensorMap, W::AbstractTensorMap, kwargs...) + return projectcomplement!(copy(X), W; kwargs...) +end function projectcomplement!(X::AbstractTensorMap, W::AbstractTensorMap; tol=10 * scalareps(X)) P = W' * X @@ -97,18 +78,6 @@ function projectcomplement!(X::AbstractTensorMap, W::AbstractTensorMap; return X end -projecthermitian(W::AbstractTensorMap) = projecthermitian!(copy(W)) -projectantihermitian(W::AbstractTensorMap) = projectantihermitian!(copy(W)) - -function projectisometric(W::AbstractTensorMap; - alg=default_svd_alg(W)) - return projectisometric!(copy(W); alg=alg) -end -function projectcomplement(X::AbstractTensorMap, W::AbstractTensorMap, - tol=10 * scalareps(X)) - return projectcomplement!(copy(X), W; tol=tol) -end - include("auxiliary.jl") include("grassmann.jl") include("stiefel.jl") diff --git a/src/auxiliary.jl b/src/auxiliary.jl index d5af541..8fea4da 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -70,20 +70,24 @@ function _subtractone!(a::AbstractMatrix) view(a, diagind(a)) .= view(a, diagind(a)) .- 1 return a end -function _polarsdd!(A::StridedMatrix) - U, S, V = svd!(A; alg=LinearAlgebra.DivideAndConquer()) - return mul!(A, U, V') -end -function _polarsvd!(A::StridedMatrix) - U, S, V = svd!(A; alg=LinearAlgebra.QRIteration()) - return mul!(A, U, V') + +# TODO: _left_polar! is more or less the same as MAK.left_polar! but doesn't compute the P +# which is not needed here. Can we unify this? +function _left_polar!(A::StridedMatrix, alg::PolarViaSVD=PolarViaSVD(LAPACK_DivideAndConquer())) + U, _, Vᴴ = svd_compact!(A, alg.svdalg) + return mul!(A, U, Vᴴ) end + +# TODO: can we move this to a dedicated MAK algorithm? +MatrixAlgebraKit.@algdef PolarNewton + +_left_polar!(A::StridedMatrix, alg::PolarNewton) = _polarnewton!(A; alg.kwargs...) function _polarnewton!(A::StridedMatrix; tol=10 * scalareps(A), maxiter=5) m, n = size(A) @assert m >= n A2 = copy(A) - Q, R = qr!(A2) - Ri = ldiv!(UpperTriangular(R)', TensorKit.MatrixAlgebra.one!(similar(R))) + Q, R = LinearAlgebra.qr!(A2) + Ri = ldiv!(UpperTriangular(R)', MatrixAlgebraKit.one!(similar(R))) R, Ri = _avgdiff!(R, Ri) i = 1 R2 = view(A, 1:n, 1:n) @@ -91,16 +95,17 @@ function _polarnewton!(A::StridedMatrix; tol=10 * scalareps(A), maxiter=5) copyto!(R2, R) while maximum(abs, Ri) > tol if i == maxiter # if not converged by now, fall back to sdd - _polarsdd!(Ri) + _left_polar!(Ri) break end - Ri = ldiv!(lu!(R2)', TensorKit.MatrixAlgebra.one!(Ri)) + Ri = ldiv!(lu!(R2)', MatrixAlgebraKit.one!(Ri)) R, Ri = _avgdiff!(R, Ri) copyto!(R2, R) i += 1 end return lmul!(Q, A) end + # in place computation of the average and difference of two arrays function _avgdiff!(A::AbstractArray, B::AbstractArray) axes(A) == axes(B) || throw(DimensionMismatch()) @@ -124,7 +129,7 @@ end function _stiefelexp(W::StridedMatrix, A::StridedMatrix, Z::StridedMatrix, α) n, p = size(W) r = min(2 * p, n) - QQ, _ = qr!([W Z]) + QQ, _ = LinearAlgebra.qr!([W Z]) Q = similar(W, n, r - p) @inbounds for j in Base.OneTo(r - p) for i in Base.OneTo(n) @@ -139,7 +144,7 @@ function _stiefelexp(W::StridedMatrix, A::StridedMatrix, Z::StridedMatrix, α) A2[1:p, (p + 1):end] .= (-α) .* (R') A2[(p + 1):end, (p + 1):end] .= 0 U = [W Q] * exp(A2) - U = _polarnewton!(U) + U = _left_polar!(U, PolarNewton()) W′ = U[:, 1:p] Q′ = U[:, (p + 1):end] R′ = R @@ -152,7 +157,7 @@ function _stiefellog(Wold::StridedMatrix, Wnew::StridedMatrix; r = min(2 * p, n) P = Wold' * Wnew dW = Wnew - Wold * P - QQ, _ = qr!([Wold dW]) + QQ, _ = LinearAlgebra.qr!([Wold dW]) Q = similar(Wold, n, r - p) @inbounds for j in Base.OneTo(r - p) for i in Base.OneTo(n) @@ -161,23 +166,17 @@ function _stiefellog(Wold::StridedMatrix, Wnew::StridedMatrix; end Q = lmul!(QQ, Q) R = Q' * dW - Wext = [Wold Q] - F = qr!([P; R]) - U = lmul!(F.Q, TensorKit.MatrixAlgebra.one!(similar(P, r, r))) + F = LinearAlgebra.qr!([P; R]) + U = lmul!(F.Q, MatrixAlgebraKit.one!(similar(P, r, r))) U[1:p, 1:p] .= P U[(p + 1):r, 1:p] .= R X = view(U, 1:p, (p + 1):r) Y = view(U, (p + 1):r, (p + 1):r) if p < n - YSVD = svd!(Y) - mul!(X, X * (YSVD.V), (YSVD.U)') - UsqrtS = YSVD.U - @inbounds for j in 1:size(UsqrtS, 2) - s = sqrt(YSVD.S[j]) - @simd for i in 1:size(UsqrtS, 1) - UsqrtS[i, j] *= s - end - end + USVᴴ = svd_compact!(Y) + mul!(X, X * USVᴴ[3]', USVᴴ[1]') + diagview(USVᴴ[2]) .= sqrt.(diagview(USVᴴ[2])) + UsqrtS = rmul!(USVᴴ[1], USVᴴ[2]) mul!(Y, UsqrtS, UsqrtS') end logU = _projectantihermitian!(log(U)) diff --git a/src/grassmann.jl b/src/grassmann.jl index ec0c938..05d8c04 100644 --- a/src/grassmann.jl +++ b/src/grassmann.jl @@ -60,7 +60,7 @@ function Base.getproperty(Δ::GrassmannTangent, sym::Symbol) elseif sym ∈ (:U, :S, :V) v = Base.getfield(Δ, sym) v !== nothing && return v - U, S, V, = tsvd(Δ.Z; alg=default_svd_alg(Δ.Z)) + U, S, V, = svd_compact(Δ.Z) Base.setfield!(Δ, :U, U) Base.setfield!(Δ, :S, S) Base.setfield!(Δ, :V, V) @@ -198,7 +198,7 @@ function invretract(Wold::AbstractTensorMap, Wnew::AbstractTensorMap; alg=nothin space(Wold) == space(Wnew) || throw(SpaceMismatch()) WodWn = Wold' * Wnew # V' * cos(S) * V * Y Wneworth = Wnew - Wold * WodWn - Vd, cS, VY = tsvd!(WodWn; alg=default_svd_alg(WodWn)) + Vd, cS, VY = svd_compact!(WodWn) Scmplx = acos(cS) # acos always returns a complex TensorMap. We cast back to real if possible. S = scalartype(WodWn) <: Real && isreal(sectortype(Scmplx)) ? real(Scmplx) : Scmplx diff --git a/src/unitary.jl b/src/unitary.jl index a5b8c2a..1e59c5a 100644 --- a/src/unitary.jl +++ b/src/unitary.jl @@ -7,6 +7,7 @@ using TensorKit import TensorKit: similarstoragetype, SectorDict using ..TensorKitManifolds: projectantihermitian!, projectisometric!, PolarNewton import ..TensorKitManifolds: base, checkbase, inner, retract, transport, transport! +import MatrixAlgebraKit as MAK struct UnitaryTangent{T<:AbstractTensorMap,TA<:AbstractTensorMap} W::T @@ -82,10 +83,10 @@ end project(X, W; metric=:euclidean) = project!(copy(X), W; metric=:euclidean) # geodesic retraction, coincides with Stiefel retraction (which is not geodesic for p < n) -function retract(W::AbstractTensorMap, Δ::UnitaryTangent, α; alg=nothing) +function retract(W::AbstractTensorMap, Δ::UnitaryTangent, α; alg=MAK.select_algorithm(left_polar!, W)) W == base(Δ) || throw(ArgumentError("not a valid tangent vector at base point")) E = exp(α * Δ.A) - W′ = projectisometric!(W * E; alg=SDD()) + W′ = projectisometric!(W * E; alg) A′ = Δ.A return W′, UnitaryTangent(W′, A′) end @@ -104,7 +105,7 @@ function transport!(Θ::UnitaryTangent, W::AbstractTensorMap, Δ::UnitaryTangent end function transport(Θ::UnitaryTangent, W::AbstractTensorMap, Δ::UnitaryTangent, α::Real, W′; alg=:stiefel) - return transport!(copy(Θ), W, Δ, α, W′; alg=alg) + return transport!(copy(Θ), W, Δ, α, W′; alg) end # transport_parallel correspondings to the torsion-free Levi-Civita connection diff --git a/test/runtests.jl b/test/runtests.jl index f9fd9c2..3001266 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,7 +8,7 @@ const α = 0.75 @testset "Grassmann with space $V" for V in spaces for T in (Float64,) - W, = leftorth(randn(T, V * V * V, V * V); alg=Polar()) + W, = left_polar(randn(T, V * V * V, V * V)) X = randn(T, space(W)) Y = randn(T, space(W)) Δ = @inferred Grassmann.project(X, W) @@ -124,7 +124,7 @@ end @testset "Unitary with space $V" for V in spaces for T in (Float64, ComplexF64) - W, = leftorth(randn(T, V * V * V, V * V); alg=Polar()) + W, = left_polar(randn(T, V * V * V, V * V)) X = randn(T, space(W)) Y = randn(T, space(W)) Δ = @inferred Unitary.project(X, W) From f3e76878ac2b79417973630dbb09d49d89f001a6 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Fri, 3 Oct 2025 18:46:07 -0400 Subject: [PATCH 2/3] format and bugfix --- src/auxiliary.jl | 3 ++- src/grassmann.jl | 3 +-- src/unitary.jl | 3 ++- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/auxiliary.jl b/src/auxiliary.jl index 8fea4da..3552553 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -73,7 +73,8 @@ end # TODO: _left_polar! is more or less the same as MAK.left_polar! but doesn't compute the P # which is not needed here. Can we unify this? -function _left_polar!(A::StridedMatrix, alg::PolarViaSVD=PolarViaSVD(LAPACK_DivideAndConquer())) +function _left_polar!(A::StridedMatrix, + alg::PolarViaSVD=PolarViaSVD(LAPACK_DivideAndConquer())) U, _, Vᴴ = svd_compact!(A, alg.svdalg) return mul!(A, U, Vᴴ) end diff --git a/src/grassmann.jl b/src/grassmann.jl index 05d8c04..8e9532e 100644 --- a/src/grassmann.jl +++ b/src/grassmann.jl @@ -7,8 +7,7 @@ module Grassmann using TensorKit using TensorKit: similarstoragetype, SectorDict using ..TensorKitManifolds: projecthermitian!, projectantihermitian!, - projectisometric!, projectcomplement!, PolarNewton, - default_svd_alg + projectisometric!, projectcomplement!, PolarNewton import ..TensorKitManifolds: base, checkbase, inner, retract, transport, transport! # special type to store tangent vectors using Z diff --git a/src/unitary.jl b/src/unitary.jl index 1e59c5a..06647b5 100644 --- a/src/unitary.jl +++ b/src/unitary.jl @@ -83,7 +83,8 @@ end project(X, W; metric=:euclidean) = project!(copy(X), W; metric=:euclidean) # geodesic retraction, coincides with Stiefel retraction (which is not geodesic for p < n) -function retract(W::AbstractTensorMap, Δ::UnitaryTangent, α; alg=MAK.select_algorithm(left_polar!, W)) +function retract(W::AbstractTensorMap, Δ::UnitaryTangent, α; + alg=MAK.select_algorithm(left_polar!, W)) W == base(Δ) || throw(ArgumentError("not a valid tangent vector at base point")) E = exp(α * Δ.A) W′ = projectisometric!(W * E; alg) From 950fcad040fc1ce6331cb398eb6f1032598b86c0 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Fri, 3 Oct 2025 18:49:57 -0400 Subject: [PATCH 3/3] remove stray `tsvd` --- src/grassmann.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/grassmann.jl b/src/grassmann.jl index 8e9532e..d6ce744 100644 --- a/src/grassmann.jl +++ b/src/grassmann.jl @@ -31,8 +31,8 @@ end # output type of U, S, V in tsvd function _tsvd_types(Z::AbstractTensorMap) - TUSV = Core.Compiler.return_type(tsvd, Tuple{typeof(Z)}) - TU, TS, TV, = TUSV.types + TUSV = Core.Compiler.return_type(svd_compact, Tuple{typeof(Z)}) + TU, TS, TV = TUSV.types return TU, TS, TV end @@ -59,7 +59,7 @@ function Base.getproperty(Δ::GrassmannTangent, sym::Symbol) elseif sym ∈ (:U, :S, :V) v = Base.getfield(Δ, sym) v !== nothing && return v - U, S, V, = svd_compact(Δ.Z) + U, S, V = svd_compact(Δ.Z) Base.setfield!(Δ, :U, U) Base.setfield!(Δ, :S, S) Base.setfield!(Δ, :V, V)