MatrixAlgebraKit decompositions (#230)

This PR is part of a large set of changes that alters the backend from `TensorKit.MatrixAlgebra` to `MatrixAlgebraKit` for the factorizations.
Concretely, this is achieved through a new submodule `TensorKit.Factorizations` that implements most of the `MatrixAlgebraKit` interface, along with some overloads and reimplementations to be somewhat backwards compatible.

## New features

- `foreachblock(f, t::AbstractTensorMap...)` is a new function that is meant to provide a uniform interface to iterate through the blocks of a collection of tensors. The main purpose is to centralize this concept to easily add multithreading over the blocks. In order to break up this PR the actual multithreading implementation is left for a follow-up PR.
- Truncated eigenvalue decompositions through `eig_trunc` and `eigh_trunc`
- The factorizations can now swap out their backends to select alternate algorithms or implementations.
- `ominus` and its unicode variant can now be used to obtain the orthogonal complement of a space, ie if `W = oplus(V1, V2)`, `V2 = ominus(W, V1)`.
 
## Breaking changes

- The `OrthogonalFactorization` structs have been removed, but their constructors are deprecated to obtain equivalent `MatrixAlgebraKit` algorithm structs.
- The factorization functions `leftorth`, `rightorth`, `tsvd`, `eig`, `eigh` have been deprecated in favor of the MatrixAlgebraKit variants.
- The truncation strategies are now directly inherited from `MatrixAlgebraKit`. This means `truncdim` is replaced by `truncrank`, `trunctol` is replaced by `truncbelow`, and there are deprecation warnings for the old functions.
- The `left_orth` and `right_orth` functions will now always output tensors with a single space connecting them, which was not the case for `Polar` decompositions previously. To retrieve the old behavior with `isposdef` `R` factors (equal domain and codomain instead of isomorphic), the functions `left_polar` and `right_polar` can be used.
- The `left_orth` and `right_orth` functions will now always have a connecting space with `isdual=false`. This is different from the previous behavior where some `N,1` or `1,N` tensors kept the `isdual` flag on the connecting leg.
- There is no more direct call to both permute and factorize in one go. This is mainly because this is incompatible with the distinction between `permute` and `braid`.

---------

Co-authored-by: Katharine Hyatt <[email protected]>
Co-authored-by: Katharine Hyatt <[email protected]>
Co-authored-by: Jutho <[email protected]>

Author: 4 people

Project.toml

@@ -1,13 +1,16 @@
 name = "TensorKit"
 uuid = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec"
 authors = ["Jutho Haegeman"]
-version = "0.14.11"
+version = "0.15.0"
 
 [deps]
 LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
+OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5"
 PackageExtensionCompat = "65ce6f38-6b18-4e1d-a461-8949797d7930"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+ScopedValues = "7e506255-f358-4e82-b7e4-beb19740aa63"
 Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67"
 TensorKitSectors = "13a9c161-d5da-41f0-bcbd-e1a08ae0647f"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
@@ -30,8 +33,11 @@ Combinatorics = "1"
 FiniteDifferences = "0.12"
 LRUCache = "1.0.2"
 LinearAlgebra = "1"
+MatrixAlgebraKit = "0.5.0"
+OhMyThreads = "0.8.0"
 PackageExtensionCompat = "1"
 Random = "1"
+ScopedValues = "1.3.0"
 Strided = "2"
 TensorKitSectors = "0.1.4, 0.2"
 TensorOperations = "5.1"

docs/Project.toml

@@ -1,10 +1,9 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec"
-TensorKitSectors = "13a9c161-d5da-41f0-bcbd-e1a08ae0647f"
 
 [compat]
 Documenter = "1"
 Random = "1"
-TensorKitSectors = "0.1"

docs/make.jl

@@ -1,6 +1,12 @@
 using Documenter
 using Random
-using TensorKit, TensorKitSectors
+using TensorKit
+using TensorKit.TensorKitSectors
+using TensorKit.MatrixAlgebraKit
+using DocumenterInterLinks
+
+links = InterLinks("MatrixAlgebraKit" => "https://quantumkithub.github.io/MatrixAlgebraKit.jl/stable/",
+                   "TensorOperations" => "https://quantumkithub.github.io/TensorOperations.jl/stable/")
 
 pages = ["Home" => "index.md",
          "Manual" => ["man/intro.md", "man/tutorial.md", "man/categories.md",
@@ -15,6 +21,7 @@ makedocs(; modules=[TensorKit, TensorKitSectors],
          format=Documenter.HTML(; prettyurls=true, mathengine=MathJax(),
                                 assets=["assets/custom.css"]),
          pages=pages,
-         pagesonly=true)
+         pagesonly=true,
+         plugins=[links])
 
 deploydocs(; repo="github.com/QuantumKitHub/TensorKit.jl.git", push_preview=true)

docs/src/lib/spaces.md

@@ -90,6 +90,7 @@ dual
 conj
 flip
 ⊕
+⊖
 zero(::ElementarySpace)
 oneunit
 supremum

docs/src/lib/tensors.md

@@ -214,26 +214,21 @@ contract!
 
 ## `TensorMap` factorizations
 
-The factorisation methods come in two flavors, namely a non-destructive version where you
-can specify an additional permutation of the domain and codomain indices before the
-factorisation is performed (provided that `sectorstyle(t)` has a symmetric braiding) as
-well as a destructive version The non-destructive methods are given first:
+The factorisation methods are powered by [MatrixAlgebraKit.jl](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl)
+and all follow the same strategy. The idea is that the `TensorMap` is interpreted as a linear
+map based on the current partition of indices between `domain` and `codomain`, and then the
+entire range of MatrixAlgebraKit functions can be called.
+Factorizing a tensor according to a different partition of the indices is possible
+by prepending the factorization step with an explicit call to [`permute`](@ref) or [`transpose`](@ref).
 
-```@docs
-leftorth
-rightorth
-leftnull
-rightnull
-tsvd
-eigh
-eig
-eigen
-isposdef
-```
+For the full list of factorizations, see [Decompositions](@extref MatrixAlgebraKit).
 
-The corresponding destructive methods have an exclamation mark at the end of their name,
-and only accept the `TensorMap` object as well as the method-specific algorithm and keyword
-arguments.
+Additionally, it is possible to obtain truncated versions of some of these factorizations
+through the [`MatrixAlgebraKit.TruncationStrategy`](@ref) objects.
 
+The exact truncation strategy can be controlled through the strategies defined in [Truncations](@extref MatrixAlgebraKit),
+but for `TensorMap`s there is also the special-purpose scheme:
 
-TODO: document svd truncation types
+```@docs
+truncspace
+```

ext/TensorKitChainRulesCoreExt/TensorKitChainRulesCoreExt.jl

@@ -3,6 +3,7 @@ module TensorKitChainRulesCoreExt
 using TensorOperations
 using VectorInterface
 using TensorKit
+using TensorKit: foreachblock
 using ChainRulesCore
 using LinearAlgebra
 using TupleTools
@@ -15,6 +16,5 @@ include("utility.jl")
 include("constructors.jl")
 include("linalg.jl")
 include("tensoroperations.jl")
-include("factorizations.jl")
 
 end