From 26bc53aed9925eb40a63d32952beeb88df59258d Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sun, 12 Feb 2023 11:57:47 +0100 Subject: [PATCH 1/4] change spacetype hierarchy to traits --- src/TensorKit.jl | 3 +- src/spaces/cartesianspace.jl | 8 +- src/spaces/complexspace.jl | 36 +++++---- src/spaces/deligne.jl | 15 ++-- src/spaces/generalspace.jl | 8 +- src/spaces/gradedspace.jl | 36 +++++---- src/spaces/homspace.jl | 6 +- src/spaces/vectorspaces.jl | 39 +++++----- src/tensors/factorizations.jl | 142 +++++++++++++++++++++------------- src/tensors/linalg.jl | 75 +++++++++++------- test/spaces.jl | 16 ++-- 11 files changed, 228 insertions(+), 156 deletions(-) diff --git a/src/TensorKit.jl b/src/TensorKit.jl index 414a0f3e..8320e6b0 100644 --- a/src/TensorKit.jl +++ b/src/TensorKit.jl @@ -17,7 +17,8 @@ export Fermion, FermionParity, FermionNumber, FermionSpin export FibonacciAnyon export IsingAnyon -export VectorSpace, Field, ElementarySpace, InnerProductSpace, EuclideanSpace # abstract vector spaces +export VectorSpace, Field, ElementarySpace # abstract vector spaces +export InnerProductStyle, NoInnerProduct, HasInnerProduct, EuclideanProduct export ComplexSpace, CartesianSpace, GeneralSpace, GradedSpace # concrete spaces export ZNSpace, Z2Space, Z3Space, Z4Space, U1Space, CU1Space, SU2Space export Vect, Rep # space constructors diff --git a/src/spaces/cartesianspace.jl b/src/spaces/cartesianspace.jl index 03bf5884..d46b9bd0 100644 --- a/src/spaces/cartesianspace.jl +++ b/src/spaces/cartesianspace.jl @@ -1,11 +1,11 @@ """ - struct CartesianSpace <: EuclideanSpace{ℝ} + struct CartesianSpace <: ElementarySpace{ℝ} A real Euclidean space `ℝ^d`, which is therefore self-dual. `CartesianSpace` has no additonal structure and is completely characterised by its dimension `d`. This is the vector space that is implicitly assumed in most of matrix algebra. """ -struct CartesianSpace <: EuclideanSpace{ℝ} +struct CartesianSpace <: ElementarySpace{ℝ} d::Int end CartesianSpace(d::Integer = 0; dual = false) = CartesianSpace(Int(d)) @@ -28,6 +28,10 @@ function CartesianSpace(dims::AbstractDict; kwargs...) end end +InnerProductStyle(::Type{CartesianSpace}) = EuclideanProduct() + +isdual(V::CartesianSpace) = false + # convenience constructor Base.getindex(::RealNumbers) = CartesianSpace Base.:^(::RealNumbers, d::Int) = CartesianSpace(d) diff --git a/src/spaces/complexspace.jl b/src/spaces/complexspace.jl index afad6586..7f443d77 100644 --- a/src/spaces/complexspace.jl +++ b/src/spaces/complexspace.jl @@ -1,13 +1,13 @@ """ - struct ComplexSpace <: EuclideanSpace{ℂ} + struct ComplexSpace <: ElementarySpace{ℂ} A standard complex vector space ℂ^d with Euclidean inner product and no additional structure. It is completely characterised by its dimension and whether its the normal space or its dual (which is canonically isomorphic to the conjugate space). """ -struct ComplexSpace <: EuclideanSpace{ℂ} - d::Int - dual::Bool +struct ComplexSpace <: ElementarySpace{ℂ} + d::Int + dual::Bool end ComplexSpace(d::Integer = 0; dual = false) = ComplexSpace(Int(d), dual) function ComplexSpace(dim::Pair; dual = false) @@ -29,6 +29,8 @@ function ComplexSpace(dims::AbstractDict; kwargs...) end end +InnerProductStyle(::Type{ComplexSpace}) = EuclideanProduct() + # convenience constructor Base.getindex(::ComplexNumbers) = ComplexSpace Base.:^(::ComplexNumbers, d::Int) = ComplexSpace(d) @@ -42,17 +44,23 @@ Base.axes(V::ComplexSpace) = Base.OneTo(dim(V)) Base.conj(V::ComplexSpace) = ComplexSpace(dim(V), !isdual(V)) Base.oneunit(::Type{ComplexSpace}) = ComplexSpace(1) -⊕(V1::ComplexSpace, V2::ComplexSpace) = isdual(V1) == isdual(V2) ? - ComplexSpace(dim(V1)+dim(V2), isdual(V1)) : - throw(SpaceMismatch("Direct sum of a vector space and its dual does not exist")) -fuse(V1::ComplexSpace, V2::ComplexSpace) = ComplexSpace(V1.d*V2.d) +function ⊕(V1::ComplexSpace, V2::ComplexSpace) + return isdual(V1) == isdual(V2) ? + ComplexSpace(dim(V1) + dim(V2), isdual(V1)) : + throw(SpaceMismatch("Direct sum of a vector space and its dual does not exist")) +end +fuse(V1::ComplexSpace, V2::ComplexSpace) = ComplexSpace(V1.d * V2.d) flip(V::ComplexSpace) = dual(V) -infimum(V1::ComplexSpace, V2::ComplexSpace) = isdual(V1) == isdual(V2) ? - ComplexSpace(min(dim(V1), dim(V2)), isdual(V1)) : - throw(SpaceMismatch("Infimum of space and dual space does not exist")) -supremum(V1::ComplexSpace, V2::ComplexSpace) = isdual(V1) == isdual(V2) ? - ComplexSpace(max(dim(V1), dim(V2)), isdual(V1)) : - throw(SpaceMismatch("Supremum of space and dual space does not exist")) +function infimum(V1::ComplexSpace, V2::ComplexSpace) + return isdual(V1) == isdual(V2) ? + ComplexSpace(min(dim(V1), dim(V2)), isdual(V1)) : + throw(SpaceMismatch("Infimum of space and dual space does not exist")) +end +function supremum(V1::ComplexSpace, V2::ComplexSpace) + return isdual(V1) == isdual(V2) ? + ComplexSpace(max(dim(V1), dim(V2)), isdual(V1)) : + throw(SpaceMismatch("Supremum of space and dual space does not exist")) +end Base.show(io::IO, V::ComplexSpace) = print(io, isdual(V) ? "(ℂ^$(V.d))'" : "ℂ^$(V.d)") diff --git a/src/spaces/deligne.jl b/src/spaces/deligne.jl index 9ff9d9e5..4a321837 100644 --- a/src/spaces/deligne.jl +++ b/src/spaces/deligne.jl @@ -15,35 +15,36 @@ natural representation spaces of the direct product of two groups. ⊠(V1::VectorSpace, V2::VectorSpace) = (V1 ⊠ one(V2)) ⊗ (one(V1) ⊠ V2) # define deligne products with empty tensor product: just add a trivial sector of the type of the empty space to each of the sectors in the non-empty space -function ⊠(V::GradedSpace, P0::ProductSpace{<:EuclideanSpace{ℂ},0}) +function ⊠(V::GradedSpace, P0::ProductSpace{<:ElementarySpace{ℂ},0}) I1 = sectortype(V) I2 = sectortype(P0) return Vect[I1 ⊠ I2](ifelse(isdual(V), dual(c), c) ⊠ one(I2) => dim(V, c) for c in sectors(V); dual = isdual(V)) end -function ⊠(P0::ProductSpace{<:EuclideanSpace{ℂ},0}, V::GradedSpace) +function ⊠(P0::ProductSpace{<:ElementarySpace{ℂ},0}, V::GradedSpace) I1 = sectortype(P0) I2 = sectortype(V) return Vect[I1 ⊠ I2](one(I1) ⊠ ifelse(isdual(V), dual(c), c) => dim(V, c) for c in sectors(V); dual = isdual(V)) end -function ⊠(V::ComplexSpace, P0::ProductSpace{<:EuclideanSpace{ℂ},0}) +function ⊠(V::ComplexSpace, P0::ProductSpace{<:ElementarySpace{ℂ},0}) I2 = sectortype(P0) return Vect[I2](one(I2) => dim(V); dual = isdual(V)) end -function ⊠(P0::ProductSpace{<:EuclideanSpace{ℂ},0}, V::ComplexSpace) +function ⊠(P0::ProductSpace{<:ElementarySpace{ℂ},0}, V::ComplexSpace) I1 = sectortype(P0) return Vect[I1](one(I1) => dim(V); dual = isdual(V)) end -function ⊠(P::ProductSpace{<:EuclideanSpace{ℂ},0}, P0::ProductSpace{<:EuclideanSpace{ℂ},0}) +function ⊠(P::ProductSpace{<:ElementarySpace{ℂ},0}, + P0::ProductSpace{<:ElementarySpace{ℂ},0}) I1 = sectortype(P) I2 = sectortype(P0) return one(Vect[I1 ⊠ I2]) end -function ⊠(P::ProductSpace{<:EuclideanSpace{ℂ}}, P0::ProductSpace{<:EuclideanSpace{ℂ},0}) +function ⊠(P::ProductSpace{<:ElementarySpace{ℂ}}, P0::ProductSpace{<:ElementarySpace{ℂ},0}) I1 = sectortype(P) I2 = sectortype(P0) S = Vect[I1 ⊠ I2] @@ -51,7 +52,7 @@ function ⊠(P::ProductSpace{<:EuclideanSpace{ℂ}}, P0::ProductSpace{<:Euclidea return ProductSpace{S,N}(map(V->V ⊠ P0, tuple(P...))) end -function ⊠(P0::ProductSpace{<:EuclideanSpace{ℂ},0}, P::ProductSpace{<:EuclideanSpace{ℂ}}) +function ⊠(P0::ProductSpace{<:ElementarySpace{ℂ},0}, P::ProductSpace{<:ElementarySpace{ℂ}}) I1 = sectortype(P0) I2 = sectortype(P) S = Vect[I1 ⊠ I2] diff --git a/src/spaces/generalspace.jl b/src/spaces/generalspace.jl index 5dd940b7..0f3268fd 100644 --- a/src/spaces/generalspace.jl +++ b/src/spaces/generalspace.jl @@ -28,10 +28,10 @@ isconj(V::GeneralSpace) = V.conj Base.axes(V::GeneralSpace) = Base.OneTo(dim(V)) -dual(V::GeneralSpace{𝕜}) where {𝕜} = - GeneralSpace{𝕜}(dim(V), !isdual(V), isconj(V)) -Base.conj(V::GeneralSpace{𝕜}) where {𝕜} = - GeneralSpace{𝕜}(dim(V), isdual(V), !isconj(V)) +InnerProductStyle(::Type{<:GeneralSpace}) = NoInnerProduct() + +dual(V::GeneralSpace{𝕜}) where {𝕜} = GeneralSpace{𝕜}(dim(V), !isdual(V), isconj(V)) +Base.conj(V::GeneralSpace{𝕜}) where {𝕜} = GeneralSpace{𝕜}(dim(V), isdual(V), !isconj(V)) function Base.show(io::IO, V::GeneralSpace{𝕜}) where {𝕜} if isconj(V) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index 67fb8da3..3489da4b 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -1,5 +1,5 @@ """ - struct GradedSpace{I<:Sector, D} <: EuclideanSpace{ℂ} + struct GradedSpace{I<:Sector, D} <: ElementarySpace{ℂ} dims::D dual::Bool end @@ -24,7 +24,7 @@ and should typically be of no concern. The concrete type `GradedSpace{I,D}` with correct `D` can be obtained as `Vect[I]`, or if `I == Irrep[G]` for some `G<:Group`, as `Rep[G]`. """ -struct GradedSpace{I<:Sector, D} <: EuclideanSpace{ℂ} +struct GradedSpace{I<:Sector,D} <: ElementarySpace{ℂ} dims::D dual::Bool end @@ -78,20 +78,26 @@ Base.hash(V::GradedSpace, h::UInt) = hash(V.dual, hash(V.dims, h)) # Corresponding methods: # properties -dim(V::GradedSpace) = - reduce(+, dim(V, c) * dim(c) for c in sectors(V); init = zero(dim(one(sectortype(V))))) - -dim(V::GradedSpace{I,<:AbstractDict}, c::I) where {I<:Sector} = - get(V.dims, isdual(V) ? dual(c) : c, 0) -dim(V::GradedSpace{I,<:Tuple}, c::I) where {I<:Sector} = - V.dims[findindex(values(I), isdual(V) ? dual(c) : c)] - -sectors(V::GradedSpace{I,<:AbstractDict}) where {I<:Sector} = - SectorSet{I}(s->isdual(V) ? dual(s) : s, keys(V.dims)) -sectors(V::GradedSpace{I,NTuple{N,Int}}) where {I<:Sector, N} = - SectorSet{I}(Iterators.filter(n->V.dims[n]!=0, 1:N)) do n - isdual(V) ? dual(values(I)[n]) : values(I)[n] +InnerProductStyle(::Type{<:GradedSpace}) = EuclideanProduct() +function dim(V::GradedSpace) + return reduce(+, dim(V, c) * dim(c) for c in sectors(V); + init=zero(dim(one(sectortype(V))))) +end +function dim(V::GradedSpace{I,<:AbstractDict}, c::I) where {I<:Sector} + return get(V.dims, isdual(V) ? dual(c) : c, 0) +end +function dim(V::GradedSpace{I,<:Tuple}, c::I) where {I<:Sector} + return V.dims[findindex(values(I), isdual(V) ? dual(c) : c)] +end + +function sectors(V::GradedSpace{I,<:AbstractDict}) where {I<:Sector} + return SectorSet{I}(s -> isdual(V) ? dual(s) : s, keys(V.dims)) +end +function sectors(V::GradedSpace{I,NTuple{N,Int}}) where {I<:Sector,N} + SectorSet{I}(Iterators.filter(n -> V.dims[n] != 0, 1:N)) do n + return isdual(V) ? dual(values(I)[n]) : values(I)[n] end +end hassector(V::GradedSpace{I}, s::I) where {I<:Sector} = dim(V, s) != 0 diff --git a/src/spaces/homspace.jl b/src/spaces/homspace.jl index 05a857d4..bf113767 100644 --- a/src/spaces/homspace.jl +++ b/src/spaces/homspace.jl @@ -16,7 +16,11 @@ codomain(W::HomSpace) = W.codomain domain(W::HomSpace) = W.domain dual(W::HomSpace) = HomSpace(dual(W.domain), dual(W.codomain)) -Base.adjoint(W::HomSpace{<:EuclideanSpace}) = HomSpace(W.domain, W.codomain) +function Base.adjoint(W::HomSpace{S}) where {S} + InnerProductStyle(S) === EuclideanProduct() || + throw(ArgumentError("adjoint requires Euclidean inner product")) + return HomSpace(W.domain, W.codomain) +end Base.hash(W::HomSpace, h::UInt) = hash(domain(W), hash(codomain(W), h)) Base.:(==)(W1::HomSpace, W2::HomSpace) = diff --git a/src/spaces/vectorspaces.jl b/src/spaces/vectorspaces.jl index 5c64e4c6..8d2fd1e6 100644 --- a/src/spaces/vectorspaces.jl +++ b/src/spaces/vectorspaces.jl @@ -115,7 +115,6 @@ that this is different from `one(V::S)`, which returns the empty product space """ Base.oneunit(V::ElementarySpace) = oneunit(typeof(V)) - """ ⊕(V1::S, V2::S, V3::S...) where {S<:ElementarySpace} -> S @@ -126,7 +125,6 @@ spaces `V1`, `V2`, ... Note that all the individual spaces should have the same function ⊕ end ⊕(V1, V2, V3, V4...) = ⊕(⊕(V1, V2), V3, V4...) - """ ⊗(V1::S, V2::S, V3::S...) where {S<:ElementarySpace} -> S @@ -174,30 +172,27 @@ For `field(V)==ℝ`, `conj(V) == V`. It is assumed that `typeof(V) == typeof(con """ Base.conj(V::ElementarySpace{ℝ}) = V -""" - abstract type InnerProductSpace{𝕜} <: ElementarySpace{𝕜} end +# trait to describe the inner product type of vector spaces +abstract type InnerProductStyle end +struct NoInnerProduct <: InnerProductStyle end # no inner product -Abstract type for denoting vector spaces with an inner product, thus a canonical mapping from -`dual(V)` to `V` (for `𝕜 ⊆ ℝ`) or from `dual(V)` to `conj(V)` (otherwise). This mapping -is provided by the metric, but no further support for working with metrics is currently -implemented. -""" -abstract type InnerProductSpace{𝕜} <: ElementarySpace{𝕜} end +abstract type HasInnerProduct <: InnerProductStyle end # inner product defined +struct EuclideanProduct <: HasInnerProduct end # euclidean inner product """ - abstract type EuclideanSpace{𝕜} <: InnerProductSpace{𝕜} end + InnerProductStyle(V::VectorSpace) -> ::InnerProductStyle + InnerProductStyle(S::Type{<:VectorSpace}) -> ::InnerProductStyle -Abstract type for denoting real or complex spaces with a standard Euclidean inner product -(i.e. orthonormal basis, and the metric is identity), such that the dual space is naturally -isomorphic to the conjugate space `dual(V) == conj(V)` (in the complex case) or even to -the space itself `dual(V) == V` (in the real case), also known as the category of -finite-dimensional Hilbert spaces ``FdHilb``. In the language of categories, this subtype -represents dagger or unitary categories, and support an adjoint operation. +Return the type of inner product for vector spaces, which can be either +* `NoInnerProduct()`: no mapping from `dual(V)` to `conj(V)`, i.e. no metric +* subtype of `HasInnerProduct`: a metric exists, but no further support is implemented. +* `EuclideanProduct()`: the metric is the identity, such that dual and conjugate spaces are isomorphic. """ -abstract type EuclideanSpace{𝕜} <: InnerProductSpace{𝕜} end # 𝕜 should be ℝ or ℂ +InnerProductStyle(V::VectorSpace) = InnerProductStyle(typeof(V)) +InnerProductStyle(::Type{<:VectorSpace}) = NoInnerProduct() -dual(V::EuclideanSpace) = conj(V) -isdual(V::EuclideanSpace{ℝ}) = false +dual(V::VectorSpace) = dual(InnerProductStyle(V), V) +dual(::EuclideanProduct, V::VectorSpace) = conj(V) """ sectortype(a) -> Type{<:Sector} @@ -248,9 +243,11 @@ vector spaces of a homogeneous type `S<:ElementarySpace{𝕜}`. """ abstract type CompositeSpace{S<:ElementarySpace} <: VectorSpace end +InnerProductStyle(::Type{<:CompositeSpace{S}}) where {S} = InnerProductStyle(S) + spacetype(S::Type{<:ElementarySpace}) = S spacetype(V::ElementarySpace) = typeof(V) # = spacetype(typeof(V)) -spacetype(::Type{<:CompositeSpace{S}}) where S = S +spacetype(::Type{<:CompositeSpace{S}}) where {S} = S spacetype(V::CompositeSpace) = spacetype(typeof(V)) # = spacetype(typeof(V)) field(P::Type{<:CompositeSpace}) = field(spacetype(P)) diff --git a/src/tensors/factorizations.jl b/src/tensors/factorizations.jl index d751d2f6..5c9e4283 100644 --- a/src/tensors/factorizations.jl +++ b/src/tensors/factorizations.jl @@ -47,8 +47,8 @@ singular values that were truncated. THe keyword `alg` can be equal to `SVD()` or `SDD()`, corresponding to the underlying LAPACK algorithm that computes the decomposition (`_gesvd` or `_gesdd`). -Orthogonality requires `spacetype(t)<:InnerProductSpace`, and `svd(!)` is currently -only implemented for `spacetype(t)<:EuclideanSpace`. +Orthogonality requires `InnerProductStyle(spacetype(t)) <: HasInnerProduct`, and `tsvd(!)` +is currently only implemented for `InnerProductStyle(spacetype(t)) === EuclideanProduct()`. """ tsvd(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; kwargs...) = tsvd!(permute(t, p1, p2; copy = true); kwargs...) @@ -71,8 +71,9 @@ standard QR decomposition such that the diagonal elements of `R` are positive. O `QRpos()` and `Polar()` are uniqe (no residual freedom) so that they always return the same result for the same input tensor `t`. -Orthogonality requires `spacetype(t)<:InnerProductSpace`, and `leftorth(!)` is currently -only implemented for `spacetype(t)<:EuclideanSpace`. +Orthogonality requires `InnerProductStyle(spacetype(t)) <: HasInnerProduct`, and +`leftorth(!)` is currently only implemented for + `InnerProductStyle(spacetype(t)) === EuclideanProduct()`. """ leftorth(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; kwargs...) = leftorth!(permute(t, p1, p2; copy = true); kwargs...) @@ -97,8 +98,9 @@ diagonal elements of `L` are positive. `Polar()` produces a Hermitian and positi semidefinite `L`. Only `LQpos()`, `RQpos()` and `Polar()` are uniqe (no residual freedom) so that they always return the same result for the same input tensor `t`. -Orthogonality requires `spacetype(t)<:InnerProductSpace`, and `rightorth(!)` is currently -only implemented for `spacetype(t)<:EuclideanSpace`. +Orthogonality requires `InnerProductStyle(spacetype(t)) <: HasInnerProduct`, and +`rightorth(!)` is currently only implemented for +`InnerProductStyle(spacetype(t)) === EuclideanProduct()`. """ rightorth(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; kwargs...) = rightorth!(permute(t, p1, p2; copy = true); kwargs...) @@ -121,8 +123,9 @@ value decomposition, and one can specify an absolute or relative tolerance for w singular values are to be considered zero, where `max(atol, norm(t)*rtol)` is used as upper bound. -Orthogonality requires `spacetype(t)<:InnerProductSpace`, and `leftnull(!)` is currently -only implemented for `spacetype(t)<:EuclideanSpace`. +Orthogonality requires `InnerProductStyle(spacetype(t)) <: HasInnerProduct`, and +`leftnull(!)` is currently only implemented for +`InnerProductStyle(spacetype(t)) === EuclideanProduct()`. """ leftnull(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; kwargs...) = leftnull!(permute(t, p1, p2; copy = true); kwargs...) @@ -147,8 +150,9 @@ value decomposition, and one can specify an absolute or relative tolerance for w singular values are to be considered zero, where `max(atol, norm(t)*rtol)` is used as upper bound. -Orthogonality requires `spacetype(t)<:InnerProductSpace`, and `rightnull(!)` is currently -only implemented for `spacetype(t)<:EuclideanSpace`. +Orthogonality requires `InnerProductStyle(spacetype(t)) <: HasInnerProduct`, and +`rightnull(!)` is currently only implemented for +`InnerProductStyle(spacetype(t)) === EuclideanProduct()`. """ rightnull(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; kwargs...) = rightnull!(permute(t, p1, p2; copy = true); kwargs...) @@ -200,13 +204,13 @@ eig(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; kwargs...) = eig!(permute(t, p1, p2; copy = true); kwargs...) """ - eigh(t::AbstractTensorMap{<:EuclideanSpace}, leftind::Tuple, rightind::Tuple) -> D, V + eigh(t::AbstractTensorMap, leftind::Tuple, rightind::Tuple) -> D, V Compute eigenvalue factorization of tensor `t` as linear map from `rightind` to `leftind`. The function `eigh` assumes that the linear map is hermitian and `D` and `V` tensors with the same `eltype` as `t`. See `eig` and `eigen` for non-hermitian tensors. Hermiticity requires that the tensor acts on inner product spaces, and the current implementation -requires `spacetyp(t) <: EuclideanSpace`. +requires `InnerProductStyle(spacetype(t)) === EuclideanProduct()`. If `leftind` and `rightind` are not specified, the current partition of left and right indices of `t` is used. In that case, less memory is allocated if one allows the data in @@ -219,11 +223,12 @@ permute(t, leftind, rightind) * V = V * D See also `eigen` and `eig`. """ -eigh(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple) = - eigh!(permute(t, p1, p2; copy = true)) +function eigh(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple) + return eigh!(permute(t, p1, p2; copy=true)) +end """ - isposdef(t::AbstractTensor{<:EuclideanSpace}, leftind::Tuple, rightind::Tuple) -> ::Bool + isposdef(t::AbstractTensor, leftind::Tuple, rightind::Tuple) -> ::Bool Test whether a tensor `t` is positive definite as linear map from `rightind` to `leftind`. @@ -258,31 +263,50 @@ LinearAlgebra.isposdef(t::AbstractTensorMap) = isposdef!(copy(t)) # Orthogonal factorizations (mutation for recycling memory): # only correct if Euclidean inner product #------------------------------------------------------------------------------------------ -leftorth!(t::AdjointTensorMap{S}; alg::OFA = QRpos()) where {S<:EuclideanSpace} = - map(adjoint, reverse(rightorth!(adjoint(t); alg = alg'))) +function leftorth!(t::AdjointTensorMap; alg::OFA=QRpos()) + InnerProductStyle(spacetype(t)) === EuclideanProduct() || + throw(ArgumentError("leftorth! only defined for Euclidean inner product spaces")) + return map(adjoint, reverse(rightorth!(adjoint(t); alg=alg'))) +end + +function rightorth!(t::AdjointTensorMap; alg::OFA=LQpos()) + InnerProductStyle(spacetype(t)) === EuclideanProduct() || + throw(ArgumentError("rightorth! only defined for Euclidean inner product spaces")) + return map(adjoint, reverse(leftorth!(adjoint(t); alg=alg'))) +end -rightorth!(t::AdjointTensorMap{S}; alg::OFA = LQpos()) where {S<:EuclideanSpace} = - map(adjoint, reverse(leftorth!(adjoint(t); alg = alg'))) +function leftnull!(t::AdjointTensorMap; alg::OFA=QR(), + kwargs...) + InnerProductStyle(spacetype(t)) === EuclideanProduct() || + throw(ArgumentError("leftnull! only defined for Euclidean inner product spaces")) + return adjoint(rightnull!(adjoint(t); alg=alg', kwargs...)) +end -leftnull!(t::AdjointTensorMap{S}; alg::OFA = QR(), kwargs...) where {S<:EuclideanSpace} = - adjoint(rightnull!(adjoint(t); alg = alg', kwargs...)) +function rightnull!(t::AdjointTensorMap; alg::OFA=LQ(), + kwargs...) + InnerProductStyle(spacetype(t)) === EuclideanProduct() || + throw(ArgumentError("rightnull! only defined for Euclidean inner product spaces")) + return adjoint(leftnull!(adjoint(t); alg=alg', kwargs...)) +end -rightnull!(t::AdjointTensorMap{S}; alg::OFA = LQ(), kwargs...) where {S<:EuclideanSpace} = - adjoint(leftnull!(adjoint(t); alg = alg', kwargs...)) +function tsvd!(t::AdjointTensorMap; + trunc::TruncationScheme=NoTruncation(), + p::Real=2, + alg::Union{SVD,SDD}=SDD()) + InnerProductStyle(spacetype(t)) === EuclideanProduct() || + throw(ArgumentError("tsvd! only defined for Euclidean inner product spaces")) -function tsvd!(t::AdjointTensorMap{S}; - trunc::TruncationScheme = NoTruncation(), - p::Real = 2, - alg::Union{SVD, SDD} = SDD()) where {S<:EuclideanSpace} - u, s, vt, err = tsvd!(adjoint(t); trunc = trunc, p = p, alg = alg) + u, s, vt, err = tsvd!(adjoint(t); trunc=trunc, p=p, alg=alg) return adjoint(vt), adjoint(s), adjoint(u), err end -function leftorth!(t::TensorMap{<:EuclideanSpace}; - alg::Union{QR, QRpos, QL, QLpos, SVD, SDD, Polar} = QRpos(), - atol::Real = zero(float(real(eltype(t)))), - rtol::Real = (alg ∉ (SVD(), SDD())) ? zero(float(real(eltype(t)))) : - eps(real(float(one(eltype(t)))))*iszero(atol)) +function leftorth!(t::TensorMap; + alg::Union{QR,QRpos,QL,QLpos,SVD,SDD,Polar}=QRpos(), + atol::Real=zero(float(real(eltype(t)))), + rtol::Real=(alg ∉ (SVD(), SDD())) ? zero(float(real(eltype(t)))) : + eps(real(float(one(eltype(t))))) * iszero(atol)) + InnerProductStyle(spacetype(t)) === EuclideanProduct() || + throw(ArgumentError("leftorth! only defined for Euclidean inner product spaces")) if !iszero(rtol) atol = max(atol, rtol*norm(t)) end @@ -313,11 +337,13 @@ function leftorth!(t::TensorMap{<:EuclideanSpace}; return TensorMap(Qdata, codomain(t)←W), TensorMap(Rdata, W←domain(t)) end -function leftnull!(t::TensorMap{<:EuclideanSpace}; - alg::Union{QR, QRpos, SVD, SDD} = QRpos(), - atol::Real = zero(float(real(eltype(t)))), - rtol::Real = (alg ∉ (SVD(), SDD())) ? zero(float(real(eltype(t)))) : - eps(real(float(one(eltype(t)))))*iszero(atol)) +function leftnull!(t::TensorMap; + alg::Union{QR,QRpos,SVD,SDD}=QRpos(), + atol::Real=zero(float(real(eltype(t)))), + rtol::Real=(alg ∉ (SVD(), SDD())) ? zero(float(real(eltype(t)))) : + eps(real(float(one(eltype(t))))) * iszero(atol)) + InnerProductStyle(spacetype(t)) === EuclideanProduct() || + throw(ArgumentError("leftnull! only defined for Euclidean inner product spaces")) if !iszero(rtol) atol = max(atol, rtol*norm(t)) end @@ -336,11 +362,13 @@ function leftnull!(t::TensorMap{<:EuclideanSpace}; return TensorMap(Ndata, V←W) end -function rightorth!(t::TensorMap{<:EuclideanSpace}; - alg::Union{LQ, LQpos, RQ, RQpos, SVD, SDD, Polar} = LQpos(), - atol::Real = zero(float(real(eltype(t)))), - rtol::Real = (alg ∉ (SVD(), SDD())) ? zero(float(real(eltype(t)))) : - eps(real(float(one(eltype(t)))))*iszero(atol)) +function rightorth!(t::TensorMap; + alg::Union{LQ,LQpos,RQ,RQpos,SVD,SDD,Polar}=LQpos(), + atol::Real=zero(float(real(eltype(t)))), + rtol::Real=(alg ∉ (SVD(), SDD())) ? zero(float(real(eltype(t)))) : + eps(real(float(one(eltype(t))))) * iszero(atol)) + InnerProductStyle(spacetype(t)) === EuclideanProduct() || + throw(ArgumentError("rightorth! only defined for Euclidean inner product spaces")) if !iszero(rtol) atol = max(atol, rtol*norm(t)) end @@ -371,11 +399,13 @@ function rightorth!(t::TensorMap{<:EuclideanSpace}; return TensorMap(Ldata, codomain(t)←W), TensorMap(Qdata, W←domain(t)) end -function rightnull!(t::TensorMap{<:EuclideanSpace}; - alg::Union{LQ, LQpos, SVD, SDD} = LQpos(), - atol::Real = zero(float(real(eltype(t)))), - rtol::Real = (alg ∉ (SVD(), SDD())) ? zero(float(real(eltype(t)))) : - eps(real(float(one(eltype(t)))))*iszero(atol)) +function rightnull!(t::TensorMap; + alg::Union{LQ,LQpos,SVD,SDD}=LQpos(), + atol::Real=zero(float(real(eltype(t)))), + rtol::Real=(alg ∉ (SVD(), SDD())) ? zero(float(real(eltype(t)))) : + eps(real(float(one(eltype(t))))) * iszero(atol)) + InnerProductStyle(spacetype(t)) === EuclideanProduct() || + throw(ArgumentError("rightnull! only defined for Euclidean inner product spaces")) if !iszero(rtol) atol = max(atol, rtol*norm(t)) end @@ -394,10 +424,12 @@ function rightnull!(t::TensorMap{<:EuclideanSpace}; return TensorMap(Ndata, W←V) end -function tsvd!(t::TensorMap{<:EuclideanSpace}; - trunc::TruncationScheme = NoTruncation(), - p::Real = 2, - alg::Union{SVD, SDD} = SDD()) +function tsvd!(t::TensorMap; + trunc::TruncationScheme=NoTruncation(), + p::Real=2, + alg::Union{SVD,SDD}=SDD()) + InnerProductStyle(spacetype(t)) === EuclideanProduct() || + throw(ArgumentError("tsvd! only defined for Euclidean inner product spaces")) S = spacetype(t) I = sectortype(t) A = storagetype(t) @@ -460,7 +492,7 @@ end function LinearAlgebra.ishermitian(t::TensorMap) domain(t) == codomain(t) || return false - spacetype(t) <: EuclideanSpace || return false # hermiticity only defined for euclidean + InnerProductStyle(spacetype(t)) === EuclideanProduct() || return false # hermiticity only defined for euclidean for (c, b) in blocks(t) ishermitian(b) || return false end @@ -469,7 +501,9 @@ end LinearAlgebra.eigen!(t::TensorMap) = ishermitian(t) ? eigh!(t) : eig!(t) -function eigh!(t::TensorMap{<:EuclideanSpace}; kwargs...) +function eigh!(t::TensorMap; kwargs...) + InnerProductStyle(spacetype(t)) === EuclideanProduct() || + throw(ArgumentError("eigh! only defined for Euclidean inner product spaces")) domain(t) == codomain(t) || throw(SpaceMismatch("`eigh!` requires domain and codomain to be the same")) S = spacetype(t) diff --git a/src/tensors/linalg.jl b/src/tensors/linalg.jl index c4c36463..29a227d7 100644 --- a/src/tensors/linalg.jl +++ b/src/tensors/linalg.jl @@ -68,7 +68,7 @@ and error will be thrown. When they are isomorphic, there is no canonical choice specific isomorphism, but the current choice is such that `isomorphism(cod, dom) == inv(isomorphism(dom, cod))`. -See also [`unitary`](@ref) when `spacetype(cod) isa EuclideanSpace`. +See also [`unitary`](@ref) when `InnerProductStyle(spacetype(cod)) === EuclideanProduct()`. """ isomorphism(cod::TensorSpace, dom::TensorSpace) = isomorphism(Matrix{Float64}, cod, dom) isomorphism(P::TensorMapSpace) = isomorphism(codomain(P), domain(P)) @@ -85,50 +85,60 @@ function isomorphism(::Type{A}, cod::ProductSpace, dom::ProductSpace) where {A<: return t end -const EuclideanTensorSpace = TensorSpace{<:EuclideanSpace} -const EuclideanTensorMapSpace = TensorMapSpace{<:EuclideanSpace} -const AbstractEuclideanTensorMap = AbstractTensorMap{<:EuclideanTensorSpace} -const EuclideanTensorMap = TensorMap{<:EuclideanTensorSpace} - """ unitary([A::Type{<:DenseMatrix} = Matrix{Float64},] cod::VectorSpace, dom::VectorSpace) -> TensorMap Return a `t::TensorMap` that implements a specific unitary isomorphism between the codomain -`cod` and the domain `dom`, for which `spacetype(dom)` (`== spacetype(cod)`) must be a -subtype of `EuclideanSpace`. Furthermore, `storagetype(t)` can optionally be chosen to be +`cod` and the domain `dom`, for which `spacetype(dom)` (`== spacetype(cod)`) must have an +inner product. Furthermore, `storagetype(t)` can optionally be chosen to be of type `A`. If the two spaces do not allow for such an isomorphism, and are thus not isomorphic, and error will be thrown. When they are isomorphic, there is no canonical choice for a specific isomorphism, but the current choice is such that `unitary(cod, dom) == inv(unitary(dom, cod)) = adjoint(unitary(dom, cod))`. """ -unitary(cod::EuclideanTensorSpace, dom::EuclideanTensorSpace) = isomorphism(cod, dom) -unitary(P::EuclideanTensorMapSpace) = isomorphism(P) -unitary(A::Type{<:DenseMatrix}, P::EuclideanTensorMapSpace) = isomorphism(A, P) -unitary(A::Type{<:DenseMatrix}, cod::EuclideanTensorSpace, dom::EuclideanTensorSpace) = - isomorphism(A, cod, dom) +function unitary(cod::TensorSpace{S}, dom::TensorSpace{S}) where {S} + InnerProductStyle(S) === EuclideanProduct() || + throw(ArgumentError("unitary requires inner product spaces")) + return isomorphism(cod, dom) +end +function unitary(P::TensorMapSpace{S}) where {S} + InnerProductStyle(S) === EuclideanProduct() || + throw(ArgumentError("unitary requires inner product spaces")) + return isomorphism(P) +end +function unitary(A::Type{<:DenseMatrix}, P::TensorMapSpace{S}) where {S} + InnerProductStyle(S) === EuclideanProduct() || + throw(ArgumentError("unitary requires inner product spaces")) + return isomorphism(A, P) +end +function unitary(A::Type{<:DenseMatrix}, cod::TensorSpace{S}, dom::TensorSpace{S}) where {S} + InnerProductStyle(S) === EuclideanProduct() || + throw(ArgumentError("unitary requires inner product spaces")) + return isomorphism(A, cod, dom) +end """ isometry([A::Type{<:DenseMatrix} = Matrix{Float64},] cod::VectorSpace, dom::VectorSpace) -> TensorMap Return a `t::TensorMap` that implements a specific isometry that embeds the domain `dom` -into the codomain `cod`, and which requires that `spacetype(dom)` (`== spacetype(cod)`) is -a subtype of `EuclideanSpace`. An isometry `t` is such that its adjoint `t'` is the left +into the codomain `cod`, and which requires that `spacetype(dom)` (`== spacetype(cod)`) has +an Euclidean inner product. An isometry `t` is such that its adjoint `t'` is the left inverse of `t`, i.e. `t'*t = id(dom)`, while `t*t'` is some idempotent endomorphism of `cod`, i.e. it squares to itself. When `dom` and `cod` do not allow for such an isometric inclusion, an error will be thrown. """ -isometry(cod::EuclideanTensorSpace, dom::EuclideanTensorSpace) = - isometry(Matrix{Float64}, cod, dom) -isometry(P::EuclideanTensorMapSpace) = isometry(codomain(P), domain(P)) -isometry(A::Type{<:DenseMatrix}, P::EuclideanTensorMapSpace) = - isometry(A, codomain(P), domain(P)) -isometry(A::Type{<:DenseMatrix}, cod::EuclideanTensorSpace, dom::EuclideanTensorSpace) = +isometry(cod::TensorSpace, dom::TensorSpace) = isometry(Matrix{Float64}, cod, dom) +isometry(P::TensorMapSpace) = isometry(codomain(P), domain(P)) +isometry(A::Type{<:DenseMatrix}, P::TensorMapSpace) = isometry(A, codomain(P), domain(P)) +isometry(A::Type{<:DenseMatrix}, cod::TensorSpace, dom::TensorSpace) = isometry(A, convert(ProductSpace, cod), convert(ProductSpace, dom)) function isometry(::Type{A}, cod::ProductSpace{S}, - dom::ProductSpace{S}) where {A<:DenseMatrix, S<:EuclideanSpace} + dom::ProductSpace{S}) where {A<:DenseMatrix, S<:ElementarySpace} + InnerProductStyle(S) === EuclideanProduct() || + throw(ArgumentError("isometries require Euclidean inner product")) dom ≾ cod || throw(SpaceMismatch("codomain $cod and domain $dom do not allow for an isometric mapping")) t = TensorMap(s->A(undef, s), cod, dom) for (c, b) in blocks(t) @@ -159,8 +169,10 @@ function Base.fill!(t::AbstractTensorMap, value::Number) end return t end -function LinearAlgebra.adjoint!(tdst::AbstractEuclideanTensorMap, - tsrc::AbstractEuclideanTensorMap) +function LinearAlgebra.adjoint!(tdst::AbstractTensorMap, + tsrc::AbstractTensorMap) + spacetype(tdst) === spacetype(tsrc) && InnerProductStyle(tdst) === EuclideanProduct() || + throw(ArgumentError("adjoint! requires Euclidean inner product spacetype")) space(tdst) == adjoint(space(tsrc)) || throw(SpaceMismatch()) for c in blocksectors(tdst) adjoint!(StridedView(block(tdst, c)), StridedView(block(tsrc, c))) @@ -203,8 +215,10 @@ function LinearAlgebra.axpby!(α::Number, t1::AbstractTensorMap, end # inner product and norm only valid for spaces with Euclidean inner product -function LinearAlgebra.dot(t1::AbstractEuclideanTensorMap, t2::AbstractEuclideanTensorMap) +function LinearAlgebra.dot(t1::AbstractTensorMap, t2::AbstractTensorMap) space(t1) == space(t2) || throw(SpaceMismatch()) + InnerProductStyle(spacetype(t1)) === EuclideanProduct() || + throw(ArgumentError("dot requires Euclidean inner product")) T = promote_type(eltype(t1), eltype(t2)) s = zero(T) for c in blocksectors(t1) @@ -213,8 +227,11 @@ function LinearAlgebra.dot(t1::AbstractEuclideanTensorMap, t2::AbstractEuclidean return s end -LinearAlgebra.norm(t::AbstractEuclideanTensorMap, p::Real = 2) = - _norm(blocks(t), p, float(zero(real(eltype(t))))) +function LinearAlgebra.norm(t::AbstractTensorMap, p::Real = 2) + InnerProductStyle(spacetype(t)) === EuclideanProduct() || + throw(ArgumentError("norm requires Euclidean inner product")) + return _norm(blocks(t), p, float(zero(real(eltype(t))))) +end function _norm(blockiter, p::Real, init::Real) if p == Inf return mapreduce(max, blockiter; init = init) do (c, b) @@ -477,8 +494,8 @@ function ⊗(t1::AbstractTensorMap{S}, t2::AbstractTensorMap{S}) where S end # deligne product of tensors -function ⊠(t1::AbstractTensorMap{<:EuclideanSpace{ℂ}}, - t2::AbstractTensorMap{<:EuclideanSpace{ℂ}}) +function ⊠(t1::AbstractTensorMap{<:ElementarySpace{ℂ}}, + t2::AbstractTensorMap{<:ElementarySpace{ℂ}}) S1 = spacetype(t1) I1 = sectortype(S1) S2 = spacetype(t2) diff --git a/test/spaces.jl b/test/spaces.jl index c74b973c..22b7e7ba 100644 --- a/test/spaces.jl +++ b/test/spaces.jl @@ -51,8 +51,8 @@ end @test eval(Meta.parse(sprint(show, typeof(V)))) == typeof(V) @test isa(V, VectorSpace) @test isa(V, ElementarySpace) - @test isa(V, InnerProductSpace) - @test isa(V, EuclideanSpace) + @test isa(InnerProductStyle(V), HasInnerProduct) + @test isa(InnerProductStyle(V), EuclideanProduct) @test isa(V, CartesianSpace) @test !isdual(V) @test !isdual(V') @@ -94,8 +94,8 @@ end @test eval(Meta.parse(sprint(show, typeof(V)))) == typeof(V) @test isa(V, VectorSpace) @test isa(V, ElementarySpace) - @test isa(V, InnerProductSpace) - @test isa(V, EuclideanSpace) + @test isa(InnerProductStyle(V), HasInnerProduct) + @test isa(InnerProductStyle(V), EuclideanProduct) @test isa(V, ComplexSpace) @test !isdual(V) @test isdual(V') @@ -151,8 +151,8 @@ end @test TensorKit.isconj(conj(V')) @test isa(V, VectorSpace) @test isa(V, ElementarySpace) - @test !isa(V, InnerProductSpace) - @test !isa(V, EuclideanSpace) + @test !isa(InnerProductStyle(V), HasInnerProduct) + @test !isa(InnerProductStyle(V), EuclideanProduct) @test @constinferred(hash(V)) == hash(deepcopy(V)) != hash(V') @test @constinferred(dual(V)) != @constinferred(conj(V)) != V @test @constinferred(field(V)) == ℂ @@ -206,8 +206,8 @@ end @test eval(Meta.parse(sprint(show, W))) == W @test isa(V, VectorSpace) @test isa(V, ElementarySpace) - @test isa(V, InnerProductSpace) - @test isa(V, EuclideanSpace) + @test isa(InnerProductStyle(V), HasInnerProduct) + @test isa(InnerProductStyle(V), EuclideanProduct) @test isa(V, GradedSpace) @test isa(V, GradedSpace{I}) @test @constinferred(dual(V)) == @constinferred(conj(V)) == @constinferred(adjoint(V)) != V From 35413dd4409afb4ce9cfee0fff134c0ab5ede20a Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sun, 12 Feb 2023 11:59:13 +0100 Subject: [PATCH 2/4] remove unused type vars --- src/tensors/tensoroperations.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/tensors/tensoroperations.jl b/src/tensors/tensoroperations.jl index 26714205..18c1efe1 100644 --- a/src/tensors/tensoroperations.jl +++ b/src/tensors/tensoroperations.jl @@ -18,9 +18,9 @@ function cached_permute(sym::Symbol, t::TensorMap{S}, end end -function cached_permute(sym::Symbol, t::AdjointTensorMap{S}, +function cached_permute(sym::Symbol, t::AdjointTensorMap, p1::IndexTuple, p2::IndexTuple=(); - copy::Bool = false) where {S, N₁, N₂} + copy::Bool = false) p1′ = adjointtensorindices(t, p2) p2′ = adjointtensorindices(t, p1) adjoint(cached_permute(sym, adjoint(t), p1′, p2′; copy = copy)) From edf21459427f5b335ebae0111a89f3cbc1926e92 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 13 Feb 2023 16:32:35 +0100 Subject: [PATCH 3/4] replace stray EuclideanSpace --- src/tensors/factorizations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tensors/factorizations.jl b/src/tensors/factorizations.jl index 5c9e4283..06995ae8 100644 --- a/src/tensors/factorizations.jl +++ b/src/tensors/factorizations.jl @@ -560,7 +560,7 @@ end function LinearAlgebra.isposdef!(t::TensorMap) domain(t) == codomain(t) || throw(SpaceMismatch("`isposdef` requires domain and codomain to be the same")) - spacetype(t) <: EuclideanSpace || return false + InnerProductStyle(spacetype(t)) === EuclideanProduct() || return false for (c, b) in blocks(t) isposdef!(b) || return false end From 61df30fb926a763923f6fb40e3c0ac8cdb59265d Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sat, 29 Apr 2023 21:36:27 +0200 Subject: [PATCH 4/4] update docs --- .gitignore | 1 + docs/src/lib/spaces.md | 2 - docs/src/man/sectors.md | 4 +- docs/src/man/spaces.md | 48 +++++++++++------------ docs/src/man/tensors.md | 74 +++++++++++++++++------------------ docs/src/man/tutorial.md | 4 +- src/tensors/abstracttensor.jl | 54 +++++++++++++++---------- src/tensors/factorizations.jl | 46 +++++++++++----------- src/tensors/linalg.jl | 6 +-- 9 files changed, 123 insertions(+), 116 deletions(-) diff --git a/.gitignore b/.gitignore index e4d7f207..ae7f128c 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,4 @@ *-old __pycache__ .ipynb* +Manifest.toml \ No newline at end of file diff --git a/docs/src/lib/spaces.md b/docs/src/lib/spaces.md index 375394e8..8c2aab18 100644 --- a/docs/src/lib/spaces.md +++ b/docs/src/lib/spaces.md @@ -11,8 +11,6 @@ Field VectorSpace ElementarySpace GeneralSpace -InnerProductSpace -EuclideanSpace CartesianSpace ComplexSpace GradedSpace diff --git a/docs/src/man/sectors.md b/docs/src/man/sectors.md index 8d018f27..ab51b601 100644 --- a/docs/src/man/sectors.md +++ b/docs/src/man/sectors.md @@ -647,9 +647,9 @@ as illustrated below. We have introduced `Sector` subtypes as a way to label the irreps or sectors in the decomposition ``V = ⨁_a ℂ^{n_a} ⊗ R_{a}``. To actually represent such spaces, we now also introduce a corresponding type `GradedSpace`, which is a subtype of -`EuclideanSpace{ℂ}`, i.e. +`ElementarySpace{ℂ}`, i.e. ```julia -struct GradedSpace{I<:Sector, D} <: EuclideanSpace{ℂ} +struct GradedSpace{I<:Sector, D} <: ElementarySpace{ℂ} dims::D dual::Bool end diff --git a/docs/src/man/spaces.md b/docs/src/man/spaces.md index 5eb770f9..75a7428f 100644 --- a/docs/src/man/spaces.md +++ b/docs/src/man/spaces.md @@ -103,40 +103,35 @@ struct GeneralSpace{𝕜} <: ElementarySpace{𝕜} end ``` -We furthermore define the abstract type +We furthermore define the trait types ```julia -abstract type InnerProductSpace{𝕜} <: ElementarySpace{𝕜} end +abstract type InnerProductStyle end +struct NoInnerProduct <: InnerProductStyle end +abstract type HasInnerProduct <: InnerProductStyle end +struct EuclideanProduct <: HasInnerProduct end ``` -to contain all vector spaces `V` which have an inner product and thus a canonical mapping -from `dual(V)` to `V` (for `𝕜 ⊆ ℝ`) or from `dual(V)` to `conj(V)` (otherwise). This -mapping is provided by the metric, but no further support for working with metrics is +to denote for a vector space `V` whether it has an inner product and thus a canonical +mapping from `dual(V)` to `V` (for `𝕜 ⊆ ℝ`) or from `dual(V)` to `conj(V)` (otherwise). +This mapping is provided by the metric, but no further support for working with metrics is currently implemented. -Finally there is -```julia -abstract type EuclideanSpace{𝕜} <: InnerProductSpace{𝕜} end -``` -to contain all spaces `V` with a standard Euclidean inner product (i.e. where the metric is -the identity). These spaces have the natural isomorphisms `dual(V) == V` (for `𝕜 == ℝ`) -or `dual(V) == conj(V)` (for ` 𝕜 == ℂ`). In the language of the previous section on -[categories](@ref s_categories), this subtype represents -[dagger or unitary categories](@ref ss_adjoints), and support an `adjoint` operation. In -particular, we have two concrete types +The `EuclideanProduct` has the natural isomorphisms `dual(V) == V` (for `𝕜 == ℝ`) +or `dual(V) == conj(V)` (for ` 𝕜 == ℂ`). +In the language of the previous section on [categories](@ref s_categories), this trait represents [dagger or unitary categories](@ref ss_adjoints), and these vector spaces support an `adjoint` operation. + +In particular, the two concrete types ```julia -struct CartesianSpace <: EuclideanSpace{ℝ} +struct CartesianSpace <: ElementarySpace{ℝ} d::Int end -struct ComplexSpace <: EuclideanSpace{ℂ} +struct ComplexSpace <: ElementarySpace{ℂ} d::Int dual::Bool end ``` -to represent the Euclidean spaces $ℝ^d$ or $ℂ^d$ without further inner structure. They can -be created using the syntax `CartesianSpace(d) == ℝ^d == ℝ[d]` and -`ComplexSpace(d) == ℂ^d == ℂ[d]`, or -`ComplexSpace(d, true) == ComplexSpace(d; dual = true) == (ℂ^d)' == ℂ[d]'` for the -dual space of the latter. Note that the brackets are required because of the precedence -rules, since `d' == d` for `d::Integer`. +represent the Euclidean spaces $ℝ^d$ or $ℂ^d$ without further inner structure. +They can be created using the syntax `CartesianSpace(d) == ℝ^d == ℝ[d]` and `ComplexSpace(d) == ℂ^d == ℂ[d]`, or `ComplexSpace(d, true) == ComplexSpace(d; dual = true) == (ℂ^d)' == ℂ[d]'` for the dual space of the latter. +Note that the brackets are required because of the precedence rules, since `d' == d` for `d::Integer`. Some examples: ```@repl tensorkit @@ -149,12 +144,15 @@ dual(ℂ^5) == (ℂ^5)' == conj(ℂ^5) == ComplexSpace(5; dual = true) typeof(ℝ^3) spacetype(ℝ^3) spacetype(ℝ[]) +InnerProductStyle(ℝ^3) +InnerProductStyle(ℂ^5) ``` + Note that `ℝ[]` and `ℂ[]` are synonyms for `CartesianSpace` and `ComplexSpace` respectively, such that yet another syntax is e.g. `ℂ[](d)`. This is not very useful in itself, and is motivated by its generalization to `GradedSpace`. We refer to the subsection on [graded spaces](@ref s_rep) on the [next page](@ref s_sectorsrepfusion) for further -information about `GradedSpace`, which is another subtype of `EuclideanSpace{ℂ}` +information about `GradedSpace`, which is another subtype of `ElementarySpace{ℂ}` with an inner structure corresponding to the irreducible representations of a group, or more generally, the simple objects of a fusion category. @@ -279,7 +277,7 @@ For completeness, we also export the strict comparison operators `≺` and `≻` ``` However, as we expect these to be less commonly used, no ASCII alternative is provided. -In the context of `spacetype(V) <: EuclideanSpace`, `V1 ≾ V2` implies that there exists +In the context of `InnerProductStyle(V) <: EuclideanProduct`, `V1 ≾ V2` implies that there exists isometries ``W:V1 → V2`` such that ``W^† ∘ W = \mathrm{id}_{V1}``, while `V1 ≅ V2` implies that there exist unitaries ``U:V1→V2`` such that ``U^† ∘ U = \mathrm{id}_{V1}`` and ``U ∘ U^† = \mathrm{id}_{V2}``. diff --git a/docs/src/man/tensors.md b/docs/src/man/tensors.md index 72b1a3d6..1b4eb452 100644 --- a/docs/src/man/tensors.md +++ b/docs/src/man/tensors.md @@ -391,17 +391,18 @@ can use the method `u = isomorphism([A::Type{<:DenseMatrix}, ] codomain, domain) will explicitly check that the domain and codomain are isomorphic, and return an error otherwise. Again, an optional first argument can be given to specify the specific type of `DenseMatrix` that is currently used to store the rather trivial data of this tensor. If -`spacetype(u) <: EuclideanSpace`, the same result can be obtained with the method `u = -unitary([A::Type{<:DenseMatrix}, ] codomain, domain)`. Note that reversing the domain and -codomain yields the inverse morphism, which in the case of `EuclideanSpace` coincides with -the adjoint morphism, i.e. `isomorphism(A, domain, codomain) == adjoint(u) == inv(u)`, where -`inv` and `adjoint` will be further discussed [below](@ref ss_tensor_linalg). Finally, if -two spaces `V1` and `V2` are such that `V2` can be embedded in `V1`, i.e. there exists an -inclusion with a left inverse, and furthermore they represent tensor products of some -`EuclideanSpace`, the function `w = isometry([A::Type{<:DenseMatrix}, ], V1, V2)` creates -one specific isometric embedding, such that `adjoint(w)*w == id(V2)` and `w*adjoint(w)` is -some hermitian idempotent (a.k.a. orthogonal projector) acting on `V1`. An error will be -thrown if such a map cannot be constructed for the given domain and codomain. +`InnerProductStyle(u) <: EuclideanProduct`, the same result can be obtained with the method +`u = unitary([A::Type{<:DenseMatrix}, ] codomain, domain)`. Note that reversing the domain +and codomain yields the inverse morphism, which in the case of `EuclideanProduct` coincides +with the adjoint morphism, i.e. `isomorphism(A, domain, codomain) == adjoint(u) == inv(u)`, +where `inv` and `adjoint` will be further discussed [below](@ref ss_tensor_linalg). +Finally, if two spaces `V1` and `V2` are such that `V2` can be embedded in `V1`, i.e. there +exists an inclusion with a left inverse, and furthermore they represent tensor products of +some `ElementarySpace` with `EuclideanProduct`, the function +`w = isometry([A::Type{<:DenseMatrix}, ], V1, V2)` creates one specific isometric +embedding, such that `adjoint(w) * w == id(V2)` and `w * adjoint(w)` is some hermitian +idempotent (a.k.a. orthogonal projector) acting on `V1`. An error will be thrown if such a +map cannot be constructed for the given domain and codomain. Let's conclude this section with some examples with `GradedSpace`. ```@repl tensors @@ -575,7 +576,7 @@ can only exist if the domain and codomain are isomorphic, which can e.g. be chec `t2`, we can use the syntax `t1\t2` or `t2/t1`. However, this syntax also accepts instances `t1` whose domain and codomain are not isomorphic, and then amounts to `pinv(t1)`, the Moore-Penrose pseudoinverse. This, however, is only really justified as minimizing the -least squares problem if `spacetype(t) <: EuclideanSpace`. +least squares problem if `InnerProductStyle(t) <: EuclideanProduct`. `AbstractTensorMap` instances behave themselves as vectors (i.e. they are `𝕜`-linear) and so they can be multiplied by scalars and, if they live in the same space, i.e. have the same @@ -588,30 +589,29 @@ operations, TensorKit.jl reexports a number of efficient in-place methods from `lmul!` and `rmul!` (for `y ← α*y` and `y ← y*α`, which is typically the same) and `mul!`, which can also be used for out-of-place scalar multiplication `y ← α*x`. -For `t::AbstractTensorMap{S}` where `S<:EuclideanSpace`, henceforth referred to as -a `(Abstract)EuclideanTensorMap`, we can compute `norm(t)`, and for two such instances, the -inner product `dot(t1, t2)`, provided `t1` and `t2` have the same domain and codomain. -Furthermore, there is `normalize(t)` and `normalize!(t)` to return a scaled version of `t` -with unit norm. These operations should also exist for `S<:InnerProductSpace`, but requires -an interface for defining a custom inner product in these spaces. Currently, there is no -concrete subtype of `InnerProductSpace` that is not a subtype of `EuclideanSpace`. In -particular, `CartesianSpace`, `ComplexSpace` and `GradedSpace` are all subtypes -of `EuclideanSpace`. - -With instances `t::AbstractEuclideanTensorMap` there is associated an adjoint operation, -given by `adjoint(t)` or simply `t'`, such that `domain(t') == codomain(t)` and -`codomain(t') == domain(t)`. Note that for an instance `t::TensorMap{S,N₁,N₂}`, `t'` is -simply stored in a wrapper called `AdjointTensorMap{S,N₂,N₁}`, which is another subtype of -`AbstractTensorMap`. This should be mostly unvisible to the user, as all methods should work -for this type as well. It can be hard to reason about the index order of `t'`, i.e. index -`i` of `t` appears in `t'` at index position `j = TensorKit.adjointtensorindex(t, i)`, -where the latter method is typically not necessary and hence unexported. There is also a -plural `TensorKit.adjointtensorindices` to convert multiple indices at once. Note that, -because the adjoint interchanges domain and codomain, we have -`space(t', j) == space(t, i)'`. +For `t::AbstractTensorMap{S}` where `InnerProductStyle(S) <: EuclideanProduct`, we can +compute `norm(t)`, and for two such instances, the inner product `dot(t1, t2)`, provided +`t1` and `t2` have the same domain and codomain. Furthermore, there is `normalize(t)` and +`normalize!(t)` to return a scaled version of `t` with unit norm. These operations should +also exist for `InnerProductStyle(S) <: HasInnerProduct`, but require an interface for +defining a custom inner product in these spaces. Currently, there is no concrete subtype of +`HasInnerProduct` that is not an `EuclideanProduct`. In particular, `CartesianSpace`, +`ComplexSpace` and `GradedSpace` all have `InnerProductStyle(V) <: EuclideanProduct`. + +With tensors that have `InnerProductStyle(t) <: EuclideanProduct` there is associated an +adjoint operation, given by `adjoint(t)` or simply `t'`, such that +`domain(t') == codomain(t)` and `codomain(t') == domain(t)`. Note that for an instance +`t::TensorMap{S,N₁,N₂}`, `t'` is simply stored in a wrapper called +`AdjointTensorMap{S,N₂,N₁}`, which is another subtype of `AbstractTensorMap`. This should +be mostly unvisible to the user, as all methods should work for this type as well. It can +be hard to reason about the index order of `t'`, i.e. index `i` of `t` appears in `t'` at +index position `j = TensorKit.adjointtensorindex(t, i)`, where the latter method is +typically not necessary and hence unexported. There is also a plural +`TensorKit.adjointtensorindices` to convert multiple indices at once. Note that, because +the adjoint interchanges domain and codomain, we have `space(t', j) == space(t, i)'`. `AbstractTensorMap` instances can furthermore be tested for exact (`t1 == t2`) or -approximate (`t1 ≈ t2`) equality, though the latter requires `norm` can be computed. +approximate (`t1 ≈ t2`) equality, though the latter requires that `norm` can be computed. When tensor map instances are endomorphisms, i.e. they have the same domain and codomain, there is a multiplicative identity which can be obtained as `one(t)` or `one!(t)`, where the @@ -798,8 +798,8 @@ inverse, to split a given index into two or more indices. For a plain tensor (i. data. However, this represents only one possibility, as there is no canonically unique way to embed the tensor product of two spaces `V₁ ⊗ V₂` in a new space `V = fuse(V₁⊗V₂)`. Such a mapping can always be accompagnied by a basis transform. However, one particular choice is -created by the function `isomorphism`, or for `EuclideanSpace` spaces, `unitary`. Hence, we -can join or fuse two indices of a tensor by first constructing +created by the function `isomorphism`, or for `EuclideanProduct` spaces, `unitary`. +Hence, we can join or fuse two indices of a tensor by first constructing `u = unitary(fuse(space(t, i) ⊗ space(t, j)), space(t, i) ⊗ space(t, j))` and then contracting this map with indices `i` and `j` of `t`, as explained in the section on [contracting tensors](@ref ss_tensor_contraction). Note, however, that a typical algorithm @@ -851,7 +851,7 @@ U, Σ, Vʰ, ϵ = tsvd(t; trunc = notrunc(), p::Real = 2, ``` This computes a (possibly truncated) singular value decomposition of -`t::TensorMap{S,N₁,N₂}` (with `S<:EuclideanSpace`), such that +`t::TensorMap{S,N₁,N₂}` (with `InnerProductStyle(t)<:EuclideanProduct`), such that `norm(t - U*Σ*Vʰ) ≈ ϵ`, where `U::TensorMap{S,N₁,1}`, `S::TensorMap{S,1,1}`, `Vʰ::TensorMap{S,1,N₂}` and `ϵ::Real`. `U` is an isometry, i.e. `U'*U` approximates the identity, whereas `U*U'` is an idempotent (squares to itself). The same holds for diff --git a/docs/src/man/tutorial.md b/docs/src/man/tutorial.md index 9a2948bc..56c3f40e 100644 --- a/docs/src/man/tutorial.md +++ b/docs/src/man/tutorial.md @@ -28,12 +28,10 @@ V = ℝ^3 typeof(V) V == CartesianSpace(3) supertype(CartesianSpace) -supertype(EuclideanSpace) -supertype(InnerProductSpace) supertype(ElementarySpace) ``` i.e. `ℝ^n` can also be created without Unicode using the longer syntax `CartesianSpace(n)`. -It is subtype of `EuclideanSpace{ℝ}`, a space with a standard (Euclidean) inner product +It is subtype of `ElementarySpace{ℝ}`, with a standard (Euclidean) inner product over the real numbers. Furthermore, ```@repl tutorial W = ℝ^3 ⊗ ℝ^2 ⊗ ℝ^4 diff --git a/src/tensors/abstracttensor.jl b/src/tensors/abstracttensor.jl index 678e6e0d..1a4abdf1 100644 --- a/src/tensors/abstracttensor.jl +++ b/src/tensors/abstracttensor.jl @@ -10,7 +10,7 @@ of vector spaces of type `S<:IndexSpace`. An `AbstractTensorMap` maps from an input space of type `ProductSpace{S, N₂}` to an output space of type `ProductSpace{S, N₁}`. """ -abstract type AbstractTensorMap{S<:IndexSpace, N₁, N₂} end +abstract type AbstractTensorMap{S<:IndexSpace,N₁,N₂} end """ AbstractTensor{S<:IndexSpace, N} = AbstractTensorMap{S, N, 0} @@ -20,18 +20,20 @@ of type `ProductSpace{S, N}`, built from elementary spaces of type `S<:IndexSpac An `AbstractTensor{S, N}` is actually a special case `AbstractTensorMap{S, N, 0}`, i.e. a tensor map with only a non-trivial output space. """ -const AbstractTensor{S<:IndexSpace, N} = AbstractTensorMap{S, N, 0} +const AbstractTensor{S<:IndexSpace,N} = AbstractTensorMap{S,N,0} # tensor characteristics Base.eltype(T::Type{<:AbstractTensorMap}) = eltype(storagetype(T)) -similarstoragetype(TT::Type{<:AbstractTensorMap}, ::Type{T}) where {T} = - Core.Compiler.return_type(similar, Tuple{storagetype(TT), Type{T}}) +function similarstoragetype(TT::Type{<:AbstractTensorMap}, ::Type{T}) where {T} + return Core.Compiler.return_type(similar, Tuple{storagetype(TT),Type{T}}) +end storagetype(t::AbstractTensorMap) = storagetype(typeof(t)) similarstoragetype(t::AbstractTensorMap, T) = similarstoragetype(typeof(t), T) Base.eltype(t::AbstractTensorMap) = eltype(typeof(t)) spacetype(t::AbstractTensorMap) = spacetype(typeof(t)) sectortype(t::AbstractTensorMap) = sectortype(typeof(t)) +InnerProductStyle(t::AbstractTensorMap) = InnerProductStyle(typeof(t)) field(t::AbstractTensorMap) = field(typeof(t)) numout(t::AbstractTensorMap) = numout(typeof(t)) numin(t::AbstractTensorMap) = numin(typeof(t)) @@ -39,10 +41,13 @@ numind(t::AbstractTensorMap) = numind(typeof(t)) spacetype(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace} = S sectortype(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace} = sectortype(S) +function InnerProductStyle(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace} + return InnerProductStyle(S) +end field(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace} = field(S) -numout(::Type{<:AbstractTensorMap{<:IndexSpace, N₁, N₂}}) where {N₁, N₂} = N₁ -numin(::Type{<:AbstractTensorMap{<:IndexSpace, N₁, N₂}}) where {N₁, N₂} = N₂ -numind(::Type{<:AbstractTensorMap{<:IndexSpace, N₁, N₂}}) where {N₁, N₂} = N₁ + N₂ +numout(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁,N₂} = N₁ +numin(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁,N₂} = N₂ +numind(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁,N₂} = N₁ + N₂ const order = numind @@ -56,22 +61,27 @@ space(t::AbstractTensorMap, i::Int) = space(t)[i] dim(t::AbstractTensorMap) = dim(space(t)) # some index manipulation utilities -codomainind(::Type{<:AbstractTensorMap{<:IndexSpace, N₁, N₂}}) where {N₁, N₂} = - ntuple(n->n, N₁) -domainind(::Type{<:AbstractTensorMap{<:IndexSpace, N₁, N₂}}) where {N₁, N₂} = - ntuple(n-> N₁+n, N₂) -allind(::Type{<:AbstractTensorMap{<:IndexSpace, N₁, N₂}}) where {N₁, N₂} = - ntuple(n->n, N₁+N₂) +function codomainind(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁,N₂} + return ntuple(n -> n, N₁) +end +function domainind(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁,N₂} + return ntuple(n -> N₁ + n, N₂) +end +function allind(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁,N₂} + return ntuple(n -> n, N₁ + N₂) +end codomainind(t::AbstractTensorMap) = codomainind(typeof(t)) domainind(t::AbstractTensorMap) = domainind(typeof(t)) allind(t::AbstractTensorMap) = allind(typeof(t)) -adjointtensorindex(t::AbstractTensorMap{<:IndexSpace, N₁, N₂}, i) where {N₁, N₂} = - ifelse(i<=N₁, N₂+i, i-N₁) +function adjointtensorindex(t::AbstractTensorMap{<:IndexSpace,N₁,N₂}, i) where {N₁,N₂} + return ifelse(i <= N₁, N₂ + i, i - N₁) +end -adjointtensorindices(t::AbstractTensorMap, indices::IndexTuple) = - map(i->adjointtensorindex(t, i), indices) +function adjointtensorindices(t::AbstractTensorMap, indices::IndexTuple) + return map(i -> adjointtensorindex(t, i), indices) +end # Equality and approximality #---------------------------- @@ -92,10 +102,11 @@ function Base.hash(t::AbstractTensorMap, h::UInt) end function Base.isapprox(t1::AbstractTensorMap, t2::AbstractTensorMap; - atol::Real=0, rtol::Real=Base.rtoldefault(eltype(t1), eltype(t2), atol)) + atol::Real=0, + rtol::Real=Base.rtoldefault(eltype(t1), eltype(t2), atol)) d = norm(t1 - t2) if isfinite(d) - return d <= max(atol, rtol*max(norm(t1), norm(t2))) + return d <= max(atol, rtol * max(norm(t1), norm(t2))) else return false end @@ -104,7 +115,7 @@ end # Conversion to Array: #---------------------- # probably not optimized for speed, only for checking purposes -function Base.convert(::Type{Array}, t::AbstractTensorMap{S, N₁, N₂}) where {S, N₁, N₂} +function Base.convert(::Type{Array}, t::AbstractTensorMap{S,N₁,N₂}) where {S,N₁,N₂} I = sectortype(t) if I === Trivial convert(Array, t[]) @@ -119,7 +130,8 @@ function Base.convert(::Type{Array}, t::AbstractTensorMap{S, N₁, N₂}) where sz2 = size(F2) d1 = TupleTools.front(sz1) d2 = TupleTools.front(sz2) - F = reshape(reshape(F1, TupleTools.prod(d1), sz1[end])*reshape(F2, TupleTools.prod(d2), sz2[end])', (d1..., d2...)) + F = reshape(reshape(F1, TupleTools.prod(d1), sz1[end]) * + reshape(F2, TupleTools.prod(d2), sz2[end])', (d1..., d2...)) if !(@isdefined A) if eltype(F) <: Complex T = complex(float(eltype(t))) diff --git a/src/tensors/factorizations.jl b/src/tensors/factorizations.jl index 06995ae8..6c0bc65a 100644 --- a/src/tensors/factorizations.jl +++ b/src/tensors/factorizations.jl @@ -47,8 +47,8 @@ singular values that were truncated. THe keyword `alg` can be equal to `SVD()` or `SDD()`, corresponding to the underlying LAPACK algorithm that computes the decomposition (`_gesvd` or `_gesdd`). -Orthogonality requires `InnerProductStyle(spacetype(t)) <: HasInnerProduct`, and `tsvd(!)` -is currently only implemented for `InnerProductStyle(spacetype(t)) === EuclideanProduct()`. +Orthogonality requires `InnerProductStyle(t) <: HasInnerProduct`, and `tsvd(!)` +is currently only implemented for `InnerProductStyle(t) === EuclideanProduct()`. """ tsvd(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; kwargs...) = tsvd!(permute(t, p1, p2; copy = true); kwargs...) @@ -71,9 +71,9 @@ standard QR decomposition such that the diagonal elements of `R` are positive. O `QRpos()` and `Polar()` are uniqe (no residual freedom) so that they always return the same result for the same input tensor `t`. -Orthogonality requires `InnerProductStyle(spacetype(t)) <: HasInnerProduct`, and +Orthogonality requires `InnerProductStyle(t) <: HasInnerProduct`, and `leftorth(!)` is currently only implemented for - `InnerProductStyle(spacetype(t)) === EuclideanProduct()`. + `InnerProductStyle(t) === EuclideanProduct()`. """ leftorth(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; kwargs...) = leftorth!(permute(t, p1, p2; copy = true); kwargs...) @@ -98,9 +98,9 @@ diagonal elements of `L` are positive. `Polar()` produces a Hermitian and positi semidefinite `L`. Only `LQpos()`, `RQpos()` and `Polar()` are uniqe (no residual freedom) so that they always return the same result for the same input tensor `t`. -Orthogonality requires `InnerProductStyle(spacetype(t)) <: HasInnerProduct`, and +Orthogonality requires `InnerProductStyle(t) <: HasInnerProduct`, and `rightorth(!)` is currently only implemented for -`InnerProductStyle(spacetype(t)) === EuclideanProduct()`. +`InnerProductStyle(t) === EuclideanProduct()`. """ rightorth(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; kwargs...) = rightorth!(permute(t, p1, p2; copy = true); kwargs...) @@ -123,9 +123,9 @@ value decomposition, and one can specify an absolute or relative tolerance for w singular values are to be considered zero, where `max(atol, norm(t)*rtol)` is used as upper bound. -Orthogonality requires `InnerProductStyle(spacetype(t)) <: HasInnerProduct`, and +Orthogonality requires `InnerProductStyle(t) <: HasInnerProduct`, and `leftnull(!)` is currently only implemented for -`InnerProductStyle(spacetype(t)) === EuclideanProduct()`. +`InnerProductStyle(t) === EuclideanProduct()`. """ leftnull(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; kwargs...) = leftnull!(permute(t, p1, p2; copy = true); kwargs...) @@ -150,9 +150,9 @@ value decomposition, and one can specify an absolute or relative tolerance for w singular values are to be considered zero, where `max(atol, norm(t)*rtol)` is used as upper bound. -Orthogonality requires `InnerProductStyle(spacetype(t)) <: HasInnerProduct`, and +Orthogonality requires `InnerProductStyle(t) <: HasInnerProduct`, and `rightnull(!)` is currently only implemented for -`InnerProductStyle(spacetype(t)) === EuclideanProduct()`. +`InnerProductStyle(t) === EuclideanProduct()`. """ rightnull(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; kwargs...) = rightnull!(permute(t, p1, p2; copy = true); kwargs...) @@ -210,7 +210,7 @@ Compute eigenvalue factorization of tensor `t` as linear map from `rightind` to The function `eigh` assumes that the linear map is hermitian and `D` and `V` tensors with the same `eltype` as `t`. See `eig` and `eigen` for non-hermitian tensors. Hermiticity requires that the tensor acts on inner product spaces, and the current implementation -requires `InnerProductStyle(spacetype(t)) === EuclideanProduct()`. +requires `InnerProductStyle(t) === EuclideanProduct()`. If `leftind` and `rightind` are not specified, the current partition of left and right indices of `t` is used. In that case, less memory is allocated if one allows the data in @@ -264,27 +264,27 @@ LinearAlgebra.isposdef(t::AbstractTensorMap) = isposdef!(copy(t)) # only correct if Euclidean inner product #------------------------------------------------------------------------------------------ function leftorth!(t::AdjointTensorMap; alg::OFA=QRpos()) - InnerProductStyle(spacetype(t)) === EuclideanProduct() || + InnerProductStyle(t) === EuclideanProduct() || throw(ArgumentError("leftorth! only defined for Euclidean inner product spaces")) return map(adjoint, reverse(rightorth!(adjoint(t); alg=alg'))) end function rightorth!(t::AdjointTensorMap; alg::OFA=LQpos()) - InnerProductStyle(spacetype(t)) === EuclideanProduct() || + InnerProductStyle(t) === EuclideanProduct() || throw(ArgumentError("rightorth! only defined for Euclidean inner product spaces")) return map(adjoint, reverse(leftorth!(adjoint(t); alg=alg'))) end function leftnull!(t::AdjointTensorMap; alg::OFA=QR(), kwargs...) - InnerProductStyle(spacetype(t)) === EuclideanProduct() || + InnerProductStyle(t) === EuclideanProduct() || throw(ArgumentError("leftnull! only defined for Euclidean inner product spaces")) return adjoint(rightnull!(adjoint(t); alg=alg', kwargs...)) end function rightnull!(t::AdjointTensorMap; alg::OFA=LQ(), kwargs...) - InnerProductStyle(spacetype(t)) === EuclideanProduct() || + InnerProductStyle(t) === EuclideanProduct() || throw(ArgumentError("rightnull! only defined for Euclidean inner product spaces")) return adjoint(leftnull!(adjoint(t); alg=alg', kwargs...)) end @@ -293,7 +293,7 @@ function tsvd!(t::AdjointTensorMap; trunc::TruncationScheme=NoTruncation(), p::Real=2, alg::Union{SVD,SDD}=SDD()) - InnerProductStyle(spacetype(t)) === EuclideanProduct() || + InnerProductStyle(t) === EuclideanProduct() || throw(ArgumentError("tsvd! only defined for Euclidean inner product spaces")) u, s, vt, err = tsvd!(adjoint(t); trunc=trunc, p=p, alg=alg) @@ -305,7 +305,7 @@ function leftorth!(t::TensorMap; atol::Real=zero(float(real(eltype(t)))), rtol::Real=(alg ∉ (SVD(), SDD())) ? zero(float(real(eltype(t)))) : eps(real(float(one(eltype(t))))) * iszero(atol)) - InnerProductStyle(spacetype(t)) === EuclideanProduct() || + InnerProductStyle(t) === EuclideanProduct() || throw(ArgumentError("leftorth! only defined for Euclidean inner product spaces")) if !iszero(rtol) atol = max(atol, rtol*norm(t)) @@ -342,7 +342,7 @@ function leftnull!(t::TensorMap; atol::Real=zero(float(real(eltype(t)))), rtol::Real=(alg ∉ (SVD(), SDD())) ? zero(float(real(eltype(t)))) : eps(real(float(one(eltype(t))))) * iszero(atol)) - InnerProductStyle(spacetype(t)) === EuclideanProduct() || + InnerProductStyle(t) === EuclideanProduct() || throw(ArgumentError("leftnull! only defined for Euclidean inner product spaces")) if !iszero(rtol) atol = max(atol, rtol*norm(t)) @@ -367,7 +367,7 @@ function rightorth!(t::TensorMap; atol::Real=zero(float(real(eltype(t)))), rtol::Real=(alg ∉ (SVD(), SDD())) ? zero(float(real(eltype(t)))) : eps(real(float(one(eltype(t))))) * iszero(atol)) - InnerProductStyle(spacetype(t)) === EuclideanProduct() || + InnerProductStyle(t) === EuclideanProduct() || throw(ArgumentError("rightorth! only defined for Euclidean inner product spaces")) if !iszero(rtol) atol = max(atol, rtol*norm(t)) @@ -404,7 +404,7 @@ function rightnull!(t::TensorMap; atol::Real=zero(float(real(eltype(t)))), rtol::Real=(alg ∉ (SVD(), SDD())) ? zero(float(real(eltype(t)))) : eps(real(float(one(eltype(t))))) * iszero(atol)) - InnerProductStyle(spacetype(t)) === EuclideanProduct() || + InnerProductStyle(t) === EuclideanProduct() || throw(ArgumentError("rightnull! only defined for Euclidean inner product spaces")) if !iszero(rtol) atol = max(atol, rtol*norm(t)) @@ -428,7 +428,7 @@ function tsvd!(t::TensorMap; trunc::TruncationScheme=NoTruncation(), p::Real=2, alg::Union{SVD,SDD}=SDD()) - InnerProductStyle(spacetype(t)) === EuclideanProduct() || + InnerProductStyle(t) === EuclideanProduct() || throw(ArgumentError("tsvd! only defined for Euclidean inner product spaces")) S = spacetype(t) I = sectortype(t) @@ -492,7 +492,7 @@ end function LinearAlgebra.ishermitian(t::TensorMap) domain(t) == codomain(t) || return false - InnerProductStyle(spacetype(t)) === EuclideanProduct() || return false # hermiticity only defined for euclidean + InnerProductStyle(t) === EuclideanProduct() || return false # hermiticity only defined for euclidean for (c, b) in blocks(t) ishermitian(b) || return false end @@ -502,7 +502,7 @@ end LinearAlgebra.eigen!(t::TensorMap) = ishermitian(t) ? eigh!(t) : eig!(t) function eigh!(t::TensorMap; kwargs...) - InnerProductStyle(spacetype(t)) === EuclideanProduct() || + InnerProductStyle(t) === EuclideanProduct() || throw(ArgumentError("eigh! only defined for Euclidean inner product spaces")) domain(t) == codomain(t) || throw(SpaceMismatch("`eigh!` requires domain and codomain to be the same")) diff --git a/src/tensors/linalg.jl b/src/tensors/linalg.jl index 29a227d7..d5cab859 100644 --- a/src/tensors/linalg.jl +++ b/src/tensors/linalg.jl @@ -68,7 +68,7 @@ and error will be thrown. When they are isomorphic, there is no canonical choice specific isomorphism, but the current choice is such that `isomorphism(cod, dom) == inv(isomorphism(dom, cod))`. -See also [`unitary`](@ref) when `InnerProductStyle(spacetype(cod)) === EuclideanProduct()`. +See also [`unitary`](@ref) when `InnerProductStyle(cod) === EuclideanProduct()`. """ isomorphism(cod::TensorSpace, dom::TensorSpace) = isomorphism(Matrix{Float64}, cod, dom) isomorphism(P::TensorMapSpace) = isomorphism(codomain(P), domain(P)) @@ -217,7 +217,7 @@ end # inner product and norm only valid for spaces with Euclidean inner product function LinearAlgebra.dot(t1::AbstractTensorMap, t2::AbstractTensorMap) space(t1) == space(t2) || throw(SpaceMismatch()) - InnerProductStyle(spacetype(t1)) === EuclideanProduct() || + InnerProductStyle(t1) === EuclideanProduct() || throw(ArgumentError("dot requires Euclidean inner product")) T = promote_type(eltype(t1), eltype(t2)) s = zero(T) @@ -228,7 +228,7 @@ function LinearAlgebra.dot(t1::AbstractTensorMap, t2::AbstractTensorMap) end function LinearAlgebra.norm(t::AbstractTensorMap, p::Real = 2) - InnerProductStyle(spacetype(t)) === EuclideanProduct() || + InnerProductStyle(t) === EuclideanProduct() || throw(ArgumentError("norm requires Euclidean inner product")) return _norm(blocks(t), p, float(zero(real(eltype(t))))) end