Updates for MatrixAlgebraKit v0.6 (#312)

* bump MatrixAlgebraKit compat to 0.6

* update `ishermitian` implementations

* update `isisometric` implementations

* add projections

* Incremental fixes

* orthnull progress

* Adjoint factorizations (AD not working)

* add missing MAK.

* remove `eig(::Adjoint)`

This doesn't actually hold for non-hermitian matrices

* add truncation error implementation

* add nullspace truncation specializations

* update orth interface in tests

* fix wrong variable

* improve adjoint support

* remove unnecessary left/rightorth functions

* pass kwargs in isisometri

* add tests for truncation error

* more nullspace tests and fixes

* remove unnecessary specialization

* add missing argument

* remove unnecessary specializations

* slight reorganization

* export new functionality

* simplify truncationerror implementation

* fix testcase

* refrain from defining left_null functions to avoid ambiguities

* fix ad tests

* add projection tests

* temporary AD fixes

* try and add back some AD tests

* apply code suggestions

* stabilize AD test

* some changes

* more updates: formatting and deprecations

* some more review updates

* handle fully-truncated blocks

* fix diagonal adjoint implementation

* fix property name for PolarViaSVD

* copy_input uses `f` instead of `f!`

* remove more implementations to fix code

* remove `check_input`

---------

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

Author: 3 people

Project.toml

@@ -35,7 +35,7 @@ Combinatorics = "1"
 FiniteDifferences = "0.12"
 LRUCache = "1.0.2"
 LinearAlgebra = "1"
-MatrixAlgebraKit = "0.5.0"
+MatrixAlgebraKit = "0.6.0"
 OhMyThreads = "0.8.0"
 PackageExtensionCompat = "1"
 Printf = "1"

src/TensorKit.jl

@@ -79,10 +79,12 @@ export left_orth, right_orth, left_null, right_null,
     qr_full!, qr_compact!, qr_null!, lq_full!, lq_compact!, lq_null!,
     svd_compact!, svd_full!, svd_trunc!, svd_compact, svd_full, svd_trunc,
     exp, exp!,
-    eigh_full!, eigh_full, eigh_trunc!, eigh_trunc, eig_full!, eig_full, eig_trunc!,
-    eig_trunc,
-    eigh_vals!, eigh_vals, eig_vals!, eig_vals,
-    isposdef, isposdef!, ishermitian, isisometry, isunitary, sylvester, rank, cond
+    eigh_full!, eigh_full, eigh_trunc!, eigh_trunc, eigh_vals!, eigh_vals,
+    eig_full!, eig_full, eig_trunc!, eig_trunc, eig_vals!, eig_vals,
+    ishermitian, project_hermitian, project_hermitian!,
+    isantihermitian, project_antihermitian, project_antihermitian!,
+    isisometric, isunitary, project_isometric, project_isometric!,
+    isposdef, isposdef!, sylvester, rank, cond
 
 export braid, braid!, permute, permute!, transpose, transpose!, twist, twist!, repartition,
     repartition!
@@ -135,7 +137,7 @@ using LinearAlgebra: norm, dot, normalize, normalize!, tr,
     adjoint, adjoint!, transpose, transpose!,
     lu, pinv, sylvester,
     eigen, eigen!, svd, svd!,
-    isposdef, isposdef!, ishermitian, rank, cond,
+    isposdef, isposdef!, rank, cond,
     Diagonal, Hermitian
 using MatrixAlgebraKit
 

src/auxiliary/deprecate.jl

@@ -186,7 +186,7 @@ function tsvd(t::AbstractTensorMap; kwargs...)
         Base.depwarn("p is a deprecated kwarg, and should be specified through the truncation strategy", :tsvd)
         kwargs = _drop_p(; kwargs...)
     end
-    return haskey(kwargs, :trunc) ? svd_trunc(t; kwargs...) : svd_compact(t; kwargs...)
+    return haskey(kwargs, :trunc) ? svd_trunc(t; kwargs...) : (svd_compact(t; kwargs...)..., abs(zero(scalartype(t))))
 end
 function tsvd!(t::AbstractTensorMap; kwargs...)
     Base.depwarn("`tsvd!` is deprecated, use `svd_compact!`, `svd_full!` or `svd_trunc!` instead", :tsvd!)

src/factorizations/adjoint.jl

@@ -6,46 +6,54 @@ _adjoint(alg::MAK.LAPACK_HouseholderQR) = MAK.LAPACK_HouseholderLQ(; alg.kwargs.
 _adjoint(alg::MAK.LAPACK_HouseholderLQ) = MAK.LAPACK_HouseholderQR(; alg.kwargs...)
 _adjoint(alg::MAK.LAPACK_HouseholderQL) = MAK.LAPACK_HouseholderRQ(; alg.kwargs...)
 _adjoint(alg::MAK.LAPACK_HouseholderRQ) = MAK.LAPACK_HouseholderQL(; alg.kwargs...)
-_adjoint(alg::MAK.PolarViaSVD) = MAK.PolarViaSVD(_adjoint(alg.svdalg))
+_adjoint(alg::MAK.PolarViaSVD) = MAK.PolarViaSVD(_adjoint(alg.svd_alg))
 _adjoint(alg::AbstractAlgorithm) = alg
 
-# 1-arg functions
-function MAK.initialize_output(::typeof(left_null!), t::AdjointTensorMap, alg::AbstractAlgorithm)
-    return adjoint(MAK.initialize_output(right_null!, adjoint(t), _adjoint(alg)))
-end
-function MAK.initialize_output(
-        ::typeof(right_null!), t::AdjointTensorMap,
-        alg::AbstractAlgorithm
-    )
-    return adjoint(MAK.initialize_output(left_null!, adjoint(t), _adjoint(alg)))
+for f in
+    [
+        :svd_compact, :svd_full, :svd_vals,
+        :qr_compact, :qr_full, :qr_null,
+        :lq_compact, :lq_full, :lq_null,
+        :eig_full, :eig_vals, :eigh_full, :eigh_trunc, :eigh_vals,
+        :left_polar, :right_polar,
+        :project_hermitian, :project_antihermitian, :project_isometric,
+    ]
+    f! = Symbol(f, :!)
+    # just return the algorithm for the parent type since we are mapping this with
+    # `_adjoint` afterwards anyways.
+    # TODO: properly handle these cases
+    @eval MAK.default_algorithm(::typeof($f!), ::Type{T}; kwargs...) where {T <: AdjointTensorMap} =
+        MAK.default_algorithm($f!, TensorKit.parenttype(T); kwargs...)
 end
 
-function MAK.left_null!(t::AdjointTensorMap, N, alg::AbstractAlgorithm)
-    right_null!(adjoint(t), adjoint(N), _adjoint(alg))
-    return N
-end
-function MAK.right_null!(t::AdjointTensorMap, N, alg::AbstractAlgorithm)
-    left_null!(adjoint(t), adjoint(N), _adjoint(alg))
-    return N
-end
+# 1-arg functions
+MAK.initialize_output(::typeof(qr_null!), t::AdjointTensorMap, alg::AbstractAlgorithm) =
+    adjoint(MAK.initialize_output(lq_null!, adjoint(t), _adjoint(alg)))
+MAK.initialize_output(::typeof(lq_null!), t::AdjointTensorMap, alg::AbstractAlgorithm) =
+    adjoint(MAK.initialize_output(qr_null!, adjoint(t), _adjoint(alg)))
 
-function MAK.is_left_isometry(t::AdjointTensorMap; kwargs...)
-    return is_right_isometry(adjoint(t); kwargs...)
-end
-function MAK.is_right_isometry(t::AdjointTensorMap; kwargs...)
-    return is_left_isometry(adjoint(t); kwargs...)
-end
+MAK.qr_null!(t::AdjointTensorMap, N, alg::AbstractAlgorithm) =
+    lq_null!(adjoint(t), adjoint(N), _adjoint(alg))
+MAK.lq_null!(t::AdjointTensorMap, N, alg::AbstractAlgorithm) =
+    qr_null!(adjoint(t), adjoint(N), _adjoint(alg))
+
+MAK.is_left_isometric(t::AdjointTensorMap; kwargs...) =
+    MAK.is_right_isometric(adjoint(t); kwargs...)
+MAK.is_right_isometric(t::AdjointTensorMap; kwargs...) =
+    MAK.is_left_isometric(adjoint(t); kwargs...)
 
 # 2-arg functions
-for (left_f!, right_f!) in zip(
-        (:qr_full!, :qr_compact!, :left_polar!, :left_orth!),
-        (:lq_full!, :lq_compact!, :right_polar!, :right_orth!)
+for (left_f, right_f) in zip(
+        (:qr_full, :qr_compact, :left_polar),
+        (:lq_full, :lq_compact, :right_polar)
     )
-    @eval function MAK.copy_input(::typeof($left_f!), t::AdjointTensorMap)
-        return adjoint(MAK.copy_input($right_f!, adjoint(t)))
+    left_f! = Symbol(left_f, :!)
+    right_f! = Symbol(right_f, :!)
+    @eval function MAK.copy_input(::typeof($left_f), t::AdjointTensorMap)
+        return adjoint(MAK.copy_input($right_f, adjoint(t)))
     end
-    @eval function MAK.copy_input(::typeof($right_f!), t::AdjointTensorMap)
-        return adjoint(MAK.copy_input($left_f!, adjoint(t)))
+    @eval function MAK.copy_input(::typeof($right_f), t::AdjointTensorMap)
+        return adjoint(MAK.copy_input($left_f, adjoint(t)))
     end
 
     @eval function MAK.initialize_output(
@@ -60,29 +68,31 @@ for (left_f!, right_f!) in zip(
     end
 
     @eval function MAK.$left_f!(t::AdjointTensorMap, F, alg::AbstractAlgorithm)
-        $right_f!(adjoint(t), reverse(adjoint.(F)), _adjoint(alg))
-        return F
+        F′ = $right_f!(adjoint(t), reverse(adjoint.(F)), _adjoint(alg))
+        return reverse(adjoint.(F′))
     end
     @eval function MAK.$right_f!(t::AdjointTensorMap, F, alg::AbstractAlgorithm)
-        $left_f!(adjoint(t), reverse(adjoint.(F)), _adjoint(alg))
-        return F
+        F′ = $left_f!(adjoint(t), reverse(adjoint.(F)), _adjoint(alg))
+        return reverse(adjoint.(F′))
     end
 end
 
 # 3-arg functions
-for f! in (:svd_full!, :svd_compact!, :svd_trunc!)
-    @eval function MAK.copy_input(::typeof($f!), t::AdjointTensorMap)
-        return adjoint(MAK.copy_input($f!, adjoint(t)))
+for f in (:svd_full, :svd_compact)
+    f! = Symbol(f, :!)
+    @eval function MAK.copy_input(::typeof($f), t::AdjointTensorMap)
+        return adjoint(MAK.copy_input($f, adjoint(t)))
     end
 
     @eval function MAK.initialize_output(
             ::typeof($f!), t::AdjointTensorMap, alg::AbstractAlgorithm
         )
         return reverse(adjoint.(MAK.initialize_output($f!, adjoint(t), _adjoint(alg))))
     end
+
     @eval function MAK.$f!(t::AdjointTensorMap, F, alg::AbstractAlgorithm)
-        $f!(adjoint(t), reverse(adjoint.(F)), _adjoint(alg))
-        return F
+        F′ = $f!(adjoint(t), reverse(adjoint.(F)), _adjoint(alg))
+        return reverse(adjoint.(F′))
     end
 
     # disambiguate by prohibition
@@ -92,17 +102,9 @@ for f! in (:svd_full!, :svd_compact!, :svd_trunc!)
         throw(MethodError($f!, (t, alg)))
     end
 end
+
 # avoid amgiguity
-function MAK.initialize_output(
-        ::typeof(svd_trunc!), t::AdjointTensorMap, alg::TruncatedAlgorithm
-    )
-    return MAK.initialize_output(svd_compact!, t, alg.alg)
-end
-# to fix ambiguity
-function MAK.svd_trunc!(t::AdjointTensorMap, USVᴴ, alg::TruncatedAlgorithm)
-    USVᴴ′ = svd_compact!(t, USVᴴ, alg.alg)
-    return MAK.truncate(svd_trunc!, USVᴴ′, alg.trunc)
-end
-function MAK.svd_compact!(t::AdjointTensorMap, USVᴴ, alg::DiagonalAlgorithm)
-    return MAK.svd_compact!(t, USVᴴ, alg.alg)
+function MAK.svd_compact!(t::AdjointTensorMap, F, alg::DiagonalAlgorithm)
+    F′ = svd_compact!(adjoint(t), reverse(adjoint.(F)), _adjoint(alg))
+    return reverse(adjoint.(F′))
 end

src/factorizations/diagonal.jl

@@ -81,55 +81,20 @@ for f! in
         :eigh_trunc!, :right_orth!, :left_orth!,
     )
     @eval function MAK.$f!(d::DiagonalTensorMap, F, alg::DiagonalAlgorithm)
-        MAK.check_input($f!, d, F, alg)
         $f!(_repack_diagonal(d), _repack_diagonal.(F), alg)
         return F
     end
 end
 
-for f! in (:qr_full!, :qr_compact!)
-    @eval function MAK.check_input(
-            ::typeof($f!), d::AbstractTensorMap, QR, ::DiagonalAlgorithm
-        )
-        Q, R = QR
-        @assert d isa DiagonalTensorMap
-        @assert Q isa DiagonalTensorMap && R isa DiagonalTensorMap
-        @check_scalar Q d
-        @check_scalar R d
-        @check_space(Q, space(d))
-        @check_space(R, space(d))
-
-        return nothing
-    end
-end
-
-for f! in (:lq_full!, :lq_compact!)
-    @eval function MAK.check_input(
-            ::typeof($f!), d::AbstractTensorMap, LQ, ::DiagonalAlgorithm
-        )
-        L, Q = LQ
-        @assert d isa DiagonalTensorMap
-        @assert Q isa DiagonalTensorMap && L isa DiagonalTensorMap
-        @check_scalar Q d
-        @check_scalar L d
-        @check_space(Q, space(d))
-        @check_space(L, space(d))
-
-        return nothing
-    end
-end
-
 # disambiguate
 function MAK.svd_compact!(t::AbstractTensorMap, USVᴴ, alg::DiagonalAlgorithm)
     return svd_full!(t, USVᴴ, alg)
 end
 
 # f_vals
 # ------
-
 for f! in (:eig_vals!, :eigh_vals!, :svd_vals!)
     @eval function MAK.$f!(d::AbstractTensorMap, V, alg::DiagonalAlgorithm)
-        MAK.check_input($f!, d, V, alg)
         $f!(_repack_diagonal(d), diagview(_repack_diagonal(V)), alg)
         return V
     end
@@ -140,64 +105,3 @@ for f! in (:eig_vals!, :eigh_vals!, :svd_vals!)
         return DiagonalTensorMap(data, d.domain)
     end
 end
-
-function MAK.check_input(::typeof(eig_full!), t::AbstractTensorMap, DV, ::DiagonalAlgorithm)
-    domain(t) == codomain(t) ||
-        throw(ArgumentError("Eigenvalue decomposition requires square input tensor"))
-
-    D, V = DV
-
-    @assert D isa DiagonalTensorMap
-    @assert V isa AbstractTensorMap
-
-    # scalartype checks
-    @check_scalar D t
-    @check_scalar V t
-
-    # space checks
-    @check_space D space(t)
-    @check_space V space(t)
-
-    return nothing
-end
-
-function MAK.check_input(::typeof(eigh_full!), t::AbstractTensorMap, DV, ::DiagonalAlgorithm)
-    domain(t) == codomain(t) ||
-        throw(ArgumentError("Eigenvalue decomposition requires square input tensor"))
-
-    D, V = DV
-
-    @assert D isa DiagonalTensorMap
-    @assert V isa AbstractTensorMap
-
-    # scalartype checks
-    @check_scalar D t real
-    @check_scalar V t
-
-    # space checks
-    @check_space D space(t)
-    @check_space V space(t)
-
-    return nothing
-end
-
-function MAK.check_input(::typeof(eig_vals!), t::AbstractTensorMap, D, ::DiagonalAlgorithm)
-    @assert D isa DiagonalTensorMap
-    @check_scalar D t
-    @check_space D space(t)
-    return nothing
-end
-
-function MAK.check_input(::typeof(eigh_vals!), t::AbstractTensorMap, D, ::DiagonalAlgorithm)
-    @assert D isa DiagonalTensorMap
-    @check_scalar D t real
-    @check_space D space(t)
-    return nothing
-end
-
-function MAK.check_input(::typeof(svd_vals!), t::AbstractTensorMap, D, ::DiagonalAlgorithm)
-    @assert D isa DiagonalTensorMap
-    @check_scalar D t real
-    @check_space D space(t)
-    return nothing
-end

src/factorizations/factorizations.jl

@@ -18,8 +18,7 @@ import MatrixAlgebraKit as MAK
 using MatrixAlgebraKit: AbstractAlgorithm, TruncatedAlgorithm, DiagonalAlgorithm
 using MatrixAlgebraKit: TruncationStrategy, NoTruncation, TruncationByValue,
     TruncationByError, TruncationIntersection, TruncationByFilter, TruncationByOrder
-using MatrixAlgebraKit: left_orth_polar!, right_orth_polar!, left_orth_svd!,
-    right_orth_svd!, left_null_svd!, right_null_svd!, diagview
+using MatrixAlgebraKit: diagview
 
 include("utility.jl")
 include("matrixalgebrakit.jl")
@@ -30,11 +29,6 @@ include("pullbacks.jl")
 
 TensorKit.one!(A::AbstractMatrix) = MatrixAlgebraKit.one!(A)
 
-function MatrixAlgebraKit.isisometry(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple)
-    t = permute(t, (p₁, p₂); copy = false)
-    return isisometry(t)
-end
-
 #------------------------------#
 # LinearAlgebra overloads
 #------------------------------#
@@ -61,38 +55,42 @@ LinearAlgebra.svdvals!(t::AbstractTensorMap) = diagview(svd_vals!(t))
 #--------------------------------------------------#
 # Checks for hermiticity and positive definiteness #
 #--------------------------------------------------#
-function LinearAlgebra.ishermitian(t::AbstractTensorMap)
-    domain(t) == codomain(t) || return false
-    InnerProductStyle(t) === EuclideanInnerProduct() || return false # hermiticity only defined for euclidean
-    for (c, b) in blocks(t)
-        ishermitian(b) || return false
+function _blockmap(f; kwargs...)
+    return function ((c, b))
+        return f(b; kwargs...)
     end
-    return true
 end
 
+function MAK.ishermitian(t::AbstractTensorMap; kwargs...)
+    return InnerProductStyle(t) === EuclideanInnerProduct() &&
+        domain(t) == codomain(t) &&
+        all(_blockmap(MAK.ishermitian; kwargs...), blocks(t))
+end
+function MAK.isantihermitian(t::AbstractTensorMap; kwargs...)
+    return InnerProductStyle(t) === EuclideanInnerProduct() &&
+        domain(t) == codomain(t) &&
+        all(_blockmap(MAK.isantihermitian; kwargs...), blocks(t))
+end
+LinearAlgebra.ishermitian(t::AbstractTensorMap) = MAK.ishermitian(t)
+
 function LinearAlgebra.isposdef(t::AbstractTensorMap)
     return isposdef!(copy_oftype(t, factorisation_scalartype(isposdef, t)))
 end
 function LinearAlgebra.isposdef!(t::AbstractTensorMap)
     domain(t) == codomain(t) ||
         throw(SpaceMismatch("`isposdef` requires domain and codomain to be the same"))
     InnerProductStyle(spacetype(t)) === EuclideanInnerProduct() || return false
-    for (c, b) in blocks(t)
-        isposdef!(b) || return false
-    end
-    return true
+    return all(_blockmap(isposdef!), blocks(t))
 end
 
 # TODO: tolerances are per-block, not global or weighted - does that matter?
-function MatrixAlgebraKit.is_left_isometry(t::AbstractTensorMap; kwargs...)
+function MAK.is_left_isometric(t::AbstractTensorMap; kwargs...)
     domain(t) ≾ codomain(t) || return false
-    f((c, b)) = MatrixAlgebraKit.is_left_isometry(b; kwargs...)
-    return all(f, blocks(t))
+    return all(_blockmap(MAK.is_left_isometric; kwargs...), blocks(t))
 end
-function MatrixAlgebraKit.is_right_isometry(t::AbstractTensorMap; kwargs...)
+function MAK.is_right_isometric(t::AbstractTensorMap; kwargs...)
     domain(t) ≿ codomain(t) || return false
-    f((c, b)) = MatrixAlgebraKit.is_right_isometry(b; kwargs...)
-    return all(f, blocks(t))
+    return all(_blockmap(MAK.is_right_isometric; kwargs...), blocks(t))
 end
 
 end