From 94ace745ba037607a36f7b34a928672b6fde4ca0 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 13 Apr 2023 09:15:39 +0200 Subject: [PATCH 1/9] bump TensorKit compat --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index c5a887d..ec25f9c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,14 +1,14 @@ name = "TensorKitManifolds" uuid = "11fa318c-39cb-4a83-b1ed-cdc7ba1e3684" authors = ["Jutho Haegeman ", "Markus Hauru "] -version = "0.6.0" +version = "0.6.1" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" [compat] -TensorKit = "0.6,0.7,0.8,0.9,1" +TensorKit = "0.6,0.7,0.8,0.9,0.10,1" julia = "1" [extras] From cb9c5bb87ca2519bc562144f7e7520639be4cbd8 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 13 Apr 2023 10:01:32 +0200 Subject: [PATCH 2/9] update CI runs on - 1.6 (lts) - 1 (release) - nightly --- .github/workflows/ci.yml | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 6e2e669..16f5faf 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -11,9 +11,8 @@ jobs: fail-fast: false matrix: version: - # - '1.0' - - '1.4' - - '1.5' + - '1.6' + - '1' - 'nightly' os: - ubuntu-latest From adf1feaf3e696d5e4b81b6b272677dd0e528bdcd Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 13 Apr 2023 10:06:40 +0200 Subject: [PATCH 3/9] format style "yas" --- .JuliaFormatter.toml | 1 + .github/workflows/format_check.yml | 41 ++++++++ src/TensorKitManifolds.jl | 45 ++++----- src/auxiliary.jl | 154 ++++++++++++++--------------- src/grassmann.jl | 95 +++++++++--------- src/stiefel.jl | 134 +++++++++++++------------ src/unitary.jl | 82 ++++++++------- test/runtests.jl | 110 +++++++++++---------- 8 files changed, 366 insertions(+), 296 deletions(-) create mode 100644 .JuliaFormatter.toml create mode 100644 .github/workflows/format_check.yml diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 0000000..9613e05 --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1 @@ +style = "yas" \ No newline at end of file diff --git a/.github/workflows/format_check.yml b/.github/workflows/format_check.yml new file mode 100644 index 0000000..45f37d0 --- /dev/null +++ b/.github/workflows/format_check.yml @@ -0,0 +1,41 @@ +name: format-check + +on: + push: + branches: + - 'master' + tags: '*' + pull_request: + +jobs: + build: + runs-on: ${{ matrix.os }} + strategy: + matrix: + julia-version: [1] + julia-arch: [x86] + os: [ubuntu-latest] + steps: + - uses: julia-actions/setup-julia@latest + with: + version: ${{ matrix.julia-version }} + + - uses: actions/checkout@v1 + - name: Install JuliaFormatter and format + # This will use the latest version by default but you can set the version like so: + # + # julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter", version="0.13.0"))' + run: | + julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter"))' + julia -e 'using JuliaFormatter; format(".", verbose=true)' + - name: Format check + run: | + julia -e ' + out = Cmd(`git diff --name-only`) |> read |> String + if out == "" + exit(0) + else + @error "Some files have not been formatted !!!" + write(stdout, out) + exit(1) + end' \ No newline at end of file diff --git a/src/TensorKitManifolds.jl b/src/TensorKitManifolds.jl index 439907a..6ca33be 100644 --- a/src/TensorKitManifolds.jl +++ b/src/TensorKitManifolds.jl @@ -28,23 +28,24 @@ 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 eleps(X) = eps(real(eltype(X))) -function isisometry(W::AbstractTensorMap; tol = 10*eleps(W)) - WdW = W'*W +function isisometry(W::AbstractTensorMap; tol=10 * eleps(W)) + WdW = W' * W s = zero(float(real(eltype(W)))) - for (c,b) in blocks(WdW) + for (c, b) in blocks(WdW) _subtractone!(b) - s += dim(c)*length(b) + s += dim(c) * length(b) end - return norm(WdW) <= tol*sqrt(s) + return norm(WdW) <= tol * sqrt(s) end -isunitary(W::AbstractTensorMap; tol = 10*eleps(W)) = - isisometry(W; tol = tol) && isisometry(W'; tol = tol) +function isunitary(W::AbstractTensorMap; tol=10 * eleps(W)) + return isisometry(W; tol=tol) && isisometry(W'; tol=tol) +end function projecthermitian!(W::AbstractTensorMap) codomain(W) == domain(W) || throw(DomainError("Tensor with distinct domain and codomain cannot be hermitian.")) - for (c,b) in blocks(W) + for (c, b) in blocks(W) _projecthermitian!(b) end return W @@ -52,7 +53,7 @@ end function projectantihermitian!(W::AbstractTensorMap) codomain(W) == domain(W) || throw(DomainError("Tensor with distinct domain and codomain cannot be anithermitian.")) - for (c,b) in blocks(W) + for (c, b) in blocks(W) _projectantihermitian!(b) end return W @@ -60,18 +61,18 @@ end struct PolarNewton <: TensorKit.OrthogonalFactorizationAlgorithm end -function projectisometric!(W::AbstractTensorMap; alg = Polar()) +function projectisometric!(W::AbstractTensorMap; alg=Polar()) if alg isa TensorKit.Polar || alg isa TensorKit.SDD - foreach(blocks(W)) do (c,b) - _polarsdd!(b) + foreach(blocks(W)) do (c, b) + return _polarsdd!(b) end elseif alg isa TensorKit.SVD - foreach(blocks(W)) do (c,b) - _polarsvd!(b) + foreach(blocks(W)) do (c, b) + return _polarsvd!(b) end elseif alg isa PolarNewton - foreach(blocks(W)) do (c,b) - _polarnewton!(b) + foreach(blocks(W)) do (c, b) + return _polarnewton!(b) end else throw(ArgumentError("unkown algorithm for projectisometric!: alg = $alg")) @@ -80,14 +81,14 @@ function projectisometric!(W::AbstractTensorMap; alg = Polar()) end function projectcomplement!(X::AbstractTensorMap, W::AbstractTensorMap; - tol = 10*eleps(X)) - P = W'*X + tol=10 * eleps(X)) + P = W' * X nP = norm(P) nX = norm(X) dP = dim(P) - while nP > tol*max(dP, nX) + while nP > tol * max(dP, nX) X = mul!(X, W, P, -1, 1) - P = W'*X + P = W' * X nP = norm(P) end return X @@ -97,11 +98,11 @@ projecthermitian(W::AbstractTensorMap) = projecthermitian!(copy(W)) projectantihermitian(W::AbstractTensorMap) = projectantihermitian!(copy(W)) function projectisometric(W::AbstractTensorMap; - alg::TensorKit.OrthogonalFactorizationAlgorithm = Polar()) + alg::TensorKit.OrthogonalFactorizationAlgorithm=Polar()) return projectisometric!(copy(W); alg=alg) end function projectcomplement(X::AbstractTensorMap, W::AbstractTensorMap, - tol = 10*eleps(X)) + tol=10 * eleps(X)) return projectcomplement!(copy(X), W; tol=tol) end diff --git a/src/auxiliary.jl b/src/auxiliary.jl index 01ac5c7..5ba8787 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -1,24 +1,24 @@ using LinearAlgebra function _projecthermitian!(a::StridedMatrix) - L = size(a,1) - for J = 1:16:L - for I = 1:16:L + L = size(a, 1) + for J in 1:16:L + for I in 1:16:L if I == J - for j = J:min(J+15,L) - for i = J:j-1 - x = (a[i,j] + conj(a[j,i]))/2 - a[i,j] = x - a[j,i] = conj(x) + for j in J:min(J + 15, L) + for i in J:(j - 1) + x = (a[i, j] + conj(a[j, i])) / 2 + a[i, j] = x + a[j, i] = conj(x) end - a[j,j] = (a[j,j] + conj(a[j,j]))/2 + a[j, j] = (a[j, j] + conj(a[j, j])) / 2 end else - for j = J:min(J+15,L) - for i = I:min(I+15,L) - x = (a[i,j] + conj(a[j,i]))/2 - a[i,j] = x - a[j,i] = conj(x) + for j in J:min(J + 15, L) + for i in I:min(I + 15, L) + x = (a[i, j] + conj(a[j, i])) / 2 + a[i, j] = x + a[j, i] = conj(x) end end end @@ -27,29 +27,29 @@ function _projecthermitian!(a::StridedMatrix) return a end function _projecthermitian!(a::AbstractMatrix) - a .= (a .+ a')./2 + a .= (a .+ a') ./ 2 return a end function _projectantihermitian!(a::StridedMatrix) - L = size(a,1) - for J = 1:16:L - for I = 1:16:L + L = size(a, 1) + for J in 1:16:L + for I in 1:16:L if I == J - for j = J:min(J+15,L) - for i = J:j-1 - x = (a[i,j] - conj(a[j,i]))/2 - a[i,j] = x - a[j,i] = -conj(x) + for j in J:min(J + 15, L) + for i in J:(j - 1) + x = (a[i, j] - conj(a[j, i])) / 2 + a[i, j] = x + a[j, i] = -conj(x) end - a[j,j] = (a[j,j] - conj(a[j,j]))/2 + a[j, j] = (a[j, j] - conj(a[j, j])) / 2 end else - for j = J:min(J+15,L) - for i = I:min(I+15,L) - x = (a[i,j] - conj(a[j,i]))/2 - a[i,j] = x - a[j,i] = -conj(x) + for j in J:min(J + 15, L) + for i in I:min(I + 15, L) + x = (a[i, j] - conj(a[j, i])) / 2 + a[i, j] = x + a[j, i] = -conj(x) end end end @@ -58,7 +58,7 @@ function _projectantihermitian!(a::StridedMatrix) return a end function _projectantihermitian!(a::AbstractMatrix) - a .= (a .- a')./2 + a .= (a .- a') ./ 2 return a end @@ -72,23 +72,23 @@ function _subtractone!(a::AbstractMatrix) end function _one!(A::AbstractMatrix) m, n = size(A) - @inbounds for j = 1:n - for i = 1:m - A[i,j] = i==j + @inbounds for j in 1:n + for i in 1:m + A[i, j] = i == j end end return A end function _polarsdd!(A::StridedMatrix) - U, S, V = svd!(A; alg = LinearAlgebra.DivideAndConquer()) + 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()) + U, S, V = svd!(A; alg=LinearAlgebra.QRIteration()) return mul!(A, U, V') end -function _polarnewton!(A::StridedMatrix; tol = 10*eleps(A), maxiter = 5) +function _polarnewton!(A::StridedMatrix; tol=10 * eleps(A), maxiter=5) m, n = size(A) @assert m >= n A2 = copy(A) @@ -97,7 +97,7 @@ function _polarnewton!(A::StridedMatrix; tol = 10*eleps(A), maxiter = 5) R, Ri = _avgdiff!(R, Ri) i = 1 R2 = view(A, 1:n, 1:n) - fill!(view(A, n+1:m, 1:n), zero(eltype(A))) + fill!(view(A, (n + 1):m, 1:n), zero(eltype(A))) copyto!(R2, R) while maximum(abs, Ri) > tol if i == maxiter # if not converged by now, fall back to sdd @@ -118,8 +118,8 @@ function _avgdiff!(A::AbstractArray, B::AbstractArray) @inbounds begin a = A[I] b = B[I] - A[I] = (a+b)/2 - B[I] = b-a + A[I] = (a + b) / 2 + B[I] = b - a end end return A, B @@ -133,84 +133,84 @@ end # Here, Q′ is a set of orthogonal columns to the colums in W′. function _stiefelexp(W::StridedMatrix, A::StridedMatrix, Z::StridedMatrix, α) n, p = size(W) - r = min(2*p, n) + r = min(2 * p, n) QQ, _ = qr!([W Z]) - Q = similar(W, n, r-p) - @inbounds for j = Base.OneTo(r-p) - for i = Base.OneTo(n) - Q[i,j] = (i == p+j) + Q = similar(W, n, r - p) + @inbounds for j in Base.OneTo(r - p) + for i in Base.OneTo(n) + Q[i, j] = (i == p + j) end end Q = lmul!(QQ, Q) - R = Q'*Z - A2 = similar(A, min(2*p, n), min(2*p, n)) + R = Q' * Z + A2 = similar(A, min(2 * p, n), min(2 * p, n)) A2[1:p, 1:p] .= α .* A - A2[p+1:end, 1:p] .= α .* R - A2[1:p, p+1:end] .= (-α) .* (R') - A2[p+1:end, p+1:end] .= 0 - U = [W Q]*exp(A2) + A2[(p + 1):end, 1:p] .= α .* R + A2[1:p, (p + 1):end] .= (-α) .* (R') + A2[(p + 1):end, (p + 1):end] .= 0 + U = [W Q] * exp(A2) U = _polarnewton!(U) - W′ = U[:,1:p] - Q′ = U[:,p+1:end] + W′ = U[:, 1:p] + Q′ = U[:, (p + 1):end] R′ = R return W′, Q, Q′, R′ end function _stiefellog(Wold::StridedMatrix, Wnew::StridedMatrix; - tol = 10*eleps(Wold), maxiter = 100) + tol=10 * eleps(Wold), maxiter=100) n, p = size(Wold) - r = min(2*p, n) - P = Wold'*Wnew - dW = Wnew - Wold*P + r = min(2 * p, n) + P = Wold' * Wnew + dW = Wnew - Wold * P QQ, _ = qr!([Wold dW]) - Q = similar(Wold, n, r-p) - @inbounds for j = Base.OneTo(r-p) - for i = Base.OneTo(n) - Q[i,j] = (i == p+j) + Q = similar(Wold, n, r - p) + @inbounds for j in Base.OneTo(r - p) + for i in Base.OneTo(n) + Q[i, j] = (i == p + j) end end Q = lmul!(QQ, Q) - R = Q'*dW - Wext = [Wold Q]; + R = Q' * dW + Wext = [Wold Q] F = qr!([P; R]) U = lmul!(F.Q, _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) + 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)') + mul!(X, X * (YSVD.V), (YSVD.U)') UsqrtS = YSVD.U - @inbounds for j = 1:size(UsqrtS, 2) + @inbounds for j in 1:size(UsqrtS, 2) s = sqrt(YSVD.S[j]) - @simd for i = 1:size(UsqrtS, 1) - UsqrtS[i,j] *= s + @simd for i in 1:size(UsqrtS, 1) + UsqrtS[i, j] *= s end end mul!(Y, UsqrtS, UsqrtS') end logU = _projectantihermitian!(log(U)) if eltype(U) <: Real - @assert mapreduce(abs ∘ imag, max, logU; init = abs(zero(eltype(logU))) ) <= tol + @assert mapreduce(abs ∘ imag, max, logU; init=abs(zero(eltype(logU)))) <= tol K = real(logU) else K = logU end - C = view(K, p+1:r, p+1:r) + C = view(K, (p + 1):r, (p + 1):r) i = 1 - τ = mapreduce(abs, max, C; init = abs(zero(eltype(C)))) + τ = mapreduce(abs, max, C; init=abs(zero(eltype(C)))) while τ > tol if i > maxiter @warn "Stiefel logarithm: not converged in $maxiter iterations, τ = $τ" break end - eC = exp(rmul!(C,-1)) - X .= X*eC - Y .= Y*eC + eC = exp(rmul!(C, -1)) + X .= X * eC + Y .= Y * eC logU = _projectantihermitian!(log(U)) if eltype(U) <: Real - @assert mapreduce(abs ∘ imag, max, logU; init = abs(zero(eltype(logU))) ) <= tol + @assert mapreduce(abs ∘ imag, max, logU; init=abs(zero(eltype(logU)))) <= tol K .= real.(logU) else K .= logU @@ -218,5 +218,5 @@ function _stiefellog(Wold::StridedMatrix, Wnew::StridedMatrix; τ = maximum(abs, C) i += 1 end - return K[1:p,1:p], Q, K[p+1:r, 1:p] + return K[1:p, 1:p], Q, K[(p + 1):r, 1:p] end diff --git a/src/grassmann.jl b/src/grassmann.jl index c7a81b1..0243fac 100644 --- a/src/grassmann.jl +++ b/src/grassmann.jl @@ -22,7 +22,7 @@ mutable struct GrassmannTangent{T<:AbstractTensorMap, S::Union{Nothing,TS} V::Union{Nothing,TV} function GrassmannTangent(W::AbstractTensorMap{S,N₁,N₂}, - Z::AbstractTensorMap{S,N₁,N₂}) where {S,N₁,N₂} + Z::AbstractTensorMap{S,N₁,N₂}) where {S,N₁,N₂} T = typeof(W) TT = promote_type(float(eltype(W)), eltype(Z)) G = sectortype(W) @@ -46,8 +46,10 @@ end Base.getindex(Δ::GrassmannTangent) = Δ.Z base(Δ::GrassmannTangent) = Δ.W -checkbase(Δ₁::GrassmannTangent, Δ₂::GrassmannTangent) = Δ₁.W == Δ₂.W ? Δ₁.W : - throw(ArgumentError("tangent vectors with different base points")) +function checkbase(Δ₁::GrassmannTangent, Δ₂::GrassmannTangent) + return Δ₁.W == Δ₂.W ? Δ₁.W : + throw(ArgumentError("tangent vectors with different base points")) +end function Base.getproperty(Δ::GrassmannTangent, sym::Symbol) if sym ∈ (:W, :Z) @@ -67,15 +69,17 @@ function Base.getproperty(Δ::GrassmannTangent, sym::Symbol) end end function Base.setproperty!(Δ::GrassmannTangent, sym::Symbol, v) - error("type GrassmannTangent does not allow to change its fields") + return error("type GrassmannTangent does not allow to change its fields") end # Basic vector space behaviour -Base.:+(Δ₁::GrassmannTangent, Δ₂::GrassmannTangent) = - GrassmannTangent(checkbase(Δ₁,Δ₂), Δ₁.Z + Δ₂.Z) -Base.:-(Δ₁::GrassmannTangent, Δ₂::GrassmannTangent) = - GrassmannTangent(checkbase(Δ₁,Δ₂), Δ₁.Z - Δ₂.Z) -Base.:-(Δ::GrassmannTangent) = (-1)*Δ +function Base.:+(Δ₁::GrassmannTangent, Δ₂::GrassmannTangent) + return GrassmannTangent(checkbase(Δ₁, Δ₂), Δ₁.Z + Δ₂.Z) +end +function Base.:-(Δ₁::GrassmannTangent, Δ₂::GrassmannTangent) + return GrassmannTangent(checkbase(Δ₁, Δ₂), Δ₁.Z - Δ₂.Z) +end +Base.:-(Δ::GrassmannTangent) = (-1) * Δ Base.:*(Δ::GrassmannTangent, α::Number) = rmul!(copy(Δ), α) Base.:*(α::Number, Δ::GrassmannTangent) = lmul!(α, copy(Δ)) @@ -125,12 +129,12 @@ function TensorKit.dot(Δ₁::GrassmannTangent, Δ₂::GrassmannTangent) checkbase(Δ₁, Δ₂) return dot(Δ₁.Z, Δ₂.Z) end -TensorKit.norm(Δ::GrassmannTangent, p::Real = 2) = norm(Δ.Z, p) +TensorKit.norm(Δ::GrassmannTangent, p::Real=2) = norm(Δ.Z, p) # tangent space methods -function project!(X::AbstractTensorMap, W::AbstractTensorMap; metric = :euclidean) +function project!(X::AbstractTensorMap, W::AbstractTensorMap; metric=:euclidean) @assert metric == :euclidean - P = W'*X + P = W' * X Z = mul!(X, W, P, -1, 1) Z = projectcomplement!(Z, W) return GrassmannTangent(W, Z) @@ -143,13 +147,14 @@ Project X onto the Grassmann tangent space at the base point `W`, which is assum isometric, i.e. `W'*W ≈ id(domain(W))`. The resulting tensor `Z` in the tangent space of `W` is given by `Z = X - W * (W'*X)` and satisfies `W'*Z = 0`. """ -project(X::AbstractTensorMap, W::AbstractTensorMap; metric = :euclidean) = - project!(copy(X), W; metric=metric) +function project(X::AbstractTensorMap, W::AbstractTensorMap; metric=:euclidean) + return project!(copy(X), W; metric=metric) +end function inner(W::AbstractTensorMap, Δ₁::GrassmannTangent, Δ₂::GrassmannTangent; - metric = :euclidean) + metric=:euclidean) @assert metric == :euclidean - Δ₁ === Δ₂ ? norm(Δ₁)^2 : real(dot(Δ₁, Δ₂)) + return Δ₁ === Δ₂ ? norm(Δ₁)^2 : real(dot(Δ₁, Δ₂)) end """ @@ -165,15 +170,15 @@ while the local tangent vector along the retraction curve is `Z′ = - W * V' * sin(α*S) * S * V + U * cos(α * S) * S * V'`. """ -function retract(W::AbstractTensorMap, Δ::GrassmannTangent, α; alg = nothing) +function retract(W::AbstractTensorMap, Δ::GrassmannTangent, α; alg=nothing) W == base(Δ) || throw(ArgumentError("not a valid tangent vector at base point")) U, S, V = Δ.U, Δ.S, Δ.V - WVd = W*V' + WVd = W * V' sSV, cSV = _sincosSV(α, S, V) # sin(S)*V, cos(S)*V - W′ = projectisometric!(WVd*cSV + U*sSV) + W′ = projectisometric!(WVd * cSV + U * sSV) sSSV = _lmul!(S, sSV) # sin(S)*S*V cSSV = _lmul!(S, cSV) # cos(S)*S*V - Z′ = projectcomplement!(-WVd*sSSV + U*cSSV, W′) + Z′ = projectcomplement!(-WVd * sSSV + U * cSSV, W′) return W′, GrassmannTangent(W′, Z′) end @@ -186,7 +191,7 @@ This is done by solving the equation `Wold * V' * cos(S) * V + U * sin(S) * V = for the isometries `U`, `V`, and `Y`, and the diagonal matrix `S`, and returning `Z = U * S * V` and `Y`. """ -function invretract(Wold::AbstractTensorMap, Wnew::AbstractTensorMap; alg = nothing) +function invretract(Wold::AbstractTensorMap, Wnew::AbstractTensorMap; alg=nothing) space(Wold) == space(Wnew) || throw(SectorMismatch()) WodWn = Wold' * Wnew # V' * cos(S) * V * Y Wneworth = Wnew - Wold * WodWn @@ -195,8 +200,8 @@ function invretract(Wold::AbstractTensorMap, Wnew::AbstractTensorMap; alg = noth # acos always returns a complex TensorMap. We cast back to real if possible. S = eltype(WodWn) <: Real && isreal(sectortype(Scmplx)) ? real(Scmplx) : Scmplx UsS = Wneworth * VY' # U * sin(S) # should be in polar decomposition form - U = projectisometric!(UsS; alg = Polar()) - Y = Vd*VY + U = projectisometric!(UsS; alg=Polar()) + Y = Vd * VY V = Vd' Z = Grassmann.GrassmannTangent(Wold, U * S * V) return Z, Y @@ -211,24 +216,24 @@ retraction. """ function relativegauge(W::AbstractTensorMap, V::AbstractTensorMap) space(W) == space(V) || throw(SectorMismatch()) - return projectisometric!(V'*W; alg = Polar()) + return projectisometric!(V' * W; alg=Polar()) end function transport!(Θ::GrassmannTangent, W::AbstractTensorMap, Δ::GrassmannTangent, α, W′; - alg = nothing) - W == checkbase(Δ,Θ) || throw(ArgumentError("not a valid tangent vector at base point")) + alg=nothing) + W == checkbase(Δ, Θ) || throw(ArgumentError("not a valid tangent vector at base point")) U, S, V = Δ.U, Δ.S, Δ.V - WVd = W*V' - UdΘ = U'*Θ.Z + WVd = W * V' + UdΘ = U' * Θ.Z sSUdθ, cSUdθ = _sincosSV(α, S, UdΘ) # sin(S)*U'*Θ, cos(S)*U'*Θ cSm1UdΘ = axpy!(-1, UdΘ, cSUdθ) # (cos(S)-1)*U'*Θ - Z′ = axpy!(true, U*cSm1UdΘ - WVd*sSUdθ, Θ.Z) + Z′ = axpy!(true, U * cSm1UdΘ - WVd * sSUdθ, Θ.Z) Z′ = projectcomplement!(Z′, W′) return GrassmannTangent(W′, Z′) end function transport(Θ::GrassmannTangent, W::AbstractTensorMap, Δ::GrassmannTangent, α, W′; - alg = nothing) - return transport!(copy(Θ), W, Δ, α, W′; alg = alg) + alg=nothing) + return transport!(copy(Θ), W, Δ, α, W′; alg=alg) end # auxiliary methods: unsafe, no checking @@ -237,16 +242,16 @@ function _sincosSV(α::Real, S::AbstractTensorMap, V::AbstractTensorMap) # S is assumed diagonal cSV = similar(V) sSV = similar(V) - @inbounds for (c,bS) in blocks(S) - bcSV = block(cSV,c) - bsSV = block(sSV,c) - bV = block(V,c) - Threads.@threads for j = 1:size(bV,2) - @simd for i = 1:size(bV, 1) - sS, cS = sincos(α*bS[i,i]) + @inbounds for (c, bS) in blocks(S) + bcSV = block(cSV, c) + bsSV = block(sSV, c) + bV = block(V, c) + Threads.@threads for j in 1:size(bV, 2) + @simd for i in 1:size(bV, 1) + sS, cS = sincos(α * bS[i, i]) # TODO: we are computing sin and cos above within the loop over j, while it is independent; moving it out the loop requires extra storage though. - bsSV[i,j] = sS*bV[i,j] - bcSV[i,j] = cS*bV[i,j] + bsSV[i, j] = sS * bV[i, j] + bcSV[i, j] = cS * bV[i, j] end end end @@ -254,11 +259,11 @@ function _sincosSV(α::Real, S::AbstractTensorMap, V::AbstractTensorMap) end # multiply V with diagonal S in place, S is assumed diagonal function _lmul!(S::AbstractTensorMap, V::AbstractTensorMap) - @inbounds for (c,bS) in blocks(S) - bV = block(V,c) - Threads.@threads for j = 1:size(bV,2) - @simd for i = 1:size(bV, 1) - bV[i,j] *= bS[i,i] + @inbounds for (c, bS) in blocks(S) + bV = block(V, c) + Threads.@threads for j in 1:size(bV, 2) + @simd for i in 1:size(bV, 1) + bV[i, j] *= bS[i, i] end end end diff --git a/src/stiefel.jl b/src/stiefel.jl index f48b338..bd28f0e 100644 --- a/src/stiefel.jl +++ b/src/stiefel.jl @@ -10,10 +10,10 @@ using ..TensorKitManifolds: projecthermitian!, projectantihermitian!, projectisometric!, projectcomplement!, PolarNewton, _stiefelexp, _stiefellog, eleps import ..TensorKitManifolds: base, checkbase, - inner, retract, transport, transport! + inner, retract, transport, transport! # special type to store tangent vectors using A and Z = Q*R, -struct StiefelTangent{T<:AbstractTensorMap, TA<:AbstractTensorMap} +struct StiefelTangent{T<:AbstractTensorMap,TA<:AbstractTensorMap} W::T A::TA Z::T @@ -25,15 +25,19 @@ end Base.getindex(Δ::StiefelTangent) = Δ.W * Δ.A + Δ.Z base(Δ::StiefelTangent) = Δ.W -checkbase(Δ₁::StiefelTangent, Δ₂::StiefelTangent) = Δ₁.W == Δ₂.W ? Δ₁.W : - throw(ArgumentError("tangent vectors with different base points")) +function checkbase(Δ₁::StiefelTangent, Δ₂::StiefelTangent) + return Δ₁.W == Δ₂.W ? Δ₁.W : + throw(ArgumentError("tangent vectors with different base points")) +end # Basic vector space behaviour -Base.:+(Δ₁::StiefelTangent, Δ₂::StiefelTangent) = - StiefelTangent(checkbase(Δ₁,Δ₂), Δ₁.A + Δ₂.A, Δ₁.Z + Δ₂.Z) -Base.:-(Δ₁::StiefelTangent, Δ₂::StiefelTangent) = - StiefelTangent(checkbase(Δ₁,Δ₂), Δ₁.A - Δ₂.A, Δ₁.Z - Δ₂.Z) -Base.:-(Δ::StiefelTangent) = (-1)*Δ +function Base.:+(Δ₁::StiefelTangent, Δ₂::StiefelTangent) + return StiefelTangent(checkbase(Δ₁, Δ₂), Δ₁.A + Δ₂.A, Δ₁.Z + Δ₂.Z) +end +function Base.:-(Δ₁::StiefelTangent, Δ₂::StiefelTangent) + return StiefelTangent(checkbase(Δ₁, Δ₂), Δ₁.A - Δ₂.A, Δ₁.Z - Δ₂.Z) +end +Base.:-(Δ::StiefelTangent) = (-1) * Δ Base.:*(Δ::StiefelTangent, α::Real) = rmul!(copy(Δ), α) Base.:*(α::Real, Δ::StiefelTangent) = lmul!(α, copy(Δ)) @@ -67,14 +71,13 @@ end function TensorKit.dot(Δ₁::StiefelTangent, Δ₂::StiefelTangent) checkbase(Δ₁, Δ₂) - dot(Δ₁.A, Δ₂.A) + dot(Δ₁.Z, Δ₂.Z) + return dot(Δ₁.A, Δ₂.A) + dot(Δ₁.Z, Δ₂.Z) end -TensorKit.norm(Δ::StiefelTangent, p::Real = 2) = - norm((norm(Δ.A, p), norm(Δ.Z, p)), p) +TensorKit.norm(Δ::StiefelTangent, p::Real=2) = norm((norm(Δ.A, p), norm(Δ.Z, p)), p) # tangent space methods function inner(W::AbstractTensorMap, Δ₁::StiefelTangent, Δ₂::StiefelTangent; - metric = :euclidean) + metric=:euclidean) if metric == :euclidean return inner_euclidean(W, Δ₁, Δ₂) elseif metric == :canonical @@ -83,7 +86,7 @@ function inner(W::AbstractTensorMap, Δ₁::StiefelTangent, Δ₂::StiefelTangen throw(ArgumentError("unknown metric: $metric")) end end -function project!(X::AbstractTensorMap, W::AbstractTensorMap; metric = :euclidean) +function project!(X::AbstractTensorMap, W::AbstractTensorMap; metric=:euclidean) if metric == :euclidean return project_euclidean!(X, W) elseif metric == :canonical @@ -92,10 +95,11 @@ function project!(X::AbstractTensorMap, W::AbstractTensorMap; metric = :euclidea throw(ArgumentError("unknown metric: `metric = $metric`")) end end -project(X::AbstractTensorMap, W::AbstractTensorMap; metric = :euclidean) = - project!(copy(X), W; metric = metric) +function project(X::AbstractTensorMap, W::AbstractTensorMap; metric=:euclidean) + return project!(copy(X), W; metric=metric) +end -function retract(W::AbstractTensorMap, Δ::StiefelTangent, α::Real; alg = :exp) +function retract(W::AbstractTensorMap, Δ::StiefelTangent, α::Real; alg=:exp) if alg == :exp return retract_exp(W, Δ, α) elseif alg == :cayley @@ -104,7 +108,7 @@ function retract(W::AbstractTensorMap, Δ::StiefelTangent, α::Real; alg = :exp) throw(ArgumentError("unknown algorithm: `alg = $alg`")) end end -function invretract(Wold::AbstractTensorMap, Wnew::AbstractTensorMap; alg = :exp) +function invretract(Wold::AbstractTensorMap, Wnew::AbstractTensorMap; alg=:exp) if alg == :exp return invretract_exp(Wold, Wnew) elseif alg == :cayley @@ -114,7 +118,7 @@ function invretract(Wold::AbstractTensorMap, Wnew::AbstractTensorMap; alg = :exp end end function transport!(Θ::StiefelTangent, W::AbstractTensorMap, Δ::StiefelTangent, α::Real, W′; - alg = :exp) + alg=:exp) if alg == :exp return transport_exp!(Θ, W, Δ, α, W′) elseif alg == :cayley @@ -124,16 +128,16 @@ function transport!(Θ::StiefelTangent, W::AbstractTensorMap, Δ::StiefelTangent end end function transport(Θ::StiefelTangent, W::AbstractTensorMap, Δ::StiefelTangent, α::Real, W′; - alg = :exp) - return transport!(copy(Θ), W, Δ, α, W′; alg = alg) + alg=:exp) + return transport!(copy(Θ), W, Δ, α, W′; alg=alg) end # euclidean metric function inner_euclidean(W::AbstractTensorMap, Δ₁::StiefelTangent, Δ₂::StiefelTangent) - Δ₁ === Δ₂ ? norm(Δ₁)^2 : real(dot(Δ₁,Δ₂)) + return Δ₁ === Δ₂ ? norm(Δ₁)^2 : real(dot(Δ₁, Δ₂)) end function project_euclidean!(X::AbstractTensorMap, W::AbstractTensorMap) - P = W'*X + P = W' * X Z = mul!(X, W, P, -1, 1) A = projectantihermitian!(P) Z = projectcomplement!(Z, W) @@ -144,13 +148,13 @@ project_euclidean(X, W) = project_euclidean!(copy(X), W) # canonical metric function inner_canonical(W::AbstractTensorMap, Δ₁::StiefelTangent, Δ₂::StiefelTangent) if Δ₁ === Δ₂ - return (norm(Δ₁.A)^2)/2 + norm(Δ₁.Z)^2 + return (norm(Δ₁.A)^2) / 2 + norm(Δ₁.Z)^2 else - return real(dot(Δ₁.A, Δ₂.A)/2 + dot(Δ₁.Z, Δ₂.Z)) + return real(dot(Δ₁.A, Δ₂.A) / 2 + dot(Δ₁.Z, Δ₂.Z)) end end function project_canonical!(X::AbstractTensorMap, W::AbstractTensorMap) - P = W'*X + P = W' * X Z = mul!(X, W, P, -1, 1) A = rmul!(projectantihermitian!(P), 2) Z = projectcomplement!(Z, W) @@ -165,11 +169,11 @@ function stiefelexp(W::AbstractTensorMap, α::Real) S = spacetype(W) G = sectortype(W) - Wdata′ = TensorKit.SectorDict{G, storagetype(W)}() - Qdata = TensorKit.SectorDict{G, storagetype(W)}() - Qdata′ = TensorKit.SectorDict{G, storagetype(W)}() - Rdata′ = TensorKit.SectorDict{G, storagetype(W)}() - dims = TensorKit.SectorDict{G, Int}() + Wdata′ = TensorKit.SectorDict{G,storagetype(W)}() + Qdata = TensorKit.SectorDict{G,storagetype(W)}() + Qdata′ = TensorKit.SectorDict{G,storagetype(W)}() + Rdata′ = TensorKit.SectorDict{G,storagetype(W)}() + dims = TensorKit.SectorDict{G,Int}() for c in blocksectors(W) w′, q, q′, r′ = _stiefelexp(block(W, c), block(A, c), block(Z, c), α) Wdata′[c] = w′ @@ -180,9 +184,9 @@ function stiefelexp(W::AbstractTensorMap, end V = S(dims) W′ = TensorMap(Wdata′, space(W)) - Q = TensorMap(Qdata, codomain(W)←V) - Q′ = TensorMap(Qdata′, codomain(W)←V) - R′ = TensorMap(Rdata′, V←domain(W)) + Q = TensorMap(Qdata, codomain(W) ← V) + Q′ = TensorMap(Qdata′, codomain(W) ← V) + R′ = TensorMap(Rdata′, V ← domain(W)) return W′, Q, Q′, R′ end @@ -191,23 +195,23 @@ function retract_exp(W::AbstractTensorMap, Δ::StiefelTangent, α::Real) W == base(Δ) || throw(ArgumentError("not a valid tangent vector at base point")) W′, Q, Q′, R′ = stiefelexp(W, Δ.A, Δ.Z, α) A′ = Δ.A - Z′ = projectcomplement!(Q′*R′, W′) # to ensure orthogonality + Z′ = projectcomplement!(Q′ * R′, W′) # to ensure orthogonality return W′, StiefelTangent(W′, A′, Z′) end function invretract_exp(Wold::AbstractTensorMap, Wnew::AbstractTensorMap; - tol = eleps(Wold)^(2/3)) + tol=eleps(Wold)^(2 / 3)) space(Wold) == space(Wnew) || throw(SectorMismatch()) S = spacetype(Wold) G = sectortype(Wold) - Adata = TensorKit.SectorDict{G, storagetype(Wold)}() - Zdata = TensorKit.SectorDict{G, storagetype(Wold)}() + Adata = TensorKit.SectorDict{G,storagetype(Wold)}() + Zdata = TensorKit.SectorDict{G,storagetype(Wold)}() for c in blocksectors(Wold) - a, q, r = _stiefellog(block(Wold, c), block(Wnew, c); tol = tol) + a, q, r = _stiefellog(block(Wold, c), block(Wnew, c); tol=tol) Adata[c] = a - Zdata[c] = q*r + Zdata[c] = q * r end - A = TensorMap(Adata, domain(Wold)←domain(Wold)) + A = TensorMap(Adata, domain(Wold) ← domain(Wold)) Z = TensorMap(Zdata, space(Wold)) return StiefelTangent(Wold, A, Z) end @@ -218,43 +222,45 @@ end # can be computed efficiently: O(np^2) + O(p^3) function transport_exp!(Θ::StiefelTangent, W::AbstractTensorMap, Δ::StiefelTangent, α::Real, W′::AbstractTensorMap) - W == checkbase(Δ,Θ) || throw(ArgumentError("not a valid tangent vector at base point")) + W == checkbase(Δ, Θ) || throw(ArgumentError("not a valid tangent vector at base point")) # TODO: stiefelexp call does not depend on Θ # cache result or find some other way not to recompute this information _W′, Q, Q′, R′ = stiefelexp(W, Δ.A, Δ.Z, α) W′ ≈ _W′ || throw(ArgumentError("not a valid tangent vector at end point")) A = Θ.A Z = Θ.Z - QZ = Q'*Θ.Z + QZ = Q' * Θ.Z A′ = Θ.A - Z′ = projectcomplement!(mul!(Z, (Q′-Q), QZ, 1, 1), W′) + Z′ = projectcomplement!(mul!(Z, (Q′ - Q), QZ, 1, 1), W′) return StiefelTangent(W′, A′, Z′) end -transport_exp(Θ::StiefelTangent, W::AbstractTensorMap, Δ::StiefelTangent, α::Real, W′) = - transport_exp!(copy(Θ), W, Δ, α, W′) +function transport_exp(Θ::StiefelTangent, W::AbstractTensorMap, Δ::StiefelTangent, α::Real, + W′) + return transport_exp!(copy(Θ), W, Δ, α, W′) +end # Cayley retraction, slightly more efficient than above? # can be computed efficiently: O(np^2) + O(p^3) function retract_cayley(W::AbstractTensorMap, Δ::StiefelTangent, α::Real) W == base(Δ) || throw(ArgumentError("not a valid tangent vector at base point")) A, Z = Δ.A, Δ.Z - ZdZ = Z'*Z - X = axpy!(α^2/4, ZdZ, axpy!(-α/2, A, one(A))) + ZdZ = Z' * Z + X = axpy!(α^2 / 4, ZdZ, axpy!(-α / 2, A, one(A))) iX = inv(X) - W′ = projectisometric!((2*W+α*Z)*iX - W; alg = PolarNewton()) - A′ = projectantihermitian!((A - (α/2)*ZdZ)*iX) - Z′ = (Z-α*(W+α/2*Z)*(iX*ZdZ)) - Z′ = projectcomplement!(Z′*projecthermitian!(iX), W′) + W′ = projectisometric!((2 * W + α * Z) * iX - W; alg=PolarNewton()) + A′ = projectantihermitian!((A - (α / 2) * ZdZ) * iX) + Z′ = (Z - α * (W + α / 2 * Z) * (iX * ZdZ)) + Z′ = projectcomplement!(Z′ * projecthermitian!(iX), W′) return W′, StiefelTangent(W′, A′, Z′) end function invretract_cayley(Wold::AbstractTensorMap, Wnew::AbstractTensorMap) space(Wnew) == space(Wold) || throw(SpaceMismatch()) - P = Wold'*Wnew - iX = rmul!(axpy!(1,P,one(P)), 1/2) + P = Wold' * Wnew + iX = rmul!(axpy!(1, P, one(P)), 1 / 2) X = inv(iX) - Z = projectcomplement!(Wnew - Wold*P, Wold)*X - A = projectantihermitian!(rmul!(axpy!(-1, X, mul!(one(X), Z', Z, 1/4, 1)), 2)) + Z = projectcomplement!(Wnew - Wold * P, Wold) * X + A = projectantihermitian!(rmul!(axpy!(-1, X, mul!(one(X), Z', Z, 1 / 4, 1)), 2)) return StiefelTangent(Wold, A, Z) end @@ -262,16 +268,18 @@ end # isometric for both euclidean and canonical metric # can be computed efficiently: O(np^2) + O(p^3) function transport_cayley!(Θ::StiefelTangent, W::AbstractTensorMap, Δ::StiefelTangent, - α::Real, W′) - W == checkbase(Δ,Θ) || throw(ArgumentError("not a valid tangent vector at base point")) + α::Real, W′) + W == checkbase(Δ, Θ) || throw(ArgumentError("not a valid tangent vector at base point")) A, Z = Δ.A, Δ.Z - X = axpy!(α^2/4, Z'*Z, axpy!(-α/2, A, one(A))) + X = axpy!(α^2 / 4, Z' * Z, axpy!(-α / 2, A, one(A))) A′ = Θ.A - ZdZ = Z'*Θ.Z - Z′ = projectcomplement!(axpy!(-α, (W+(α/2)*Z)*(X\ZdZ), Θ.Z), W′) + ZdZ = Z' * Θ.Z + Z′ = projectcomplement!(axpy!(-α, (W + (α / 2) * Z) * (X \ ZdZ), Θ.Z), W′) return StiefelTangent(W′, A′, Z′) end -transport_cayley(Θ::StiefelTangent, W::AbstractTensorMap, Δ::StiefelTangent, α::Real, W′) = - transport_cayley!(copy(Θ), W, Δ, α, W′) +function transport_cayley(Θ::StiefelTangent, W::AbstractTensorMap, Δ::StiefelTangent, + α::Real, W′) + return transport_cayley!(copy(Θ), W, Δ, α, W′) +end end diff --git a/src/unitary.jl b/src/unitary.jl index 5073994..ea74042 100644 --- a/src/unitary.jl +++ b/src/unitary.jl @@ -8,28 +8,32 @@ import TensorKit: similarstoragetype, SectorDict using ..TensorKitManifolds: projectantihermitian!, projectisometric!, PolarNewton import ..TensorKitManifolds: base, checkbase, inner, retract, transport, transport! -struct UnitaryTangent{T<:AbstractTensorMap, TA<:AbstractTensorMap} +struct UnitaryTangent{T<:AbstractTensorMap,TA<:AbstractTensorMap} W::T A::TA function UnitaryTangent(W::AbstractTensorMap{S,N₁,N₂}, A::AbstractTensorMap{S,N₂,N₂}) where {S,N₁,N₂} T = typeof(W) TA = typeof(A) - return new{T,TA}(W,A) + return new{T,TA}(W, A) end end Base.copy(Δ::UnitaryTangent) = UnitaryTangent(Δ.W, copy(Δ.A)) Base.getindex(Δ::UnitaryTangent) = Δ.W * Δ.A base(Δ::UnitaryTangent) = Δ.W -checkbase(Δ₁::UnitaryTangent, Δ₂::UnitaryTangent) = Δ₁.W == Δ₂.W ? Δ₁.W : - throw(ArgumentError("tangent vectors with different base points")) +function checkbase(Δ₁::UnitaryTangent, Δ₂::UnitaryTangent) + return Δ₁.W == Δ₂.W ? Δ₁.W : + throw(ArgumentError("tangent vectors with different base points")) +end # Basic vector space behaviour -Base.:+(Δ₁::UnitaryTangent, Δ₂::UnitaryTangent) = - UnitaryTangent(checkbase(Δ₁,Δ₂), Δ₁.A + Δ₂.A) -Base.:-(Δ₁::UnitaryTangent, Δ₂::UnitaryTangent) = - UnitaryTangent(checkbase(Δ₁,Δ₂), Δ₁.A - Δ₂.A) -Base.:-(Δ::UnitaryTangent) = (-1)*Δ +function Base.:+(Δ₁::UnitaryTangent, Δ₂::UnitaryTangent) + return UnitaryTangent(checkbase(Δ₁, Δ₂), Δ₁.A + Δ₂.A) +end +function Base.:-(Δ₁::UnitaryTangent, Δ₂::UnitaryTangent) + return UnitaryTangent(checkbase(Δ₁, Δ₂), Δ₁.A - Δ₂.A) +end +Base.:-(Δ::UnitaryTangent) = (-1) * Δ Base.:*(Δ::UnitaryTangent, α::Real) = rmul!(copy(Δ), α) Base.:*(α::Real, Δ::UnitaryTangent) = lmul!(α, copy(Δ)) @@ -61,34 +65,35 @@ function TensorKit.dot(Δ₁::UnitaryTangent, Δ₂::UnitaryTangent) checkbase(Δ₁, Δ₂) return dot(Δ₁.A, Δ₂.A) end -TensorKit.norm(Δ::UnitaryTangent, p::Real = 2) = norm(Δ.A, p) +TensorKit.norm(Δ::UnitaryTangent, p::Real=2) = norm(Δ.A, p) # tangent space methods function inner(W::AbstractTensorMap, Δ₁::UnitaryTangent, Δ₂::UnitaryTangent; - metric = :euclidean) + metric=:euclidean) @assert metric == :euclidean - Δ₁ === Δ₂ ? norm(Δ₁)^2 : real(dot(Δ₁,Δ₂)) + return Δ₁ === Δ₂ ? norm(Δ₁)^2 : real(dot(Δ₁, Δ₂)) end -function project!(X::AbstractTensorMap, W::AbstractTensorMap; metric = :euclidean) +function project!(X::AbstractTensorMap, W::AbstractTensorMap; metric=:euclidean) @assert metric == :euclidean - P = W'*X + P = W' * X A = projectantihermitian!(P) return UnitaryTangent(W, A) end -project(X, W; metric = :euclidean) = project!(copy(X), W; metric = :euclidean) +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=nothing) W == base(Δ) || throw(ArgumentError("not a valid tangent vector at base point")) - E = exp(α*Δ.A) - W′ = projectisometric!(W*E; alg = SDD()) + E = exp(α * Δ.A) + W′ = projectisometric!(W * E; alg=SDD()) A′ = Δ.A return W′, UnitaryTangent(W′, A′) end -# isometric vector transports compatible with above retraction (also with differential of retraction) +# isometric vector transports compatible with above retraction +# (also with differential of retraction) function transport!(Θ::UnitaryTangent, W::AbstractTensorMap, Δ::UnitaryTangent, α::Real, W′; - alg = :stiefel) + alg=:stiefel) if alg == :parallel return transport_parallel!(Θ, W, Δ, α, W′) elseif alg == :stiefel @@ -98,31 +103,34 @@ function transport!(Θ::UnitaryTangent, W::AbstractTensorMap, Δ::UnitaryTangent end end function transport(Θ::UnitaryTangent, W::AbstractTensorMap, Δ::UnitaryTangent, α::Real, W′; - alg = :stiefel) - return transport!(copy(Θ), W, Δ, α, W′; alg = alg) + alg=:stiefel) + return transport!(copy(Θ), W, Δ, α, W′; alg=alg) end # transport_parallel correspondings to the torsion-free Levi-Civita connection -# transport_stiefel is compatible to Stiefel.transport and corresponds to a non-torsion-free connection -function transport_parallel!(Θ::UnitaryTangent, - W::AbstractTensorMap, - Δ::UnitaryTangent, α, W′) - W == checkbase(Δ,Θ) || throw(ArgumentError("not a valid tangent vector at base point")) - E = exp((α/2)*Δ.A) - A′ = projectantihermitian!(E'*Θ.A*E) # exra projection for stability +# transport_stiefel is compatible to Stiefel.transport and corresponds to a non-torsion-free +# connection +function transport_parallel!(Θ::UnitaryTangent, W::AbstractTensorMap, Δ::UnitaryTangent, α, + W′) + W == checkbase(Δ, Θ) || throw(ArgumentError("not a valid tangent vector at base point")) + E = exp((α / 2) * Δ.A) + A′ = projectantihermitian!(E' * Θ.A * E) # exra projection for stability return UnitaryTangent(W′, A′) end -transport_parallel(Θ::UnitaryTangent, W::AbstractTensorMap, Δ::UnitaryTangent, α, W′) = - transport_parallel!(copy(Θ), W, Δ, α, W′) +function transport_parallel(Θ::UnitaryTangent, W::AbstractTensorMap, Δ::UnitaryTangent, α, + W′) + return transport_parallel!(copy(Θ), W, Δ, α, W′) +end -function transport_stiefel!(Θ::UnitaryTangent, - W::AbstractTensorMap, - Δ::UnitaryTangent, α, W′) - W == checkbase(Δ,Θ) || throw(ArgumentError("not a valid tangent vector at base point")) +function transport_stiefel!(Θ::UnitaryTangent, W::AbstractTensorMap, Δ::UnitaryTangent, + α, W′) + W == checkbase(Δ, Θ) || throw(ArgumentError("not a valid tangent vector at base point")) A′ = Θ.A return UnitaryTangent(W′, A′) end -transport_stiefel(Θ::UnitaryTangent, W::AbstractTensorMap, Δ::UnitaryTangent, α, W′) = - transport_stiefel!(copy(Θ), W, Δ, α, W′) +function transport_stiefel(Θ::UnitaryTangent, W::AbstractTensorMap, Δ::UnitaryTangent, + α, W′) + return transport_stiefel!(copy(Θ), W, Δ, α, W′) +end end diff --git a/test/runtests.jl b/test/runtests.jl index 91e7517..ab9d77f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,41 +1,43 @@ using TensorKit, TensorKitManifolds using Test -spaces = (ℂ^4, ℤ₂Space(2,2), U₁Space(0=>2,1=>1,-1=>1), SU₂Space(0=>2,1/2=>1)) +spaces = (ℂ^4, ℤ₂Space(2, 2), U₁Space(0 => 2, 1 => 1, -1 => 1), + SU₂Space(0 => 2, 1 / 2 => 1)) const ϵ = 1e-7 const α = 0.75 @testset "Grassmann with space $V" for V in spaces for T in (Float64,) - W, = leftorth(TensorMap(randn, T, V*V*V, V*V); alg = Polar()) + W, = leftorth(TensorMap(randn, T, V * V * V, V * V); alg=Polar()) X = TensorMap(randn, T, space(W)) Y = TensorMap(randn, T, space(W)) Δ = @inferred Grassmann.project(X, W) Θ = Grassmann.project(Y, W) γ = randn(T) - Ξ = -Δ + γ*Θ - @test norm(W'*Δ[]) <= sqrt(eps(real(T)))*dim(domain(W)) - @test norm(W'*Θ[]) <= sqrt(eps(real(T)))*dim(domain(W)) - @test norm(W'*Ξ[]) <= sqrt(eps(real(T)))*dim(domain(W)) + Ξ = -Δ + γ * Θ + @test norm(W' * Δ[]) <= sqrt(eps(real(T))) * dim(domain(W)) + @test norm(W' * Θ[]) <= sqrt(eps(real(T))) * dim(domain(W)) + @test norm(W' * Ξ[]) <= sqrt(eps(real(T))) * dim(domain(W)) @test norm(zero(W)) == 0 @test (@inferred Grassmann.inner(W, Δ, Θ)) ≈ real(dot(Δ[], Θ[])) @test Grassmann.inner(W, Δ, Θ) ≈ real(dot(X, Θ[])) - @test Grassmann.inner(W, Δ, Θ) ≈ real(dot(Δ[],Y)) + @test Grassmann.inner(W, Δ, Θ) ≈ real(dot(Δ[], Y)) @test Grassmann.inner(W, Δ, Δ) ≈ norm(Δ[])^2 W2, = @inferred Grassmann.retract(W, Δ, ϵ) - @test W2 ≈ W + ϵ*Δ[] + @test W2 ≈ W + ϵ * Δ[] W2, Δ2′ = Grassmann.retract(W, Δ, α) - @test norm(W2'*Δ2′[]) <= sqrt(eps(real(T)))*dim(domain(W)) - @test Δ2′[] ≈ (first(Grassmann.retract(W, Δ, α + ϵ/2)) - - first(Grassmann.retract(W, Δ, α - ϵ/2)))/(ϵ) + @test norm(W2' * Δ2′[]) <= sqrt(eps(real(T))) * dim(domain(W)) + @test Δ2′[] ≈ + (first(Grassmann.retract(W, Δ, α + ϵ / 2)) - + first(Grassmann.retract(W, Δ, α - ϵ / 2))) / (ϵ) Δ2 = @inferred Grassmann.transport(Δ, W, Δ, α, W2) Θ2 = Grassmann.transport(Θ, W, Δ, α, W2) Ξ2 = Grassmann.transport(Ξ, W, Δ, α, W2) @test Δ2[] ≈ Δ2′[] - @test norm(W2'*Δ2[]) <= sqrt(eps(real(T)))*dim(domain(W)) - @test norm(W2'*Θ2[]) <= sqrt(eps(real(T)))*dim(domain(W)) - @test norm(W2'*Ξ2[]) <= sqrt(eps(real(T)))*dim(domain(W)) + @test norm(W2' * Δ2[]) <= sqrt(eps(real(T))) * dim(domain(W)) + @test norm(W2' * Θ2[]) <= sqrt(eps(real(T))) * dim(domain(W)) + @test norm(W2' * Ξ2[]) <= sqrt(eps(real(T))) * dim(domain(W)) @test Ξ2[] ≈ -Δ2[] + γ * Θ2[] @test Grassmann.inner(W2, Δ2, Θ2) ≈ Grassmann.inner(W, Δ, Θ) @test Grassmann.inner(W2, Ξ2, Θ2) ≈ Grassmann.inner(W, Ξ, Θ) @@ -51,40 +53,42 @@ end @testset "Stiefel with space $V" for V in spaces for T in (Float64, ComplexF64) - W = TensorMap(randhaar, T, V*V*V, V*V) + W = TensorMap(randhaar, T, V * V * V, V * V) X = TensorMap(randn, T, space(W)) Y = TensorMap(randn, T, space(W)) Δ = @inferred Stiefel.project_euclidean(X, W) Θ = Stiefel.project_canonical(Y, W) γ = rand() - Ξ = -Δ + γ*Θ - @test norm(W'*Δ[] + Δ[]'*W) <= sqrt(eps(real(T)))*dim(domain(W)) - @test norm(W'*Θ[] + Θ[]'*W) <= sqrt(eps(real(T)))*dim(domain(W)) - @test norm(W'*Ξ[] + Ξ[]'*W) <= sqrt(eps(real(T)))*dim(domain(W)) + Ξ = -Δ + γ * Θ + @test norm(W' * Δ[] + Δ[]' * W) <= sqrt(eps(real(T))) * dim(domain(W)) + @test norm(W' * Θ[] + Θ[]' * W) <= sqrt(eps(real(T))) * dim(domain(W)) + @test norm(W' * Ξ[] + Ξ[]' * W) <= sqrt(eps(real(T))) * dim(domain(W)) @test norm(zero(W)) == 0 @test (@inferred Stiefel.inner_euclidean(W, Δ, Θ)) ≈ real(dot(Δ[], Θ[])) @test (@inferred Stiefel.inner_canonical(W, Δ, Θ)) ≈ - real(dot(Δ[], Θ[] - W*(W'*Θ[])/2)) + real(dot(Δ[], Θ[] - W * (W' * Θ[]) / 2)) @test Stiefel.inner_euclidean(W, Δ, Θ) ≈ real(dot(X, Θ[])) - @test !(Stiefel.inner_euclidean(W, Δ, Θ) ≈ real(dot(Δ[],Y))) + @test !(Stiefel.inner_euclidean(W, Δ, Θ) ≈ real(dot(Δ[], Y))) @test !(Stiefel.inner_canonical(W, Δ, Θ) ≈ real(dot(X, Θ[]))) - @test Stiefel.inner_canonical(W, Δ, Θ) ≈ real(dot(Δ[],Y)) + @test Stiefel.inner_canonical(W, Δ, Θ) ≈ real(dot(Δ[], Y)) @test Stiefel.inner_euclidean(W, Δ, Δ) ≈ norm(Δ[])^2 - @test Stiefel.inner_canonical(W, Δ, Δ) ≈ (1//2)*norm(W'*Δ[])^2 + norm(Δ[]-W*(W'Δ[]))^2 + @test Stiefel.inner_canonical(W, Δ, Δ) ≈ + (1 // 2) * norm(W' * Δ[])^2 + norm(Δ[] - W * (W'Δ[]))^2 W2, = @inferred Stiefel.retract_exp(W, Δ, ϵ) - @test W2 ≈ W + ϵ*Δ[] + @test W2 ≈ W + ϵ * Δ[] W2, Δ2′ = Stiefel.retract_exp(W, Δ, α) - @test norm(W2'*Δ2′[] + Δ2′[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W)) - @test Δ2′[] ≈ (first(Stiefel.retract_exp(W, Δ, α+ϵ/2)) - - first(Stiefel.retract_exp(W, Δ, α-ϵ/2)))/(ϵ) + @test norm(W2' * Δ2′[] + Δ2′[]' * W2) <= sqrt(eps(real(T))) * dim(domain(W)) + @test Δ2′[] ≈ + (first(Stiefel.retract_exp(W, Δ, α + ϵ / 2)) - + first(Stiefel.retract_exp(W, Δ, α - ϵ / 2))) / (ϵ) Δ2 = @inferred Stiefel.transport_exp(Δ, W, Δ, α, W2) Θ2 = Stiefel.transport_exp(Θ, W, Δ, α, W2) Ξ2 = Stiefel.transport_exp(Ξ, W, Δ, α, W2) @test Δ2′[] ≈ Δ2[] - @test norm(W2'*Δ2[] + Δ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W)) - @test norm(W2'*Θ2[] + Θ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W)) - @test norm(W2'*Ξ2[] + Ξ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W)) + @test norm(W2' * Δ2[] + Δ2[]' * W2) <= sqrt(eps(real(T))) * dim(domain(W)) + @test norm(W2' * Θ2[] + Θ2[]' * W2) <= sqrt(eps(real(T))) * dim(domain(W)) + @test norm(W2' * Ξ2[] + Ξ2[]' * W2) <= sqrt(eps(real(T))) * dim(domain(W)) @test Ξ2[] ≈ -Δ2[] + γ * Θ2[] @test Stiefel.inner_euclidean(W2, Δ2, Θ2) ≈ Stiefel.inner_euclidean(W, Δ, Θ) @test Stiefel.inner_euclidean(W2, Ξ2, Θ2) ≈ Stiefel.inner_euclidean(W, Ξ, Θ) @@ -92,19 +96,20 @@ end @test Stiefel.inner_canonical(W2, Ξ2, Θ2) ≈ Stiefel.inner_canonical(W, Ξ, Θ) W2, = @inferred Stiefel.retract_cayley(W, Δ, ϵ) - @test W2 ≈ W + ϵ*Δ[] + @test W2 ≈ W + ϵ * Δ[] W2, Δ2′ = Stiefel.retract_cayley(W, Δ, α) - @test norm(W2'*Δ2′[] + Δ2′[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W)) - @test Δ2′[] ≈ (first(Stiefel.retract_cayley(W, Δ, α+ϵ/2)) - - first(Stiefel.retract_cayley(W, Δ, α-ϵ/2)))/(ϵ) + @test norm(W2' * Δ2′[] + Δ2′[]' * W2) <= sqrt(eps(real(T))) * dim(domain(W)) + @test Δ2′[] ≈ + (first(Stiefel.retract_cayley(W, Δ, α + ϵ / 2)) - + first(Stiefel.retract_cayley(W, Δ, α - ϵ / 2))) / (ϵ) @test norm(Δ2′) <= norm(Δ) Δ2 = @inferred Stiefel.transport_cayley(Δ, W, Δ, α, W2) Θ2 = Stiefel.transport_cayley(Θ, W, Δ, α, W2) Ξ2 = Stiefel.transport_cayley(Ξ, W, Δ, α, W2) @test !(Δ2′[] ≈ Δ2[]) - @test norm(W2'*Δ2[] + Δ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W)) - @test norm(W2'*Θ2[] + Θ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W)) - @test norm(W2'*Ξ2[] + Ξ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W)) + @test norm(W2' * Δ2[] + Δ2[]' * W2) <= sqrt(eps(real(T))) * dim(domain(W)) + @test norm(W2' * Θ2[] + Θ2[]' * W2) <= sqrt(eps(real(T))) * dim(domain(W)) + @test norm(W2' * Ξ2[] + Ξ2[]' * W2) <= sqrt(eps(real(T))) * dim(domain(W)) @test Ξ2[] ≈ -Δ2[] + γ * Θ2[] @test Stiefel.inner_euclidean(W2, Δ2, Θ2) ≈ Stiefel.inner_euclidean(W, Δ, Θ) @test Stiefel.inner_euclidean(W2, Ξ2, Θ2) ≈ Stiefel.inner_euclidean(W, Ξ, Θ) @@ -119,35 +124,36 @@ end @testset "Unitary with space $V" for V in spaces for T in (Float64, ComplexF64) - W, = leftorth(TensorMap(randn, T, V*V*V, V*V); alg = Polar()) + W, = leftorth(TensorMap(randn, T, V * V * V, V * V); alg=Polar()) X = TensorMap(randn, T, space(W)) Y = TensorMap(randn, T, space(W)) Δ = @inferred Unitary.project(X, W) Θ = Unitary.project(Y, W) γ = randn() - Ξ = -Δ + γ*Θ - @test norm(W'*Δ[] + Δ[]'*W) <= sqrt(eps(real(T)))*dim(domain(W)) - @test norm(W'*Θ[] + Θ[]'*W) <= sqrt(eps(real(T)))*dim(domain(W)) - @test norm(W'*Ξ[] + Ξ[]'*W) <= sqrt(eps(real(T)))*dim(domain(W)) + Ξ = -Δ + γ * Θ + @test norm(W' * Δ[] + Δ[]' * W) <= sqrt(eps(real(T))) * dim(domain(W)) + @test norm(W' * Θ[] + Θ[]' * W) <= sqrt(eps(real(T))) * dim(domain(W)) + @test norm(W' * Ξ[] + Ξ[]' * W) <= sqrt(eps(real(T))) * dim(domain(W)) @test norm(zero(W)) == 0 @test (@inferred Unitary.inner(W, Δ, Θ)) ≈ real(dot(Δ[], Θ[])) @test Unitary.inner(W, Δ, Θ) ≈ real(dot(X, Θ[])) - @test Unitary.inner(W, Δ, Θ) ≈ real(dot(Δ[],Y)) + @test Unitary.inner(W, Δ, Θ) ≈ real(dot(Δ[], Y)) @test Unitary.inner(W, Δ, Δ) ≈ norm(Δ[])^2 W2, = @inferred Unitary.retract(W, Δ, ϵ) - @test W2 ≈ W + ϵ*Δ[] + @test W2 ≈ W + ϵ * Δ[] W2, Δ2′ = Unitary.retract(W, Δ, α) - @test norm(W2'*Δ2′[] + Δ2′[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W)) - @test Δ2′[] ≈ (first(Unitary.retract(W, Δ, α+ϵ/2)) - - first(Unitary.retract(W, Δ, α-ϵ/2)))/(ϵ) + @test norm(W2' * Δ2′[] + Δ2′[]' * W2) <= sqrt(eps(real(T))) * dim(domain(W)) + @test Δ2′[] ≈ + (first(Unitary.retract(W, Δ, α + ϵ / 2)) - + first(Unitary.retract(W, Δ, α - ϵ / 2))) / (ϵ) Δ2 = @inferred Unitary.transport_parallel(Δ, W, Δ, α, W2) Θ2 = Unitary.transport_parallel(Θ, W, Δ, α, W2) Ξ2 = Unitary.transport_parallel(Ξ, W, Δ, α, W2) @test Δ2′[] ≈ Δ2[] - @test norm(W2'*Θ2[] + Θ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W)) - @test norm(W2'*Ξ2[] + Ξ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W)) + @test norm(W2' * Θ2[] + Θ2[]' * W2) <= sqrt(eps(real(T))) * dim(domain(W)) + @test norm(W2' * Ξ2[] + Ξ2[]' * W2) <= sqrt(eps(real(T))) * dim(domain(W)) @test Ξ2[] ≈ -Δ2[] + γ * Θ2[] @test Unitary.inner(W2, Δ2, Θ2) ≈ Unitary.inner(W, Δ, Θ) @test Unitary.inner(W2, Ξ2, Θ2) ≈ Unitary.inner(W, Ξ, Θ) @@ -156,9 +162,9 @@ end Θ2 = Unitary.transport_stiefel(Θ, W, Δ, α, W2) Ξ2 = Unitary.transport_stiefel(Ξ, W, Δ, α, W2) @test Δ2′[] ≈ Δ2[] - @test norm(W2'*Δ2[] + Δ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W)) - @test norm(W2'*Θ2[] + Θ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W)) - @test norm(W2'*Ξ2[] + Ξ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W)) + @test norm(W2' * Δ2[] + Δ2[]' * W2) <= sqrt(eps(real(T))) * dim(domain(W)) + @test norm(W2' * Θ2[] + Θ2[]' * W2) <= sqrt(eps(real(T))) * dim(domain(W)) + @test norm(W2' * Ξ2[] + Ξ2[]' * W2) <= sqrt(eps(real(T))) * dim(domain(W)) @test Ξ2[] ≈ -Δ2[] + γ * Θ2[] @test Unitary.inner(W2, Δ2, Θ2) ≈ Unitary.inner(W, Δ, Θ) @test Unitary.inner(W2, Ξ2, Θ2) ≈ Unitary.inner(W, Ξ, Θ) From 76acaed3c5f238033800817f9eb55adcb5523f89 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Fri, 1 Sep 2023 14:14:24 +0200 Subject: [PATCH 4/9] Update to TensorKit v0.11 --- Project.toml | 2 +- src/TensorKitManifolds.jl | 4 ++-- src/grassmann.jl | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index ec25f9c..42d4b52 100644 --- a/Project.toml +++ b/Project.toml @@ -8,7 +8,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" [compat] -TensorKit = "0.6,0.7,0.8,0.9,0.10,1" +TensorKit = "0.11" julia = "1" [extras] diff --git a/src/TensorKitManifolds.jl b/src/TensorKitManifolds.jl index 6ca33be..0316413 100644 --- a/src/TensorKitManifolds.jl +++ b/src/TensorKitManifolds.jl @@ -26,11 +26,11 @@ function checkbase end 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 -eleps(X) = eps(real(eltype(X))) +eleps(X) = eps(real(scalartype(X))) function isisometry(W::AbstractTensorMap; tol=10 * eleps(W)) WdW = W' * W - s = zero(float(real(eltype(W)))) + s = zero(float(real(scalartype(W)))) for (c, b) in blocks(WdW) _subtractone!(b) s += dim(c) * length(b) diff --git a/src/grassmann.jl b/src/grassmann.jl index 0243fac..c96cabb 100644 --- a/src/grassmann.jl +++ b/src/grassmann.jl @@ -24,7 +24,7 @@ mutable struct GrassmannTangent{T<:AbstractTensorMap, function GrassmannTangent(W::AbstractTensorMap{S,N₁,N₂}, Z::AbstractTensorMap{S,N₁,N₂}) where {S,N₁,N₂} T = typeof(W) - TT = promote_type(float(eltype(W)), eltype(Z)) + TT = promote_type(float(scalartype(W)), scalartype(Z)) G = sectortype(W) M = similarstoragetype(W, TT) Mr = similarstoragetype(W, real(TT)) @@ -198,7 +198,7 @@ function invretract(Wold::AbstractTensorMap, Wnew::AbstractTensorMap; alg=nothin Vd, cS, VY = tsvd!(WodWn) Scmplx = acos(cS) # acos always returns a complex TensorMap. We cast back to real if possible. - S = eltype(WodWn) <: Real && isreal(sectortype(Scmplx)) ? real(Scmplx) : Scmplx + S = scalartype(WodWn) <: Real && isreal(sectortype(Scmplx)) ? real(Scmplx) : Scmplx UsS = Wneworth * VY' # U * sin(S) # should be in polar decomposition form U = projectisometric!(UsS; alg=Polar()) Y = Vd * VY From 0c77e11666eed7dd27b1c7ac8fc454d8732f766f Mon Sep 17 00:00:00 2001 From: lkdvos Date: Fri, 1 Sep 2023 14:15:48 +0200 Subject: [PATCH 5/9] remove unused variable --- src/grassmann.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/grassmann.jl b/src/grassmann.jl index c96cabb..0f63c6f 100644 --- a/src/grassmann.jl +++ b/src/grassmann.jl @@ -25,7 +25,6 @@ mutable struct GrassmannTangent{T<:AbstractTensorMap, Z::AbstractTensorMap{S,N₁,N₂}) where {S,N₁,N₂} T = typeof(W) TT = promote_type(float(scalartype(W)), scalartype(Z)) - G = sectortype(W) M = similarstoragetype(W, TT) Mr = similarstoragetype(W, real(TT)) TU = tensormaptype(S, N₁, 1, M) From 3578b7dbe15c1f1f953a19ef4a851affa86e19f2 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Fri, 1 Sep 2023 14:31:02 +0200 Subject: [PATCH 6/9] Update CI actions --- .github/workflows/ci.yml | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 16f5faf..fdb3e17 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -21,15 +21,17 @@ jobs: arch: - x64 steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v3 - uses: julia-actions/setup-julia@v1 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} + - uses: julia-actions/cache@v1 - uses: julia-actions/julia-buildpkg@latest - uses: julia-actions/julia-runtest@latest env: JULIA_NUM_THREADS: 4 - - uses: julia-actions/julia-uploadcodecov@latest - env: - CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v3 + with: + file: lcov.info From 1e9f43594893e4637c378c5376f40c3a9acf6947 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sat, 2 Sep 2023 10:20:16 +0200 Subject: [PATCH 7/9] Set minimal Julia version to 1.6 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 42d4b52..83db012 100644 --- a/Project.toml +++ b/Project.toml @@ -9,7 +9,7 @@ TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" [compat] TensorKit = "0.11" -julia = "1" +julia = "1.6" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From cffbf312faf31ec3561c3805d237180814b04339 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sat, 2 Sep 2023 10:24:27 +0200 Subject: [PATCH 8/9] `eleps` -> `scalareps` --- src/TensorKitManifolds.jl | 10 +++++----- src/auxiliary.jl | 4 ++-- src/stiefel.jl | 4 ++-- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/TensorKitManifolds.jl b/src/TensorKitManifolds.jl index 0316413..793de7e 100644 --- a/src/TensorKitManifolds.jl +++ b/src/TensorKitManifolds.jl @@ -26,9 +26,9 @@ function checkbase end 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 -eleps(X) = eps(real(scalartype(X))) +scalareps(X) = eps(real(scalartype(X))) -function isisometry(W::AbstractTensorMap; tol=10 * eleps(W)) +function isisometry(W::AbstractTensorMap; tol=10 * scalareps(W)) WdW = W' * W s = zero(float(real(scalartype(W)))) for (c, b) in blocks(WdW) @@ -38,7 +38,7 @@ function isisometry(W::AbstractTensorMap; tol=10 * eleps(W)) return norm(WdW) <= tol * sqrt(s) end -function isunitary(W::AbstractTensorMap; tol=10 * eleps(W)) +function isunitary(W::AbstractTensorMap; tol=10 * scalareps(W)) return isisometry(W; tol=tol) && isisometry(W'; tol=tol) end @@ -81,7 +81,7 @@ function projectisometric!(W::AbstractTensorMap; alg=Polar()) end function projectcomplement!(X::AbstractTensorMap, W::AbstractTensorMap; - tol=10 * eleps(X)) + tol=10 * scalareps(X)) P = W' * X nP = norm(P) nX = norm(X) @@ -102,7 +102,7 @@ function projectisometric(W::AbstractTensorMap; return projectisometric!(copy(W); alg=alg) end function projectcomplement(X::AbstractTensorMap, W::AbstractTensorMap, - tol=10 * eleps(X)) + tol=10 * scalareps(X)) return projectcomplement!(copy(X), W; tol=tol) end diff --git a/src/auxiliary.jl b/src/auxiliary.jl index 5ba8787..f0dea33 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -88,7 +88,7 @@ function _polarsvd!(A::StridedMatrix) U, S, V = svd!(A; alg=LinearAlgebra.QRIteration()) return mul!(A, U, V') end -function _polarnewton!(A::StridedMatrix; tol=10 * eleps(A), maxiter=5) +function _polarnewton!(A::StridedMatrix; tol=10 * scalareps(A), maxiter=5) m, n = size(A) @assert m >= n A2 = copy(A) @@ -157,7 +157,7 @@ function _stiefelexp(W::StridedMatrix, A::StridedMatrix, Z::StridedMatrix, α) end function _stiefellog(Wold::StridedMatrix, Wnew::StridedMatrix; - tol=10 * eleps(Wold), maxiter=100) + tol=10 * scalareps(Wold), maxiter=100) n, p = size(Wold) r = min(2 * p, n) P = Wold' * Wnew diff --git a/src/stiefel.jl b/src/stiefel.jl index bd28f0e..16c216a 100644 --- a/src/stiefel.jl +++ b/src/stiefel.jl @@ -8,7 +8,7 @@ using TensorKit using TensorKit: similarstoragetype, SectorDict using ..TensorKitManifolds: projecthermitian!, projectantihermitian!, projectisometric!, projectcomplement!, PolarNewton, - _stiefelexp, _stiefellog, eleps + _stiefelexp, _stiefellog, scalareps import ..TensorKitManifolds: base, checkbase, inner, retract, transport, transport! @@ -199,7 +199,7 @@ function retract_exp(W::AbstractTensorMap, Δ::StiefelTangent, α::Real) return W′, StiefelTangent(W′, A′, Z′) end function invretract_exp(Wold::AbstractTensorMap, Wnew::AbstractTensorMap; - tol=eleps(Wold)^(2 / 3)) + tol=scalareps(Wold)^(2 / 3)) space(Wold) == space(Wnew) || throw(SectorMismatch()) S = spacetype(Wold) From 08bce532c4d3cad74a96cea4a68e8cc9cd087c37 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sat, 2 Sep 2023 10:25:00 +0200 Subject: [PATCH 9/9] Recycle `_one!` from TensorKit --- src/auxiliary.jl | 16 +++------------- 1 file changed, 3 insertions(+), 13 deletions(-) diff --git a/src/auxiliary.jl b/src/auxiliary.jl index f0dea33..8f391fb 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -70,16 +70,6 @@ function _subtractone!(a::AbstractMatrix) view(a, diagind(a)) .= view(a, diagind(a)) .- 1 return a end -function _one!(A::AbstractMatrix) - m, n = size(A) - @inbounds for j in 1:n - for i in 1:m - A[i, j] = i == j - end - end - return A -end - function _polarsdd!(A::StridedMatrix) U, S, V = svd!(A; alg=LinearAlgebra.DivideAndConquer()) return mul!(A, U, V') @@ -93,7 +83,7 @@ function _polarnewton!(A::StridedMatrix; tol=10 * scalareps(A), maxiter=5) @assert m >= n A2 = copy(A) Q, R = qr!(A2) - Ri = ldiv!(UpperTriangular(R)', _one!(similar(R))) + Ri = ldiv!(UpperTriangular(R)', TensorKit._one!(similar(R))) R, Ri = _avgdiff!(R, Ri) i = 1 R2 = view(A, 1:n, 1:n) @@ -104,7 +94,7 @@ function _polarnewton!(A::StridedMatrix; tol=10 * scalareps(A), maxiter=5) _polarsdd!(Ri) break end - Ri = ldiv!(lu!(R2)', _one!(Ri)) + Ri = ldiv!(lu!(R2)', TensorKit._one!(Ri)) R, Ri = _avgdiff!(R, Ri) copyto!(R2, R) i += 1 @@ -173,7 +163,7 @@ function _stiefellog(Wold::StridedMatrix, Wnew::StridedMatrix; R = Q' * dW Wext = [Wold Q] F = qr!([P; R]) - U = lmul!(F.Q, _one!(similar(P, r, r))) + U = lmul!(F.Q, TensorKit._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)