Implement Eigendecomposition

src/Numerics.jl

@@ -119,6 +119,108 @@ function factorinds(tensor, left_inds, right_inds)
  return left_inds, right_inds
 end
 
+# TODO is this an `AbstractTensorNetwork`?
+# TODO add fancier `show` method
+struct TensorEigen{T,V,Nᵣ,S<:AbstractVector{V},U<:AbstractArray{T,Nᵣ}} <: Factorization{T}
+ values::Tensor{V,1,S}
+ vectors::Tensor{T,Nᵣ,U}
+ right_inds::Vector{Symbol}
+end
+
+function Base.getproperty(obj::TensorEigen, name::Symbol)
+ if name === :U
+ return obj.vectors
+ elseif name === :Λ
+ return obj.values
+ elseif name ∈ [:Uinv, :U⁻¹]
+ U = reshape(parent(obj.vectors), prod(size(obj.vectors)[1:(end - 1)]), size(obj.vectors)[end])
+ Uinv = inv(U)
+ return Tensor(Uinv, [only(inds(obj.values)), obj.right_inds...])
+ end
+ return getfield(obj, name)
+end
+
+function Base.inv(F::TensorEigen)
+ U = reshape(parent(F.vectors), prod(size(F.vectors)[1:(end - 1)]), size(F.vectors)[end])
+ left_inds = inds(F.vectors)[1:(end - 1)]
+ return Tensor(U * inv(Diagonal(F.values)) / U, [left_inds..., F.right_inds...])
+end
+LinearAlgebra.det(x::TensorEigen) = prod(x.values)
+
+Base.iterate(x::TensorEigen) = (x.values, :vectors)
+Base.iterate(x::TensorEigen, state) = state == :vectors ? (x.vectors, nothing) : nothing
+
+LinearAlgebra.eigen(t::Tensor{<:Any,2}; kwargs...) = @invoke eigen(t::Tensor; left_inds=(first(inds(t)),), kwargs...)
+
+"""
+ LinearAlgebra.eigen(tensor::Tensor; left_inds, right_inds, virtualind, kwargs...)
+
+Perform Eigendecomposition on a tensor.
+
+# Keyword arguments
+
+ - `left_inds`: left indices to be used in the eigendecomposition. Defaults to all indices of `t` except `right_inds`.
+ - `right_inds`: right indices to be used in the eigendecomposition. Defaults to all indices of `t` except `left_inds`.
+ - `virtualind`: name of the virtual bond. Defaults to a random `Symbol`.
+"""
+function LinearAlgebra.eigen(tensor::Tensor; left_inds=(), right_inds=(), virtualind=Symbol(uuid4()), kwargs...)
+ left_inds, right_inds = factorinds(tensor, left_inds, right_inds)
+
+ virtualind ∉ inds(tensor) ||
+ throw(ArgumentError("new virtual bond name ($virtualind) cannot be already be present"))
+
+ # permute array
+ left_sizes = map(Base.Fix1(size, tensor), left_inds)
+ right_sizes = map(Base.Fix1(size, tensor), right_inds)
+ tensor = permutedims(tensor, [left_inds..., right_inds...])
+ data = reshape(parent(tensor), prod(left_sizes), prod(right_sizes))
+
+ # compute eigendecomposition
+ Λ, U = eigen(data; kwargs...)
+
+ # tensorify results
+ Λ = Tensor(Λ, [virtualind])
+ U = Tensor(reshape(U, left_sizes..., size(U, 2)), [left_inds..., virtualind])
+
+ return TensorEigen(Λ, U, right_inds)
+end
+
+LinearAlgebra.eigvals(t::Tensor{<:Any,2}; kwargs...) = eigvals(parent(t); kwargs...)
+
+# TODO document when it returns a `Tensor` and when returns an `Array`
+"""
+ LinearAlgebra.eigvals(tensor::Tensor; left_inds, right_inds, kwargs...)
+
+Perform Eigendecomposition on a tensor and return the eigenvalues.
+
+# Keyword arguments
+
+ - `left_inds`: left indices to be used in the eigendecomposition. Defaults to all indices of `t` except `right_inds`.
+ - `right_inds`: right indices to be used in the eigendecomposition. Defaults to all indices of `t` except `left_inds`.
+"""
+function LinearAlgebra.eigvals(tensor::Tensor; left_inds=(), right_inds=(), kwargs...)
+ F = eigen(tensor; left_inds, right_inds, kwargs...)
+ return parent(F.values)
+end
+
+LinearAlgebra.eigvecs(t::Tensor{<:Any,2}; kwargs...) = eigvecs(parent(t); kwargs...)
+
+"""
+ LinearAlgebra.eigvecs(tensor::Tensor; left_inds, right_inds, kwargs...)
+
+Perform Eigendecomposition on a tensor and return the eigenvectors.
+
+# Keyword arguments
+
+ - `left_inds`: left indices to be used in the eigendecomposition. Defaults to all indices of `t` except `right_inds`.
+ - `right_inds`: right indices to be used in the eigendecomposition. Defaults to all indices of `t` except `left_inds`.
+ - `virtualind`: name of the virtual bond. Defaults to a random `Symbol`.
+"""
+function LinearAlgebra.eigvecs(tensor::Tensor; left_inds=(), right_inds=(), kwargs...)
+ F = eigen(tensor; left_inds, right_inds, kwargs...)
+ return F.vectors
+end
+
 LinearAlgebra.svd(t::Tensor{<:Any,2}; kwargs...) = Base.@invoke svd(t::Tensor; left_inds=(first(inds(t)),), kwargs...)
 
 """

src/Tensor.jl

@@ -243,3 +243,7 @@ function __expand_repeat(array, axis, size)
 end
 
 LinearAlgebra.opnorm(x::Tensor, p::Real) = opnorm(parent(x), p)
+
+LinearAlgebra.det(x::Tensor{T,2}) where {T} = det(parent(x))
+LinearAlgebra.logdet(x::Tensor{T,2}) where {T} = logdet(parent(x))
+LinearAlgebra.tr(x::Tensor{T,2}) where {T} = tr(parent(x))

src/TensorNetwork.jl

@@ -641,6 +641,13 @@ contract(t::Tensor, tn::AbstractTensorNetwork; kwargs...) = contract(tn, t; kwar
  return contract(intermediates...; dims=suminds(path))
 end
 
+function LinearAlgebra.eigen!(tn::AbstractTensorNetwork; left_inds=Symbol[], right_inds=Symbol[], kwargs...)
+ tensor = tn[left_inds ∪ right_inds...]
+ (; U, Λ, U⁻¹) = eigen(tensor; left_inds, right_inds, kwargs...)
+ replace!(tn, tensor => TensorNetwork([U, Λ, U⁻¹]))
+ return tn
+end
+
 function LinearAlgebra.svd!(tn::AbstractTensorNetwork; left_inds=Symbol[], right_inds=Symbol[], kwargs...)
  tensor = tn[left_inds ∪ right_inds...]
  U, s, Vt = svd(tensor; left_inds, right_inds, kwargs...)

test/Numerics_test.jl

@@ -138,6 +138,43 @@
  end
  end
 
+ @testset "eigen" begin
+ A = Tensor(rand(2, 2, 2, 2), (:i, :j, :k, :l))
+
+ # throw if index is not present
+ @test_throws ArgumentError eigen(A, left_inds=[:z])
+ @test_throws ArgumentError eigen(A, right_inds=[:z])
+
+ # throw if no inds left
+ @test_throws ArgumentError eigen(A, left_inds=(:i, :j, :k, :l))
+ @test_throws ArgumentError eigen(A, right_inds=(:i, :j, :k, :l))
+
+ # throw if chosen virtual index already present
+ @test_throws ArgumentError eigen(A, left_inds=(:i,), virtualind=:j)
+
+ F = eigen(A; left_inds=[:i, :j], virtualind=:x)
+
+ # iteration deconstruction support
+ @test try
+ Λ, U = F
+ true
+ catch
+ false
+ end
+
+ (; Λ, U, U⁻¹) = F
+
+ @test inds(U) == [:i, :j, :x]
+ @test inds(Λ) == [:x]
+ @test inds(U⁻¹) == [:k, :l, :x]
+
+ @test size(U) == (2, 2, 2)
+ @test size(Λ) == (2,)
+ @test size(U⁻¹) == (2, 2, 2)
+
+ @test isapprox(contract(contract(U, Λ; out=inds(U)), U⁻¹), A)
+ end
+
  @testset "svd" begin
  data = rand(ComplexF64, 2, 4, 6, 8)
  tensor = Tensor(data, (:i, :j, :k, :l))