Tensor factorizations

As tensors are linear maps, they support various kinds of factorizations. These functions all interpret the provided AbstractTensorMap instances as a map from domain to codomain, which can be thought of as reshaping the tensor into a matrix according to the current bipartition of the indices.

TensorKit's factorizations are provided by MatrixAlgebraKit.jl, which is used to supply both the interface, as well as the implementation of the various operations on the blocks of data. For specific details on the provided functionality, we refer to its documentation page.

Finally, note that each of the factorizations takes the current partition of domain and codomain as the axis along which to matricize and perform the factorization. In order to obtain factorizations according to a different bipartition of the indices, we can use any of the previously mentioned index manipulations before the factorization.

Some examples to conclude this section:

julia> V1 = SU₂Space(0 => 2, 1/2 => 1)Rep[SU₂](…) of dim 4:
   0 => 2
 1/2 => 1
julia> V2 = SU₂Space(0 => 1, 1/2 => 1, 1 => 1)Rep[SU₂](…) of dim 6:
   0 => 1
 1/2 => 1
   1 => 1
julia> t = randn(V1 ⊗ V1, V2);
julia> U, S, Vh = svd_compact(t);
julia> t ≈ U * S * Vhtrue
julia> D, V = eigh_full(t' * t);
julia> D ≈ S * Strue
julia> U' * U ≈ id(domain(U))true
julia> S6←6 DiagonalTensorMap{Float64, Rep[SU₂], Vector{Float64}}:
 codomain: ⊗(Rep[SU₂](0 => 1, 1/2 => 1, 1 => 1))
 domain: ⊗(Rep[SU₂](0 => 1, 1/2 => 1, 1 => 1))
 blocks:
 * Irrep[SU₂](0) => 1×1 LinearAlgebra.Diagonal{Float64, SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}}:
 2.4026867079596546

 * Irrep[SU₂](1/2) => 1×1 LinearAlgebra.Diagonal{Float64, SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}}:
 1.052650270259043

 * Irrep[SU₂](1) => 1×1 LinearAlgebra.Diagonal{Float64, SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}}:
 1.0104361741133534
julia> Q, R = left_orth(t; alg = :svd);
julia> Q' * Q ≈ id(domain(Q))true
julia> t ≈ Q * Rtrue
julia> U2, S2, Vh2, ε = svd_trunc(t; trunc = truncspace(V1));
julia> Vh2 * Vh2' ≈ id(codomain(Vh2))true
julia> S23←3 DiagonalTensorMap{Float64, Rep[SU₂], Vector{Float64}}:
 codomain: ⊗(Rep[SU₂](0 => 1, 1/2 => 1))
 domain: ⊗(Rep[SU₂](0 => 1, 1/2 => 1))
 blocks:
 * Irrep[SU₂](0) => 1×1 LinearAlgebra.Diagonal{Float64, SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}}:
 2.4026867079596546

 * Irrep[SU₂](1/2) => 1×1 LinearAlgebra.Diagonal{Float64, SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}}:
 1.052650270259043
julia> ε ≈ norm(block(S, Irrep[SU₂](1))) * sqrt(dim(Irrep[SU₂](1)))true
julia> L, Q = right_orth(permute(t, ((1,), (2, 3))));
julia> codomain(L), domain(L), domain(Q)(⊗(Rep[SU₂](0 => 2, 1/2 => 1)), ⊗(Rep[SU₂](0 => 2, 1/2 => 1)), (Rep[SU₂](0 => 2, 1/2 => 1)' ⊗ Rep[SU₂](0 => 1, 1/2 => 1, 1 => 1)))
julia> Q * Q'4←4 TensorMap{Float64, Rep[SU₂], 1, 1, Vector{Float64}}:
 codomain: ⊗(Rep[SU₂](0 => 2, 1/2 => 1))
 domain: ⊗(Rep[SU₂](0 => 2, 1/2 => 1))
 blocks:
 * Irrep[SU₂](0) => 2×2 reshape(view(::Vector{Float64}, 1:4), 2, 2) with eltype Float64:
 1.0          2.91074e-17
 2.91074e-17  1.0

 * Irrep[SU₂](1/2) => 1×1 reshape(view(::Vector{Float64}, 5:5), 1, 1) with eltype Float64:
 1.0000000000000002
julia> P = Q' * Q;
julia> P ≈ P * Ptrue
julia> t′ = permute(t, ((1,), (2, 3)));
julia> t′ ≈ t′ * Ptrue

To compute only the R or L factor with QR/LQ-based left_orth! and right_orth!, pass a zero-dimensional tensor map as the basis output.

V = infimum(fuse(codomain(t)), fuse(domain(t)))
V0 = zerospace(spacetype(t))
Q0 = similar(t, V0 ← V0)

R = similar(t, V ← domain(t))
_, R = left_orth!(copy(t), (Q0, R))

L = similar(t, codomain(t) ← V)
L, _ = right_orth!(copy(t), (L, Q0))

For a DiagonalTensorMap, use an empty DiagonalTensorMap as the basis output instead. Differentiation of these R/L-only paths is not yet supported.