Skip to content

Commit

Permalink
Merge pull request #72 from lkdvos/spacetypes
Browse files Browse the repository at this point in the history
Refactor: abstract spacetypes to traits

- Removes `InnerProductSpace` and `EuclideanSpace` as abstract types.
- Adds `InnerProductStyle` as a trait based mechanism to query if a space has an inner product.
- Implements `NoInnerProduct`, `HasInnerProduct` and `EuclideanInnerProduct` as traits
  • Loading branch information
lkdvos authored May 17, 2023
2 parents c5aa8cb + 61df30f commit 820f092
Show file tree
Hide file tree
Showing 19 changed files with 328 additions and 249 deletions.
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -5,3 +5,4 @@
*-old
__pycache__
.ipynb*
Manifest.toml
2 changes: 0 additions & 2 deletions docs/src/lib/spaces.md
Original file line number Diff line number Diff line change
Expand Up @@ -11,8 +11,6 @@ Field
VectorSpace
ElementarySpace
GeneralSpace
InnerProductSpace
EuclideanSpace
CartesianSpace
ComplexSpace
GradedSpace
Expand Down
4 changes: 2 additions & 2 deletions docs/src/man/sectors.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
48 changes: 23 additions & 25 deletions docs/src/man/spaces.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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.

Expand Down Expand Up @@ -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}``.
Expand Down
74 changes: 37 additions & 37 deletions docs/src/man/tensors.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
4 changes: 1 addition & 3 deletions docs/src/man/tutorial.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
3 changes: 2 additions & 1 deletion src/TensorKit.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
8 changes: 6 additions & 2 deletions src/spaces/cartesianspace.jl
Original file line number Diff line number Diff line change
@@ -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))
Expand All @@ -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)
Expand Down
36 changes: 22 additions & 14 deletions src/spaces/complexspace.jl
Original file line number Diff line number Diff line change
@@ -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)
Expand All @@ -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)
Expand All @@ -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)")
Loading

0 comments on commit 820f092

Please sign in to comment.