Return the eigendecomposition spectrum as an unmatricized array

Project.toml

@@ -1,6 +1,6 @@
 name = "TensorAlgebra"
 uuid = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
-version = "0.19.3"
+version = "0.20.0"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
 
 [workspace]

docs/Project.toml

@@ -11,4 +11,4 @@ path = ".."
 Documenter = "1.8.1"
 ITensorFormatter = "0.2.27"
 Literate = "2.20.1"
-TensorAlgebra = "0.19"
+TensorAlgebra = "0.20"

examples/Project.toml

@@ -5,4 +5,4 @@ TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
 path = ".."
 
 [compat]
-TensorAlgebra = "0.19"
+TensorAlgebra = "0.20"

src/factorizations.jl

@@ -286,15 +286,16 @@ function svd_trunc(A, ndims_codomain::Val; kwargs...)
     return svd_trunc(MatricizeStyle(A), A, ndims_codomain; kwargs...)
 end
 
-# Eigendecomposition: `D` is the rank × rank spectrum, left as a matrix, while `V` carries
-# the codomain axes plus a trailing rank axis.
+# Eigendecomposition: `D` is the rank × rank spectrum and `V` carries the codomain axes plus a
+# trailing rank axis. Both are unmatricized back to the array type, as in `svd_*`.
 for f in (:eigh_full, :eig_full, :eigh_trunc, :eig_trunc)
     @eval begin
         function $f(style::MatricizeStyle, A, ndims_codomain::Val; kwargs...)
             A_mat = matricize(style, A, ndims_codomain)
             D, V = MatrixAlgebraKit.$f(A_mat; kwargs...)
             axes_codomain = first(bipartition(axes(A), ndims_codomain))
-            return D, unmatricize(style, V, axes_codomain, (conj(axes(V, ndims(V))),))
+            return unmatricize(style, D, (axes(D, 1),), (conj(axes(D, 2)),)),
+                unmatricize(style, V, axes_codomain, (conj(axes(V, ndims(V))),))
         end
         function $f(A, ndims_codomain::Val; kwargs...)
             return $f(MatricizeStyle(A), A, ndims_codomain; kwargs...)

test/Project.toml

@@ -37,7 +37,7 @@ Random = "1.10"
 SafeTestsets = "0.1"
 StableRNGs = "1.0.2"
 Suppressor = "0.2"
-TensorAlgebra = "0.19"
+TensorAlgebra = "0.20"
 TensorKit = "0.17"
 TensorOperations = "5.1.4"
 Test = "1.10"

test/test_factorizations.jl

@@ -1,4 +1,4 @@
-using LinearAlgebra: LinearAlgebra, I, diag, norm
+using LinearAlgebra: LinearAlgebra, Diagonal, I, diag, norm
 using MatrixAlgebraKit: truncrank
 using TensorAlgebra: TensorAlgebra, contract, eig_full, eig_vals, eigh_full, eigh_vals,
     gram_eigh_full, gram_eigh_full_with_pinv, left_null, left_orth, left_polar, lq_compact,
@@ -90,6 +90,8 @@ end
     D, V = eig_full(A, labels_A, labels_V, labels_V′)
     @test A == Acopy # should not have altered initial array
     @test eltype(D) == eltype(V) && eltype(D) <: Complex
+    # `D` is unmatricized like the other factors; the `(1, 1)` unmatricize preserves `Diagonal`.
+    @test D isa Diagonal
 
     AV = contract((:a, :b, :D), A, labels_A, V, (labels_V′..., :D))
     VD = contract((:a, :b, :D), V, (labels_V..., :D′), D, (:D′, :D))
@@ -112,6 +114,7 @@ end
     @test A == Acopy # should not have altered initial array
     @test eltype(D) <: Real
     @test eltype(V) == eltype(A)
+    @test D isa Diagonal
 
     AV = contract((:a, :b, :D), A, labels_A, V, (labels_V′..., :D))
     VD = contract((:a, :b, :D), V, (labels_V..., :D′), D, (:D′, :D))
@@ -167,6 +170,7 @@ end
     Acopy = copy(A)
     U, S, Vᴴ = @constinferred svd_compact(A, labels_A, labels_U, labels_Vᴴ)
     @test A == Acopy # should not have altered initial array
+    @test S isa Diagonal
     US, labels_US = contract(U, (labels_U..., :u), S, (:u, :v))
     A′ = contract(labels_A, US, labels_US, Vᴴ, (:v, labels_Vᴴ...))
     @test A ≈ A′