diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 7f94f9373..e7be5d8e6 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -28,7 +28,7 @@ jobs: - macOS-latest - windows-latest steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} @@ -57,7 +57,7 @@ jobs: - windows-latest continue-on-error: true steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index 4908988c9..82540cbdf 100644 --- a/.github/workflows/Documentation.yml +++ b/.github/workflows/Documentation.yml @@ -20,7 +20,7 @@ jobs: os: - ubuntu-latest steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - uses: julia-actions/setup-julia@latest with: version: ${{ matrix.version }} diff --git a/.github/workflows/DocumentationCleanup.yml b/.github/workflows/DocumentationCleanup.yml index 5be23b964..f8a1e525e 100644 --- a/.github/workflows/DocumentationCleanup.yml +++ b/.github/workflows/DocumentationCleanup.yml @@ -16,7 +16,7 @@ jobs: contents: write steps: - name: Checkout gh-pages branch - uses: actions/checkout@v4 + uses: actions/checkout@v5 with: ref: gh-pages - name: Delete preview and history + push changes diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml index aa23cd6a2..e112b3028 100644 --- a/.github/workflows/FormatCheck.yml +++ b/.github/workflows/FormatCheck.yml @@ -23,7 +23,7 @@ jobs: with: version: ${{ matrix.version }} - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - name: Install JuliaFormatter and format # This will use the latest version by default but you can set the version like so: # diff --git a/CITATION.cff b/CITATION.cff new file mode 100644 index 000000000..628a9c587 --- /dev/null +++ b/CITATION.cff @@ -0,0 +1,27 @@ +cff-version: 1.2.0 +message: "If you use this software, please cite it as below." +authors: +- family-names: "Devos" + given-names: "Lukas" + orcid: "https://orcid.org/0000-0002-0256-4200" +- family-names: "Haegeman" + given-names: "Jutho" + orcid: "https://orcid.org/0000-0002-0858-291X" +title: "TensorKit.jl" +version: "0.14.9" +doi: "10.5281/zenodo.8421339" +date-released: "2025-07-18" +url: "https://github.com/Jutho/TensorKit.jl" +preferred-citation: + type: article + authors: + - family-names: "Devos" + given-names: "Lukas" + orcid: "https://orcid.org/0000-0002-0256-4200" + - family-names: "Haegeman" + given-names: "Jutho" + orcid: "https://orcid.org/0000-0002-0858-291X" + doi: "10.48550/arXiv.2508.10076" + journal: "arXiv" + title: "TensorKit.jl: A Julia package for large-scale tensor computations, with a hint of category theory" + year: 2025 diff --git a/Project.toml b/Project.toml index b063816f7..a1847b6dd 100644 --- a/Project.toml +++ b/Project.toml @@ -6,8 +6,11 @@ version = "0.14.9" [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.3" +OhMyThreads = "0.8.0" PackageExtensionCompat = "1" Random = "1" +ScopedValues = "1.3.0" Strided = "2" TensorKitSectors = "0.1" TensorOperations = "5.1" diff --git a/README.md b/README.md index 835982c6a..e3acc2c9c 100644 --- a/README.md +++ b/README.md @@ -153,6 +153,7 @@ Check out the [tutorial](https://jutho.github.io/TensorKit.jl/stable/man/tutoria full [documentation](https://jutho.github.io/TensorKit.jl/stable). ## Installation + `TensorKit.jl` can be installed with the Julia package manager. From the Julia REPL, type `]` to enter the Pkg REPL mode and run: ``` @@ -180,3 +181,9 @@ platforms with a 64-bit architecture. Contributions are very welcome, as are feature requests and suggestions. Please open an [issue][issues-url] if you encounter any problems. [issues-url]: https://github.com/Jutho/TensorKit.jl/issues + +## Acknowledgements + +The design and development of the TensorKit.jl package have benefited from countless discussions with many people, including most current and former members of the Quantum Group at Ghent University. +Being an open-source software project developed over the course of many years, we also thank all past, current and future contributors, including the bug reports and feature requests that have shaped this package. +In particular, we like to thank [Maarten Van Damme](@maartenvd), who initiated the MPSKit.jl package on top of TensorKit.jl early-on, and has as such had a strong influence on the development and design decisions of the TensorKit.jl package. diff --git a/docs/src/lib/spaces.md b/docs/src/lib/spaces.md index 83350156d..205301d83 100644 --- a/docs/src/lib/spaces.md +++ b/docs/src/lib/spaces.md @@ -90,6 +90,7 @@ dual conj flip ⊕ +⊖ zero(::ElementarySpace) oneunit supremum diff --git a/ext/TensorKitChainRulesCoreExt/TensorKitChainRulesCoreExt.jl b/ext/TensorKitChainRulesCoreExt/TensorKitChainRulesCoreExt.jl index 16c7583d1..36bb4108d 100644 --- a/ext/TensorKitChainRulesCoreExt/TensorKitChainRulesCoreExt.jl +++ b/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 @@ -11,6 +12,11 @@ import TensorOperations as TO using TensorOperations: promote_contract, tensoralloc_add, tensoralloc_contract using VectorInterface: promote_scale, promote_add +using MatrixAlgebraKit +using MatrixAlgebraKit: TruncationStrategy, + svd_compact_pullback!, eig_full_pullback!, eigh_full_pullback!, + qr_compact_pullback!, lq_compact_pullback! + include("utility.jl") include("constructors.jl") include("linalg.jl") diff --git a/ext/TensorKitChainRulesCoreExt/factorizations.jl b/ext/TensorKitChainRulesCoreExt/factorizations.jl index abfb724bc..0c6924411 100644 --- a/ext/TensorKitChainRulesCoreExt/factorizations.jl +++ b/ext/TensorKitChainRulesCoreExt/factorizations.jl @@ -1,42 +1,38 @@ # Factorizations rules # -------------------- -function ChainRulesCore.rrule(::typeof(TensorKit.tsvd!), t::AbstractTensorMap; - trunc::TensorKit.TruncationScheme=TensorKit.NoTruncation(), - p::Real=2, - alg::Union{TensorKit.SVD,TensorKit.SDD}=TensorKit.SDD()) - U, Σ, V⁺, truncerr = tsvd(t; trunc=TensorKit.NoTruncation(), p=p, alg=alg) - - if !(trunc isa TensorKit.NoTruncation) && !isempty(blocksectors(t)) - Σdata = TensorKit.SectorDict(c => diag(b) for (c, b) in blocks(Σ)) - - truncdim = TensorKit._compute_truncdim(Σdata, trunc, p) - truncerr = TensorKit._compute_truncerr(Σdata, truncdim, p) - - SVDdata = TensorKit.SectorDict(c => (block(U, c), Σc, block(V⁺, c)) - for (c, Σc) in Σdata) - - Ũ, Σ̃, Ṽ⁺ = TensorKit._create_svdtensors(t, SVDdata, truncdim) - else - Ũ, Σ̃, Ṽ⁺ = U, Σ, V⁺ - end +for f in (:tsvd, :eig, :eigh) + f! = Symbol(f, :!) + f_trunc! = f == :tsvd ? :svd_trunc! : Symbol(f, :_trunc!) + f_pullback = Symbol(f, :_pullback) + f_pullback! = f == :tsvd ? :svd_compact_pullback! : Symbol(f, :_full_pullback!) + @eval function ChainRulesCore.rrule(::typeof(TensorKit.$f!), t::AbstractTensorMap; + trunc::TruncationStrategy=TensorKit.notrunc(), + kwargs...) + # TODO: I think we can use f! here without issues because we don't actually require + # the data of `t` anymore. + F = $f(t; trunc=TensorKit.notrunc(), kwargs...) + + if trunc != TensorKit.notrunc() && !isempty(blocksectors(t)) + F′ = MatrixAlgebraKit.truncate!($f_trunc!, F, trunc) + else + F′ = F + end - function tsvd!_pullback(ΔUSVϵ) - ΔU, ΔΣ, ΔV⁺, = unthunk.(ΔUSVϵ) - Δt = similar(t) - for (c, b) in blocks(Δt) - Uc, Σc, V⁺c = block(U, c), block(Σ, c), block(V⁺, c) - ΔUc, ΔΣc, ΔV⁺c = block(ΔU, c), block(ΔΣ, c), block(ΔV⁺, c) - Σdc = view(Σc, diagind(Σc)) - ΔΣdc = (ΔΣc isa AbstractZero) ? ΔΣc : view(ΔΣc, diagind(ΔΣc)) - svd_pullback!(b, Uc, Σdc, V⁺c, ΔUc, ΔΣdc, ΔV⁺c) + function $f_pullback(ΔF′) + ΔF = unthunk.(ΔF′) + Δt = zerovector(t) + foreachblock(Δt) do c, (b,) + Fc = block.(F, Ref(c)) + ΔFc = block.(ΔF, Ref(c)) + $f_pullback!(b, Fc, ΔFc) + return nothing + end + return NoTangent(), Δt end - return NoTangent(), Δt - end - function tsvd!_pullback(::Tuple{ZeroTangent,ZeroTangent,ZeroTangent}) - return NoTangent(), ZeroTangent() - end + $f_pullback(::Tuple{ZeroTangent,Vararg{ZeroTangent}}) = NoTangent(), ZeroTangent() - return (Ũ, Σ̃, Ṽ⁺, truncerr), tsvd!_pullback + return F′, $f_pullback + end end function ChainRulesCore.rrule(::typeof(LinearAlgebra.svdvals!), t::AbstractTensorMap) @@ -53,50 +49,6 @@ function ChainRulesCore.rrule(::typeof(LinearAlgebra.svdvals!), t::AbstractTenso return s, svdvals_pullback end -function ChainRulesCore.rrule(::typeof(TensorKit.eig!), t::AbstractTensorMap; kwargs...) - D, V = eig(t; kwargs...) - - function eig!_pullback((_ΔD, _ΔV)) - ΔD, ΔV = unthunk(_ΔD), unthunk(_ΔV) - Δt = similar(t) - for (c, b) in blocks(Δt) - Dc, Vc = block(D, c), block(V, c) - ΔDc, ΔVc = block(ΔD, c), block(ΔV, c) - Ddc = view(Dc, diagind(Dc)) - ΔDdc = (ΔDc isa AbstractZero) ? ΔDc : view(ΔDc, diagind(ΔDc)) - eig_pullback!(b, Ddc, Vc, ΔDdc, ΔVc) - end - return NoTangent(), Δt - end - function eig!_pullback(::Tuple{ZeroTangent,ZeroTangent}) - return NoTangent(), ZeroTangent() - end - - return (D, V), eig!_pullback -end - -function ChainRulesCore.rrule(::typeof(TensorKit.eigh!), t::AbstractTensorMap; kwargs...) - D, V = eigh(t; kwargs...) - - function eigh!_pullback((_ΔD, _ΔV)) - ΔD, ΔV = unthunk(_ΔD), unthunk(_ΔV) - Δt = similar(t) - for (c, b) in blocks(Δt) - Dc, Vc = block(D, c), block(V, c) - ΔDc, ΔVc = block(ΔD, c), block(ΔV, c) - Ddc = view(Dc, diagind(Dc)) - ΔDdc = (ΔDc isa AbstractZero) ? ΔDc : view(ΔDc, diagind(ΔDc)) - eigh_pullback!(b, Ddc, Vc, ΔDdc, ΔVc) - end - return NoTangent(), Δt - end - function eigh!_pullback(::Tuple{ZeroTangent,ZeroTangent}) - return NoTangent(), ZeroTangent() - end - - return (D, V), eigh!_pullback -end - function ChainRulesCore.rrule(::typeof(LinearAlgebra.eigvals!), t::AbstractTensorMap; sortby=nothing, kwargs...) @assert sortby === nothing "only `sortby=nothing` is supported" @@ -115,367 +67,38 @@ end function ChainRulesCore.rrule(::typeof(leftorth!), t::AbstractTensorMap; alg=QRpos()) alg isa TensorKit.QR || alg isa TensorKit.QRpos || error("only `alg=QR()` and `alg=QRpos()` are supported") - Q, R = leftorth(t; alg) - function leftorth!_pullback((_ΔQ, _ΔR)) - ΔQ, ΔR = unthunk(_ΔQ), unthunk(_ΔR) - Δt = similar(t) - for (c, b) in blocks(Δt) - qr_pullback!(b, block(Q, c), block(R, c), block(ΔQ, c), block(ΔR, c)) + QR = leftorth(t; alg) + function leftorth!_pullback(ΔQR′) + ΔQR = unthunk.(ΔQR′) + Δt = zerovector(t) + foreachblock(Δt) do c, (b,) + QRc = block.(QR, Ref(c)) + ΔQRc = block.(ΔQR, Ref(c)) + qr_compact_pullback!(b, QRc, ΔQRc) + return nothing end return NoTangent(), Δt end - leftorth!_pullback(::Tuple{ZeroTangent,ZeroTangent}) = NoTangent(), ZeroTangent() - return (Q, R), leftorth!_pullback + leftorth!_pullback(::NTuple{2,ZeroTangent}) = NoTangent(), ZeroTangent() + + return QR, leftorth!_pullback end function ChainRulesCore.rrule(::typeof(rightorth!), t::AbstractTensorMap; alg=LQpos()) alg isa TensorKit.LQ || alg isa TensorKit.LQpos || error("only `alg=LQ()` and `alg=LQpos()` are supported") - L, Q = rightorth(t; alg) - function rightorth!_pullback((_ΔL, _ΔQ)) - ΔL, ΔQ = unthunk(_ΔL), unthunk(_ΔQ) - Δt = similar(t) - for (c, b) in blocks(Δt) - lq_pullback!(b, block(L, c), block(Q, c), block(ΔL, c), block(ΔQ, c)) + LQ = rightorth(t; alg) + function rightorth!_pullback(ΔLQ′) + ΔLQ = unthunk(ΔLQ′) + Δt = zerovector(t) + foreachblock(Δt) do c, (b,) + LQc = block.(LQ, Ref(c)) + ΔLQc = block.(ΔLQ, Ref(c)) + lq_compact_pullback!(b, LQc, ΔLQc) + return nothing end return NoTangent(), Δt end - rightorth!_pullback(::Tuple{ZeroTangent,ZeroTangent}) = NoTangent(), ZeroTangent() - return (L, Q), rightorth!_pullback -end - -# Corresponding matrix factorisations: implemented as mutating methods -# --------------------------------------------------------------------- -# helper routines -safe_inv(a, tol) = abs(a) < tol ? zero(a) : inv(a) - -function lowertriangularind(A::AbstractMatrix) - m, n = size(A) - I = Vector{Int}(undef, div(m * (m - 1), 2) + m * (n - m)) - offset = 0 - for j in 1:n - r = (j + 1):m - I[offset .- j .+ r] = (j - 1) * m .+ r - offset += length(r) - end - return I -end - -function uppertriangularind(A::AbstractMatrix) - m, n = size(A) - I = Vector{Int}(undef, div(m * (m - 1), 2) + m * (n - m)) - offset = 0 - for i in 1:m - r = (i + 1):n - I[offset .- i .+ r] = i .+ m .* (r .- 1) - offset += length(r) - end - return I -end - -# SVD_pullback: pullback implementation for general (possibly truncated) SVD -# -# Arguments are U, S and Vd of full (non-truncated, but still thin) SVD, as well as -# cotangent ΔU, ΔS, ΔVd variables of truncated SVD -# -# Checks whether the cotangent variables are such that they would couple to gauge-dependent -# degrees of freedom (phases of singular vectors), and prints a warning if this is the case -# -# An implementation that only uses U, S, and Vd from truncated SVD is also possible, but -# requires solving a Sylvester equation, which does not seem to be supported on GPUs. -# -# Other implementation considerations for GPU compatibility: -# no scalar indexing, lots of broadcasting and views -# -function svd_pullback!(ΔA::AbstractMatrix, U::AbstractMatrix, S::AbstractVector, - Vd::AbstractMatrix, ΔU, ΔS, ΔVd; - tol::Real=default_pullback_gaugetol(S)) - - # Basic size checks and determination - m, n = size(U, 1), size(Vd, 2) - size(U, 2) == size(Vd, 1) == length(S) == min(m, n) || throw(DimensionMismatch()) - p = -1 - if !(ΔU isa AbstractZero) - m == size(ΔU, 1) || throw(DimensionMismatch()) - p = size(ΔU, 2) - end - if !(ΔVd isa AbstractZero) - n == size(ΔVd, 2) || throw(DimensionMismatch()) - if p == -1 - p = size(ΔVd, 1) - else - p == size(ΔVd, 1) || throw(DimensionMismatch()) - end - end - if !(ΔS isa AbstractZero) - if p == -1 - p = length(ΔS) - else - p == length(ΔS) || throw(DimensionMismatch()) - end - end - Up = view(U, :, 1:p) - Vp = view(Vd, 1:p, :)' - Sp = view(S, 1:p) - - # rank - r = searchsortedlast(S, tol; rev=true) - - # compute antihermitian part of projection of ΔU and ΔV onto U and V - # also already subtract this projection from ΔU and ΔV - if !(ΔU isa AbstractZero) - UΔU = Up' * ΔU - aUΔU = rmul!(UΔU - UΔU', 1 / 2) - if m > p - ΔU -= Up * UΔU - end - else - aUΔU = fill!(similar(U, (p, p)), 0) - end - if !(ΔVd isa AbstractZero) - VΔV = Vp' * ΔVd' - aVΔV = rmul!(VΔV - VΔV', 1 / 2) - if n > p - ΔVd -= VΔV' * Vp' - end - else - aVΔV = fill!(similar(Vd, (p, p)), 0) - end - - # check whether cotangents arise from gauge-invariance objective function - mask = abs.(Sp' .- Sp) .< tol - Δgauge = norm(view(aUΔU, mask) + view(aVΔV, mask), Inf) - if p > r - rprange = (r + 1):p - Δgauge = max(Δgauge, norm(view(aUΔU, rprange, rprange), Inf)) - Δgauge = max(Δgauge, norm(view(aVΔV, rprange, rprange), Inf)) - end - Δgauge < tol || - @warn "`svd` cotangents sensitive to gauge choice: (|Δgauge| = $Δgauge)" - - UdΔAV = (aUΔU .+ aVΔV) .* safe_inv.(Sp' .- Sp, tol) .+ - (aUΔU .- aVΔV) .* safe_inv.(Sp' .+ Sp, tol) - if !(ΔS isa ZeroTangent) - UdΔAV[diagind(UdΔAV)] .+= real.(ΔS) - # in principle, ΔS is real, but maybe not if coming from an anyonic tensor - end - mul!(ΔA, Up, UdΔAV * Vp') - - if r > p # contribution from truncation - Ur = view(U, :, (p + 1):r) - Vr = view(Vd, (p + 1):r, :)' - Sr = view(S, (p + 1):r) - - if !(ΔU isa AbstractZero) - UrΔU = Ur' * ΔU - if m > r - ΔU -= Ur * UrΔU # subtract this part from ΔU - end - else - UrΔU = fill!(similar(U, (r - p, p)), 0) - end - if !(ΔVd isa AbstractZero) - VrΔV = Vr' * ΔVd' - if n > r - ΔVd -= VrΔV' * Vr' # subtract this part from ΔV - end - else - VrΔV = fill!(similar(Vd, (r - p, p)), 0) - end - - X = (1 // 2) .* ((UrΔU .+ VrΔV) .* safe_inv.(Sp' .- Sr, tol) .+ - (UrΔU .- VrΔV) .* safe_inv.(Sp' .+ Sr, tol)) - Y = (1 // 2) .* ((UrΔU .+ VrΔV) .* safe_inv.(Sp' .- Sr, tol) .- - (UrΔU .- VrΔV) .* safe_inv.(Sp' .+ Sr, tol)) - - # ΔA += Ur * X * Vp' + Up * Y' * Vr' - mul!(ΔA, Ur, X * Vp', 1, 1) - mul!(ΔA, Up * Y', Vr', 1, 1) - end - - if m > max(r, p) && !(ΔU isa AbstractZero) # remaining ΔU is already orthogonal to U[:,1:max(p,r)] - # ΔA += (ΔU .* safe_inv.(Sp', tol)) * Vp' - mul!(ΔA, ΔU .* safe_inv.(Sp', tol), Vp', 1, 1) - end - if n > max(r, p) && !(ΔVd isa AbstractZero) # remaining ΔV is already orthogonal to V[:,1:max(p,r)] - # ΔA += U * (safe_inv.(Sp, tol) .* ΔVd) - mul!(ΔA, Up, safe_inv.(Sp, tol) .* ΔVd, 1, 1) - end - return ΔA -end - -function eig_pullback!(ΔA::AbstractMatrix, D::AbstractVector, V::AbstractMatrix, ΔD, ΔV; - tol::Real=default_pullback_gaugetol(D)) - - # Basic size checks and determination - n = LinearAlgebra.checksquare(V) - n == length(D) || throw(DimensionMismatch()) - - if !(ΔV isa AbstractZero) - VdΔV = V' * ΔV - - mask = abs.(transpose(D) .- D) .< tol - Δgauge = norm(view(VdΔV, mask), Inf) - Δgauge < tol || - @warn "`eig` cotangents sensitive to gauge choice: (|Δgauge| = $Δgauge)" - - VdΔV .*= conj.(safe_inv.(transpose(D) .- D, tol)) - - if !(ΔD isa AbstractZero) - view(VdΔV, diagind(VdΔV)) .+= ΔD - end - PΔV = V' \ VdΔV - if eltype(ΔA) <: Real - ΔAc = mul!(VdΔV, PΔV, V') # recycle VdΔV memory - ΔA .= real.(ΔAc) - else - mul!(ΔA, PΔV, V') - end - else - PΔV = V' \ Diagonal(ΔD) - if eltype(ΔA) <: Real - ΔAc = PΔV * V' - ΔA .= real.(ΔAc) - else - mul!(ΔA, PΔV, V') - end - end - return ΔA -end - -function eigh_pullback!(ΔA::AbstractMatrix, D::AbstractVector, V::AbstractMatrix, ΔD, ΔV; - tol::Real=default_pullback_gaugetol(D)) - - # Basic size checks and determination - n = LinearAlgebra.checksquare(V) - n == length(D) || throw(DimensionMismatch()) - - if !(ΔV isa AbstractZero) - VdΔV = V' * ΔV - aVdΔV = rmul!(VdΔV - VdΔV', 1 / 2) - - mask = abs.(D' .- D) .< tol - Δgauge = norm(view(aVdΔV, mask)) - Δgauge < tol || - @warn "`eigh` cotangents sensitive to gauge choice: (|Δgauge| = $Δgauge)" - - aVdΔV .*= safe_inv.(D' .- D, tol) - - if !(ΔD isa AbstractZero) - view(aVdΔV, diagind(aVdΔV)) .+= real.(ΔD) - # in principle, ΔD is real, but maybe not if coming from an anyonic tensor - end - # recylce VdΔV space - mul!(ΔA, mul!(VdΔV, V, aVdΔV), V') - else - mul!(ΔA, V * Diagonal(ΔD), V') - end - return ΔA -end - -function qr_pullback!(ΔA::AbstractMatrix, Q::AbstractMatrix, R::AbstractMatrix, ΔQ, ΔR; - tol::Real=default_pullback_gaugetol(R)) - Rd = view(R, diagind(R)) - p = something(findlast(≥(tol) ∘ abs, Rd), 0) - m, n = size(R) - - Q1 = view(Q, :, 1:p) - R1 = view(R, 1:p, :) - R11 = view(R, 1:p, 1:p) - - ΔA1 = view(ΔA, :, 1:p) - ΔQ1 = view(ΔQ, :, 1:p) - ΔR1 = view(ΔR, 1:p, :) - - M = similar(R, (p, p)) - ΔR isa AbstractZero || mul!(M, ΔR1, R1') - ΔQ isa AbstractZero || mul!(M, Q1', ΔQ1, -1, !(ΔR isa AbstractZero)) - view(M, lowertriangularind(M)) .= conj.(view(M, uppertriangularind(M))) - if eltype(M) <: Complex - Md = view(M, diagind(M)) - Md .= real.(Md) - end - - ΔA1 .= ΔQ1 - mul!(ΔA1, Q1, M, +1, 1) - - if n > p - R12 = view(R, 1:p, (p + 1):n) - ΔA2 = view(ΔA, :, (p + 1):n) - ΔR12 = view(ΔR, 1:p, (p + 1):n) - - if ΔR isa AbstractZero - ΔA2 .= zero(eltype(ΔA)) - else - mul!(ΔA2, Q1, ΔR12) - mul!(ΔA1, ΔA2, R12', -1, 1) - end - end - if m > p && !(ΔQ isa AbstractZero) # case where R is not full rank - Q2 = view(Q, :, (p + 1):m) - ΔQ2 = view(ΔQ, :, (p + 1):m) - Q1dΔQ2 = Q1' * ΔQ2 - Δgauge = norm(mul!(copy(ΔQ2), Q1, Q1dΔQ2, -1, 1), Inf) - Δgauge < tol || - @warn "`qr` cotangents sensitive to gauge choice: (|Δgauge| = $Δgauge)" - mul!(ΔA1, Q2, Q1dΔQ2', -1, 1) - end - rdiv!(ΔA1, UpperTriangular(R11)') - return ΔA -end - -function lq_pullback!(ΔA::AbstractMatrix, L::AbstractMatrix, Q::AbstractMatrix, ΔL, ΔQ; - tol::Real=default_pullback_gaugetol(L)) - Ld = view(L, diagind(L)) - p = something(findlast(≥(tol) ∘ abs, Ld), 0) - m, n = size(L) - - L1 = view(L, :, 1:p) - L11 = view(L, 1:p, 1:p) - Q1 = view(Q, 1:p, :) - - ΔA1 = view(ΔA, 1:p, :) - ΔQ1 = view(ΔQ, 1:p, :) - ΔL1 = view(ΔL, :, 1:p) - - M = similar(L, (p, p)) - ΔL isa AbstractZero || mul!(M, L1', ΔL1) - ΔQ isa AbstractZero || mul!(M, ΔQ1, Q1', -1, !(ΔL isa AbstractZero)) - view(M, uppertriangularind(M)) .= conj.(view(M, lowertriangularind(M))) - if eltype(M) <: Complex - Md = view(M, diagind(M)) - Md .= real.(Md) - end - - ΔA1 .= ΔQ1 - mul!(ΔA1, M, Q1, +1, 1) - - if m > p - L21 = view(L, (p + 1):m, 1:p) - ΔA2 = view(ΔA, (p + 1):m, :) - ΔL21 = view(ΔL, (p + 1):m, 1:p) - - if ΔL isa AbstractZero - ΔA2 .= zero(eltype(ΔA)) - else - mul!(ΔA2, ΔL21, Q1) - mul!(ΔA1, L21', ΔA2, -1, 1) - end - end - if n > p && !(ΔQ isa AbstractZero) # case where R is not full rank - Q2 = view(Q, (p + 1):n, :) - ΔQ2 = view(ΔQ, (p + 1):n, :) - ΔQ2Q1d = ΔQ2 * Q1' - Δgauge = norm(mul!(copy(ΔQ2), ΔQ2Q1d, Q1, -1, 1)) - Δgauge < tol || - @warn "`lq` cotangents sensitive to gauge choice: (|Δgauge| = $Δgauge)" - mul!(ΔA1, ΔQ2Q1d', Q2, -1, 1) - end - ldiv!(LowerTriangular(L11)', ΔA1) - return ΔA -end - -function default_pullback_gaugetol(a) - n = norm(a, Inf) - return eps(eltype(n))^(3 / 4) * max(n, one(n)) + rightorth!_pullback(::NTuple{2,ZeroTangent}) = NoTangent(), ZeroTangent() + return LQ, rightorth!_pullback end diff --git a/src/TensorKit.jl b/src/TensorKit.jl index 6f7467e37..4301fe88b 100644 --- a/src/TensorKit.jl +++ b/src/TensorKit.jl @@ -31,7 +31,7 @@ export TruncationScheme export SpaceMismatch, SectorMismatch, IndexError # error types # general vector space methods -export space, field, dual, dim, reduceddim, dims, fuse, flip, isdual, oplus, +export space, field, dual, dim, reduceddim, dims, fuse, flip, isdual, oplus, ominus, insertleftunit, insertrightunit, removeunit # partial order for vector spaces @@ -47,7 +47,7 @@ export ZNSpace, SU2Irrep, U1Irrep, CU1Irrep # bendleft, bendright, foldleft, foldright, cycleclockwise, cycleanticlockwise # some unicode -export ⊕, ⊗, ×, ⊠, ℂ, ℝ, ℤ, ←, →, ≾, ≿, ≅, ≺, ≻ +export ⊕, ⊗, ⊖, ×, ⊠, ℂ, ℝ, ℤ, ←, →, ≾, ≿, ≅, ≺, ≻ export ℤ₂, ℤ₃, ℤ₄, U₁, SU, SU₂, CU₁ export fℤ₂, fU₁, fSU₂ export ℤ₂Space, ℤ₃Space, ℤ₄Space, U₁Space, CU₁Space, SU₂Space @@ -72,8 +72,9 @@ export inner, dot, norm, normalize, normalize!, tr export mul!, lmul!, rmul!, adjoint!, pinv, axpy!, axpby! export leftorth, rightorth, leftnull, rightnull, leftorth!, rightorth!, leftnull!, rightnull!, + left_polar, left_polar!, right_polar, right_polar!, tsvd!, tsvd, eigen, eigen!, eig, eig!, eigh, eigh!, exp, exp!, - isposdef, isposdef!, ishermitian, sylvester, rank, cond + isposdef, isposdef!, ishermitian, isisometry, isunitary, sylvester, rank, cond export braid, braid!, permute, permute!, transpose, transpose!, twist, twist!, repartition, repartition! export catdomain, catcodomain @@ -104,7 +105,11 @@ using TensorOperations: TensorOperations, @tensor, @tensoropt, @ncon, ncon using TensorOperations: IndexTuple, Index2Tuple, linearize, AbstractBackend const TO = TensorOperations +using MatrixAlgebraKit: MatrixAlgebraKit as MAK + using LRUCache +using OhMyThreads +using ScopedValues using TensorKitSectors import TensorKitSectors: dim, BraidingStyle, FusionStyle, ⊠, ⊗ @@ -117,7 +122,7 @@ using Base: @boundscheck, @propagate_inbounds, @constprop, SizeUnknown, HasLength, HasShape, IsInfinite, EltypeUnknown, HasEltype using Base.Iterators: product, filter -using LinearAlgebra: LinearAlgebra +using LinearAlgebra: LinearAlgebra, BlasFloat using LinearAlgebra: norm, dot, normalize, normalize!, tr, axpy!, axpby!, lmul!, rmul!, mul!, ldiv!, rdiv!, adjoint, adjoint!, transpose, transpose!, @@ -125,6 +130,7 @@ using LinearAlgebra: norm, dot, normalize, normalize!, tr, eigen, eigen!, svd, svd!, isposdef, isposdef!, ishermitian, rank, cond, Diagonal, Hermitian +using MatrixAlgebraKit import Base.Meta @@ -202,6 +208,7 @@ end #------------------------------------- # general definitions include("tensors/abstracttensor.jl") +include("tensors/backends.jl") include("tensors/blockiterator.jl") include("tensors/tensor.jl") include("tensors/adjoint.jl") @@ -211,10 +218,13 @@ include("tensors/tensoroperations.jl") include("tensors/treetransformers.jl") include("tensors/indexmanipulations.jl") include("tensors/diagonal.jl") -include("tensors/truncation.jl") -include("tensors/factorizations.jl") include("tensors/braidingtensor.jl") +include("tensors/factorizations/factorizations.jl") +using .Factorizations +# include("tensors/factorizations/matrixalgebrakit.jl") +# include("tensors/truncation.jl") + # # Planar macros and related functionality # #----------------------------------------- @nospecialize diff --git a/src/auxiliary/deprecate.jl b/src/auxiliary/deprecate.jl index fa7667b2b..b235cbd7c 100644 --- a/src/auxiliary/deprecate.jl +++ b/src/auxiliary/deprecate.jl @@ -1,29 +1,29 @@ import Base: transpose #! format: off -Base.@deprecate(permute(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; copy::Bool=false), - permute(t, (p1, p2); copy=copy)) -Base.@deprecate(transpose(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; copy::Bool=false), - transpose(t, (p1, p2); copy=copy)) -Base.@deprecate(braid(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple, levels; copy::Bool=false), - braid(t, (p1, p2), levels; copy=copy)) - -Base.@deprecate(tsvd(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...), - tsvd(t, (p₁, p₂); kwargs...)) -Base.@deprecate(leftorth(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...), - leftorth(t, (p₁, p₂); kwargs...)) -Base.@deprecate(rightorth(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...), - rightorth(t, (p₁, p₂); kwargs...)) -Base.@deprecate(leftnull(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...), - leftnull(t, (p₁, p₂); kwargs...)) -Base.@deprecate(rightnull(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...), - rightnull(t, (p₁, p₂); kwargs...)) -Base.@deprecate(LinearAlgebra.eigen(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...), - LinearAlgebra.eigen(t, (p₁, p₂); kwargs...), false) -Base.@deprecate(eig(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...), - eig(t, (p₁, p₂); kwargs...)) -Base.@deprecate(eigh(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...), - eigh(t, (p₁, p₂); kwargs...)) +# Base.@deprecate(permute(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; copy::Bool=false), +# permute(t, (p1, p2); copy=copy)) +# Base.@deprecate(transpose(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; copy::Bool=false), +# transpose(t, (p1, p2); copy=copy)) +# Base.@deprecate(braid(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple, levels; copy::Bool=false), +# braid(t, (p1, p2), levels; copy=copy)) + +# Base.@deprecate(tsvd(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...), +# tsvd(t, (p₁, p₂); kwargs...)) +# Base.@deprecate(leftorth(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...), +# leftorth(t, (p₁, p₂); kwargs...)) +# Base.@deprecate(rightorth(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...), +# rightorth(t, (p₁, p₂); kwargs...)) +# Base.@deprecate(leftnull(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...), +# leftnull(t, (p₁, p₂); kwargs...)) +# Base.@deprecate(rightnull(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...), +# rightnull(t, (p₁, p₂); kwargs...)) +# Base.@deprecate(LinearAlgebra.eigen(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...), +# LinearAlgebra.eigen(t, (p₁, p₂); kwargs...), false) +# Base.@deprecate(eig(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...), +# eig(t, (p₁, p₂); kwargs...)) +# Base.@deprecate(eigh(t::AbstractTensorMap, p₁::IndexTuple, p₂::IndexTuple; kwargs...), +# eigh(t, (p₁, p₂); kwargs...)) for f in (:rand, :randn, :zeros, :ones) @eval begin diff --git a/src/auxiliary/linalg.jl b/src/auxiliary/linalg.jl index 4277be67f..8a20eb8ac 100644 --- a/src/auxiliary/linalg.jl +++ b/src/auxiliary/linalg.jl @@ -46,341 +46,5 @@ Base.adjoint(alg::Union{SVD,SDD,Polar}) = alg const OFA = OrthogonalFactorizationAlgorithm const SVDAlg = Union{SVD,SDD} -# Matrix algebra: entrypoint for calling matrix methods from within tensor implementations -#------------------------------------------------------------------------------------------ -module MatrixAlgebra -# TODO: all methods tha twe define here will need an extended version for CuMatrix in the -# CUDA package extension. - -# TODO: other methods to include here: -# mul! (possibly call matmul! instead) -# adjoint! -# sylvester -# exp! -# schur!? -# - -using LinearAlgebra -using LinearAlgebra: BlasFloat, BlasReal, BlasComplex, checksquare - -using ..TensorKit: OrthogonalFactorizationAlgorithm, - QL, QLpos, QR, QRpos, LQ, LQpos, RQ, RQpos, SVD, SDD, Polar - -# only defined in >v1.7 -@static if VERSION < v"1.7-" - _rf_findmax((fm, im), (fx, ix)) = isless(fm, fx) ? (fx, ix) : (fm, im) - _argmax(f, domain) = mapfoldl(x -> (f(x), x), _rf_findmax, domain)[2] -else - _argmax(f, domain) = argmax(f, domain) -end - -# TODO: define for CuMatrix if we support this -function one!(A::StridedMatrix) - length(A) > 0 || return A - copyto!(A, LinearAlgebra.I) - return A -end - safesign(s::Real) = ifelse(s < zero(s), -one(s), +one(s)) safesign(s::Complex) = ifelse(iszero(s), one(s), s / abs(s)) - -function leftorth!(A::StridedMatrix{<:BlasFloat}, alg::Union{QR,QRpos}, atol::Real) - iszero(atol) || throw(ArgumentError("nonzero atol not supported by $alg")) - m, n = size(A) - k = min(m, n) - A, T = LAPACK.geqrt!(A, min(minimum(size(A)), 36)) - Q = similar(A, m, k) - for j in 1:k - for i in 1:m - Q[i, j] = i == j - end - end - Q = LAPACK.gemqrt!('L', 'N', A, T, Q) - R = triu!(A[1:k, :]) - - if isa(alg, QRpos) - @inbounds for j in 1:k - s = safesign(R[j, j]) - @simd for i in 1:m - Q[i, j] *= s - end - end - @inbounds for j in size(R, 2):-1:1 - for i in 1:min(k, j) - R[i, j] = R[i, j] * conj(safesign(R[i, i])) - end - end - end - return Q, R -end - -function leftorth!(A::StridedMatrix{<:BlasFloat}, alg::Union{QL,QLpos}, atol::Real) - iszero(atol) || throw(ArgumentError("nonzero atol not supported by $alg")) - m, n = size(A) - @assert m >= n - - nhalf = div(n, 2) - #swap columns in A - @inbounds for j in 1:nhalf, i in 1:m - A[i, j], A[i, n + 1 - j] = A[i, n + 1 - j], A[i, j] - end - Q, R = leftorth!(A, isa(alg, QL) ? QR() : QRpos(), atol) - - #swap columns in Q - @inbounds for j in 1:nhalf, i in 1:m - Q[i, j], Q[i, n + 1 - j] = Q[i, n + 1 - j], Q[i, j] - end - #swap rows and columns in R - @inbounds for j in 1:nhalf, i in 1:n - R[i, j], R[n + 1 - i, n + 1 - j] = R[n + 1 - i, n + 1 - j], R[i, j] - end - if isodd(n) - j = nhalf + 1 - @inbounds for i in 1:nhalf - R[i, j], R[n + 1 - i, j] = R[n + 1 - i, j], R[i, j] - end - end - return Q, R -end - -function leftorth!(A::StridedMatrix{<:BlasFloat}, alg::Union{SVD,SDD,Polar}, atol::Real) - U, S, V = alg isa SVD ? LAPACK.gesvd!('S', 'S', A) : LAPACK.gesdd!('S', A) - if isa(alg, Union{SVD,SDD}) - n = count(s -> s .> atol, S) - if n != length(S) - return U[:, 1:n], lmul!(Diagonal(S[1:n]), V[1:n, :]) - else - return U, lmul!(Diagonal(S), V) - end - else - iszero(atol) || throw(ArgumentError("nonzero atol not supported by $alg")) - # TODO: check Lapack to see if we can recycle memory of A - Q = mul!(A, U, V) - Sq = map!(sqrt, S, S) - SqV = lmul!(Diagonal(Sq), V) - R = SqV' * SqV - return Q, R - end -end - -function leftnull!(A::StridedMatrix{<:BlasFloat}, alg::Union{QR,QRpos}, atol::Real) - iszero(atol) || throw(ArgumentError("nonzero atol not supported by $alg")) - m, n = size(A) - m >= n || throw(ArgumentError("no null space if less rows than columns")) - - A, T = LAPACK.geqrt!(A, min(minimum(size(A)), 36)) - N = similar(A, m, max(0, m - n)) - fill!(N, 0) - for k in 1:(m - n) - N[n + k, k] = 1 - end - return N = LAPACK.gemqrt!('L', 'N', A, T, N) -end - -function leftnull!(A::StridedMatrix{<:BlasFloat}, alg::Union{SVD,SDD}, atol::Real) - size(A, 2) == 0 && return one!(similar(A, (size(A, 1), size(A, 1)))) - U, S, V = alg isa SVD ? LAPACK.gesvd!('A', 'N', A) : LAPACK.gesdd!('A', A) - indstart = count(>(atol), S) + 1 - return U[:, indstart:end] -end - -function rightorth!(A::StridedMatrix{<:BlasFloat}, alg::Union{LQ,LQpos,RQ,RQpos}, - atol::Real) - iszero(atol) || throw(ArgumentError("nonzero atol not supported by $alg")) - # TODO: geqrfp seems a bit slower than geqrt in the intermediate region around - # matrix size 100, which is the interesting region. => Investigate and fix - m, n = size(A) - k = min(m, n) - At = transpose!(similar(A, n, m), A) - - if isa(alg, RQ) || isa(alg, RQpos) - @assert m <= n - - mhalf = div(m, 2) - # swap columns in At - @inbounds for j in 1:mhalf, i in 1:n - At[i, j], At[i, m + 1 - j] = At[i, m + 1 - j], At[i, j] - end - Qt, Rt = leftorth!(At, isa(alg, RQ) ? QR() : QRpos(), atol) - - @inbounds for j in 1:mhalf, i in 1:n - Qt[i, j], Qt[i, m + 1 - j] = Qt[i, m + 1 - j], Qt[i, j] - end - @inbounds for j in 1:mhalf, i in 1:m - Rt[i, j], Rt[m + 1 - i, m + 1 - j] = Rt[m + 1 - i, m + 1 - j], Rt[i, j] - end - if isodd(m) - j = mhalf + 1 - @inbounds for i in 1:mhalf - Rt[i, j], Rt[m + 1 - i, j] = Rt[m + 1 - i, j], Rt[i, j] - end - end - Q = transpose!(A, Qt) - R = transpose!(similar(A, (m, m)), Rt) # TODO: efficient in place - return R, Q - else - Qt, Lt = leftorth!(At, alg', atol) - if m > n - L = transpose!(A, Lt) - Q = transpose!(similar(A, (n, n)), Qt) # TODO: efficient in place - else - Q = transpose!(A, Qt) - L = transpose!(similar(A, (m, m)), Lt) # TODO: efficient in place - end - return L, Q - end -end - -function rightorth!(A::StridedMatrix{<:BlasFloat}, alg::Union{SVD,SDD,Polar}, atol::Real) - U, S, V = alg isa SVD ? LAPACK.gesvd!('S', 'S', A) : LAPACK.gesdd!('S', A) - if isa(alg, Union{SVD,SDD}) - n = count(s -> s .> atol, S) - if n != length(S) - return rmul!(U[:, 1:n], Diagonal(S[1:n])), V[1:n, :] - else - return rmul!(U, Diagonal(S)), V - end - else - iszero(atol) || throw(ArgumentError("nonzero atol not supported by $alg")) - Q = mul!(A, U, V) - Sq = map!(sqrt, S, S) - USq = rmul!(U, Diagonal(Sq)) - L = USq * USq' - return L, Q - end -end - -function rightnull!(A::StridedMatrix{<:BlasFloat}, alg::Union{LQ,LQpos}, atol::Real) - iszero(atol) || throw(ArgumentError("nonzero atol not supported by $alg")) - m, n = size(A) - k = min(m, n) - At = adjoint!(similar(A, n, m), A) - At, T = LAPACK.geqrt!(At, min(k, 36)) - N = similar(A, max(n - m, 0), n) - fill!(N, 0) - for k in 1:(n - m) - N[k, m + k] = 1 - end - return N = LAPACK.gemqrt!('R', eltype(At) <: Real ? 'T' : 'C', At, T, N) -end - -function rightnull!(A::StridedMatrix{<:BlasFloat}, alg::Union{SVD,SDD}, atol::Real) - size(A, 1) == 0 && return one!(similar(A, (size(A, 2), size(A, 2)))) - U, S, V = alg isa SVD ? LAPACK.gesvd!('N', 'A', A) : LAPACK.gesdd!('A', A) - indstart = count(>(atol), S) + 1 - return V[indstart:end, :] -end - -function svd!(A::StridedMatrix{T}, alg::Union{SVD,SDD}) where {T<:BlasFloat} - # fix another type instability in LAPACK wrappers - TT = Tuple{Matrix{T},Vector{real(T)},Matrix{T}} - U, S, V = alg isa SVD ? LAPACK.gesvd!('S', 'S', A)::TT : LAPACK.gesdd!('S', A)::TT - return U, S, V -end - -function eig!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) where {T<:BlasReal} - n = checksquare(A) - n == 0 && return zeros(Complex{T}, 0), zeros(Complex{T}, 0, 0) - - A, DR, DI, VL, VR, _ = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : - (scale ? 'S' : 'N'), 'N', 'V', 'N', A) - D = complex.(DR, DI) - V = zeros(Complex{T}, n, n) - j = 1 - while j <= n - if DI[j] == 0 - vr = view(VR, :, j) - s = conj(sign(_argmax(abs, vr))) - V[:, j] .= s .* vr - else - vr = view(VR, :, j) - vi = view(VR, :, j + 1) - s = conj(sign(_argmax(abs, vr))) # vectors coming from lapack have already real absmax component - V[:, j] .= s .* (vr .+ im .* vi) - V[:, j + 1] .= s .* (vr .- im .* vi) - j += 1 - end - j += 1 - end - return D, V -end - -function eig!(A::StridedMatrix{T}; permute::Bool=true, - scale::Bool=true) where {T<:BlasComplex} - n = checksquare(A) - n == 0 && return zeros(T, 0), zeros(T, 0, 0) - D, V = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'V', 'N', - A)[[2, 4]] - for j in 1:n - v = view(V, :, j) - s = conj(sign(_argmax(abs, v))) - v .*= s - end - return D, V -end - -function eigh!(A::StridedMatrix{T}) where {T<:BlasFloat} - n = checksquare(A) - n == 0 && return zeros(real(T), 0), zeros(T, 0, 0) - D, V = LAPACK.syevr!('V', 'A', 'U', A, 0.0, 0.0, 0, 0, -1.0) - for j in 1:n - v = view(V, :, j) - s = conj(sign(_argmax(abs, v))) - v .*= s - end - return D, V -end - -## Old stuff and experiments - -# using LinearAlgebra: BlasFloat, Char, BlasInt, LAPACK, LAPACKException, -# DimensionMismatch, SingularException, PosDefException, chkstride1, -# checksquare, -# triu! - -# TODO: reconsider the following implementation -# Unfortunately, geqrfp seems a bit slower than geqrt in the intermediate region -# around matrix size 100, which is the interesting region. => Investigate and maybe fix -# function _leftorth!(A::StridedMatrix{<:BlasFloat}) -# m, n = size(A) -# A, τ = geqrfp!(A) -# Q = LAPACK.ormqr!('L', 'N', A, τ, eye(eltype(A), m, min(m, n))) -# R = triu!(A[1:min(m, n), :]) -# return Q, R -# end - -# geqrfp!: computes qrpos factorization, missing in Base -# geqrfp!(A::StridedMatrix{<:BlasFloat}) = -# ((m, n) = size(A); geqrfp!(A, similar(A, min(m, n)))) -# -# for (geqrfp, elty, relty) in -# ((:dgeqrfp_, :Float64, :Float64), (:sgeqrfp_, :Float32, :Float32), -# (:zgeqrfp_, :ComplexF64, :Float64), (:cgeqrfp_, :ComplexF32, :Float32)) -# @eval begin -# function geqrfp!(A::StridedMatrix{$elty}, tau::StridedVector{$elty}) -# chkstride1(A, tau) -# m, n = size(A) -# if length(tau) != min(m, n) -# throw(DimensionMismatch("tau has length $(length(tau)), but needs length $(min(m, n))")) -# end -# work = Vector{$elty}(1) -# lwork = BlasInt(-1) -# info = Ref{BlasInt}() -# for i = 1:2 # first call returns lwork as work[1] -# ccall((@blasfunc($geqrfp), liblapack), Nothing, -# (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, -# Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), -# Ref(m), Ref(n), A, Ref(max(1, stride(A, 2))), -# tau, work, Ref(lwork), info) -# chklapackerror(info[]) -# if i == 1 -# lwork = BlasInt(real(work[1])) -# resize!(work, lwork) -# end -# end -# A, tau -# end -# end -# end - -end diff --git a/src/auxiliary/random.jl b/src/auxiliary/random.jl index bc9df6f65..3289cdc3c 100644 --- a/src/auxiliary/random.jl +++ b/src/auxiliary/random.jl @@ -20,6 +20,6 @@ function randisometry!(rng::Random.AbstractRNG, A::AbstractMatrix) dims = size(A) dims[1] >= dims[2] || throw(DimensionMismatch("cannot create isometric matrix with dimensions $dims; isometry needs to be tall or square")) - Q, = MatrixAlgebra.leftorth!(Random.randn!(rng, A), QRpos(), 0) + Q, = leftorth!(Random.randn!(rng, A); alg=QRpos()) return copy!(A, Q) end diff --git a/src/spaces/cartesianspace.jl b/src/spaces/cartesianspace.jl index f29ed263d..fe12f3dc6 100644 --- a/src/spaces/cartesianspace.jl +++ b/src/spaces/cartesianspace.jl @@ -49,7 +49,13 @@ sectortype(::Type{CartesianSpace}) = Trivial Base.oneunit(::Type{CartesianSpace}) = CartesianSpace(1) Base.zero(::Type{CartesianSpace}) = CartesianSpace(0) + ⊕(V₁::CartesianSpace, V₂::CartesianSpace) = CartesianSpace(V₁.d + V₂.d) +function ⊖(V::CartesianSpace, W::CartesianSpace) + V ≿ W || throw(ArgumentError("$(W) is not a subspace of $(V)")) + return CartesianSpace(dim(V) - dim(W)) +end + fuse(V₁::CartesianSpace, V₂::CartesianSpace) = CartesianSpace(V₁.d * V₂.d) flip(V::CartesianSpace) = V diff --git a/src/spaces/complexspace.jl b/src/spaces/complexspace.jl index ff05888b8..1031db614 100644 --- a/src/spaces/complexspace.jl +++ b/src/spaces/complexspace.jl @@ -50,11 +50,18 @@ Base.conj(V::ComplexSpace) = ComplexSpace(dim(V), !isdual(V)) Base.oneunit(::Type{ComplexSpace}) = ComplexSpace(1) Base.zero(::Type{ComplexSpace}) = ComplexSpace(0) + function ⊕(V₁::ComplexSpace, V₂::ComplexSpace) return isdual(V₁) == isdual(V₂) ? ComplexSpace(dim(V₁) + dim(V₂), isdual(V₁)) : throw(SpaceMismatch("Direct sum of a vector space and its dual does not exist")) end +function ⊖(V::ComplexSpace, W::ComplexSpace) + (V ≿ W && isdual(V) == isdual(W)) || + throw(ArgumentError("$(W) is not a subspace of $(V)")) + return ComplexSpace(dim(V) - dim(W), isdual(V)) +end + fuse(V₁::ComplexSpace, V₂::ComplexSpace) = ComplexSpace(V₁.d * V₂.d) flip(V::ComplexSpace) = dual(V) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index 7aa746076..4d8389dcd 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -149,6 +149,13 @@ function ⊕(V₁::GradedSpace{I}, V₂::GradedSpace{I}) where {I<:Sector} return typeof(V₁)(dims; dual=dual1) end +function ⊖(V::GradedSpace{I}, W::GradedSpace{I}) where {I<:Sector} + dual = isdual(V) + V ≿ W && dual == isdual(W) || + throw(SpaceMismatch("$(W) is not a subspace of $(V)")) + return typeof(V)(c => dim(V, c) - dim(W, c) for c in sectors(V); dual) +end + function flip(V::GradedSpace{I}) where {I<:Sector} if isdual(V) typeof(V)(c => dim(V, c) for c in sectors(V)) diff --git a/src/spaces/vectorspaces.jl b/src/spaces/vectorspaces.jl index 844da081f..847182536 100644 --- a/src/spaces/vectorspaces.jl +++ b/src/spaces/vectorspaces.jl @@ -150,6 +150,17 @@ function ⊕ end ⊕(V::Vararg{ElementarySpace}) = foldl(⊕, V) const oplus = ⊕ +""" + ⊖(V::ElementarySpace, W::ElementarySpace) -> X::ElementarySpace + ominus(V::ElementarySpace, W::ElementarySpace) -> X::ElementarySpace + +Return a space that is equivalent to the orthogonal complement of `W` in `V`, +i.e. an instance `X::ElementarySpace` such that `V = W ⊕ X`. +""" +⊖(V₁::S, V₂::S) where {S<:ElementarySpace} +⊖(V₁::VectorSpace, V₂::VectorSpace) = ⊖(promote(V₁, V₂)...) +const ominus = ⊖ + """ ⊗(V₁::S, V₂::S, V₃::S...) where {S<:ElementarySpace} -> S diff --git a/src/tensors/backends.jl b/src/tensors/backends.jl new file mode 100644 index 000000000..0fc8f99f6 --- /dev/null +++ b/src/tensors/backends.jl @@ -0,0 +1,32 @@ +# Scheduler implementation +# ------------------------ +function select_scheduler(scheduler=OhMyThreads.Implementation.NotGiven(); kwargs...) + return if scheduler == OhMyThreads.Implementation.NotGiven() && isempty(kwargs) + Threads.nthreads() == 1 ? SerialScheduler() : DynamicScheduler() + else + OhMyThreads.Implementation._scheduler_from_userinput(scheduler; kwargs...) + end +end + +""" + const blockscheduler = ScopedValue{Scheduler}(SerialScheduler()) + +The default scheduler used when looping over different blocks in the matrix representation of a +tensor. +For controlling this value, see also [`set_blockscheduler`](@ref) and [`with_blockscheduler`](@ref). +""" +const blockscheduler = ScopedValue{Scheduler}(SerialScheduler()) + +""" + with_blockscheduler(f, [scheduler]; kwargs...) + +Run `f` in a scope where the `blockscheduler` is determined by `scheduler' and `kwargs...`. +""" +@inline function with_blockscheduler(f, scheduler=OhMyThreads.Implementation.NotGiven(); + kwargs...) + @with blockscheduler => select_scheduler(scheduler; kwargs...) f() +end + +# TODO: disable for trivial symmetry or small tensors? +default_blockscheduler(t::AbstractTensorMap) = default_blockscheduler(typeof(t)) +default_blockscheduler(::Type{T}) where {T<:AbstractTensorMap} = blockscheduler[] diff --git a/src/tensors/blockiterator.jl b/src/tensors/blockiterator.jl index 7929b5d19..2938de824 100644 --- a/src/tensors/blockiterator.jl +++ b/src/tensors/blockiterator.jl @@ -12,3 +12,34 @@ Base.IteratorSize(::BlockIterator) = Base.HasLength() Base.IteratorEltype(::BlockIterator) = Base.HasEltype() Base.eltype(::Type{<:BlockIterator{T}}) where {T} = Pair{sectortype(T),blocktype(T)} Base.length(iter::BlockIterator) = length(iter.structure) + +# TODO: fast-path when structures are the same? +# TODO: implement scheduler +""" + foreachblock(f, ts::AbstractTensorMap...; [scheduler]) + +Apply `f` to each block of `t` and the corresponding blocks of `ts`. +Optionally, `scheduler` can be used to parallelize the computation. +This function is equivalent to the following loop: + +```julia +for c in union(blocksectors.(ts)...) + bs = map(t -> block(t, c), ts) + f(c, bs) +end +``` +""" +function foreachblock(f, t::AbstractTensorMap, ts::AbstractTensorMap...; scheduler=nothing) + tensors = (t, ts...) + allsectors = union(blocksectors.(tensors)...) + foreach(allsectors) do c + return f(c, block.(tensors, Ref(c))) + end + return nothing +end +function foreachblock(f, t::AbstractTensorMap; scheduler=nothing) + foreach(blocks(t)) do (c, b) + return f(c, (b,)) + end + return nothing +end diff --git a/src/tensors/diagonal.jl b/src/tensors/diagonal.jl index 88d0c3b25..5a3840f1b 100644 --- a/src/tensors/diagonal.jl +++ b/src/tensors/diagonal.jl @@ -317,65 +317,6 @@ function LinearAlgebra.isposdef(d::DiagonalTensorMap) return all(isposdef, d.data) end -function eig!(d::DiagonalTensorMap) - return d, one(d) -end -function eigh!(d::DiagonalTensorMap{<:Real}) - return d, one(d) -end -function eigh!(d::DiagonalTensorMap{<:Complex}) - # TODO: should this test for hermiticity? `eigh!(::TensorMap)` also does not do this. - return DiagonalTensorMap(real(d.data), d.domain), one(d) -end - -function leftorth!(d::DiagonalTensorMap; alg=QR(), kwargs...) - @assert alg isa Union{QR,QL} - return one(d), d # TODO: this is only correct for `alg = QR()` or `alg = QL()` -end -function rightorth!(d::DiagonalTensorMap; alg=LQ(), kwargs...) - @assert alg isa Union{LQ,RQ} - return d, one(d) # TODO: this is only correct for `alg = LQ()` or `alg = RQ()` -end -# not much to do here: -leftnull!(d::DiagonalTensorMap; kwargs...) = leftnull!(TensorMap(d); kwargs...) -rightnull!(d::DiagonalTensorMap; kwargs...) = rightnull!(TensorMap(d); kwargs...) - -function tsvd!(d::DiagonalTensorMap; trunc=NoTruncation(), p::Real=2, alg=SDD()) - return _tsvd!(d, alg, trunc, p) -end -# helper function -function _compute_svddata!(d::DiagonalTensorMap, alg::Union{SVD,SDD}) - InnerProductStyle(d) === EuclideanInnerProduct() || throw_invalid_innerproduct(:tsvd!) - I = sectortype(d) - dims = SectorDict{I,Int}() - generator = Base.Iterators.map(blocks(d)) do (c, b) - lb = length(b.diag) - U = zerovector!(similar(b.diag, lb, lb)) - V = zerovector!(similar(b.diag, lb, lb)) - p = sortperm(b.diag; by=abs, rev=true) - for (i, pi) in enumerate(p) - U[pi, i] = MatrixAlgebra.safesign(b.diag[pi]) - V[i, pi] = 1 - end - Σ = abs.(view(b.diag, p)) - dims[c] = lb - return c => (U, Σ, V) - end - SVDdata = SectorDict(generator) - return SVDdata, dims -end - -function LinearAlgebra.svdvals(d::DiagonalTensorMap) - return SectorDict(c => LinearAlgebra.svdvals(b) for (c, b) in blocks(d)) -end -function LinearAlgebra.eigvals(d::DiagonalTensorMap) - return SectorDict(c => LinearAlgebra.eigvals(b) for (c, b) in blocks(d)) -end - -function LinearAlgebra.cond(d::DiagonalTensorMap, p::Real=2) - return LinearAlgebra.cond(Diagonal(d.data), p) -end - # matrix functions for f in (:exp, :cos, :sin, :tan, :cot, :cosh, :sinh, :tanh, :coth, :atan, :acot, :asinh, :sqrt, diff --git a/src/tensors/factorizations.jl b/src/tensors/factorizations.jl deleted file mode 100644 index 2292b6c17..000000000 --- a/src/tensors/factorizations.jl +++ /dev/null @@ -1,685 +0,0 @@ -# Tensor factorization -#---------------------- -function factorisation_scalartype(t::AbstractTensorMap) - T = scalartype(t) - return promote_type(Float32, typeof(zero(T) / sqrt(abs2(one(T))))) -end -factorisation_scalartype(f, t) = factorisation_scalartype(t) - -function permutedcopy_oftype(t::AbstractTensorMap, T::Type{<:Number}, p::Index2Tuple) - return permute!(similar(t, T, permute(space(t), p)), t, p) -end -function copy_oftype(t::AbstractTensorMap, T::Type{<:Number}) - return copy!(similar(t, T), t) -end - -""" - tsvd(t::AbstractTensorMap, (leftind, rightind)::Index2Tuple; - trunc::TruncationScheme = notrunc(), p::Real = 2, alg::Union{SVD, SDD} = SDD()) - -> U, S, V, ϵ - -Compute the (possibly truncated) singular value decomposition such that -`norm(permute(t, (leftind, rightind)) - U * S * V) ≈ ϵ`, where `ϵ` thus represents the truncation error. - -If `leftind` and `rightind` are not specified, the current partition of left and right -indices of `t` is used. In that case, less memory is allocated if one allows the data in -`t` to be destroyed/overwritten, by using `tsvd!(t, trunc = notrunc(), p = 2)`. - -A truncation parameter `trunc` can be specified for the new internal dimension, in which -case a truncated singular value decomposition will be computed. Choices are: -* `notrunc()`: no truncation (default); -* `truncerr(η::Real)`: truncates such that the p-norm of the truncated singular values is - smaller than `η`; -* `truncdim(χ::Int)`: truncates such that the equivalent total dimension of the internal - vector space is no larger than `χ`; -* `truncspace(V)`: truncates such that the dimension of the internal vector space is - smaller than that of `V` in any sector. -* `truncbelow(η::Real)`: truncates such that every singular value is larger then `η` ; - -Truncation options can also be combined using `&`, i.e. `truncbelow(η) & truncdim(χ)` will -choose the truncation space such that every singular value is larger than `η`, and the -equivalent total dimension of the internal vector space is no larger than `χ`. - -The method `tsvd` also returns the truncation error `ϵ`, computed as the `p` norm of the -singular values that were truncated. - -THe keyword `alg` can be equal to `SVD()` or `SDD()`, corresponding to the underlying LAPACK -algorithm that computes the decomposition (`_gesvd` or `_gesdd`). - -Orthogonality requires `InnerProductStyle(t) <: HasInnerProduct`, and `tsvd(!)` -is currently only implemented for `InnerProductStyle(t) === EuclideanInnerProduct()`. -""" -function tsvd(t::AbstractTensorMap, p::Index2Tuple; kwargs...) - tcopy = permutedcopy_oftype(t, factorisation_scalartype(tsvd, t), p) - return tsvd!(tcopy; kwargs...) -end - -function LinearAlgebra.svdvals(t::AbstractTensorMap) - tcopy = copy_oftype(t, factorisation_scalartype(tsvd, t)) - return LinearAlgebra.svdvals!(tcopy) -end - -""" - leftorth(t::AbstractTensorMap, (leftind, rightind)::Index2Tuple; - alg::OrthogonalFactorizationAlgorithm = QRpos()) -> Q, R - -Create orthonormal basis `Q` for indices in `leftind`, and remainder `R` such that -`permute(t, (leftind, rightind)) = Q*R`. - -If `leftind` and `rightind` are not specified, the current partition of left and right -indices of `t` is used. In that case, less memory is allocated if one allows the data in `t` -to be destroyed/overwritten, by using `leftorth!(t, alg = QRpos())`. - -Different algorithms are available, namely `QR()`, `QRpos()`, `SVD()` and `Polar()`. `QR()` -and `QRpos()` use a standard QR decomposition, producing an upper triangular matrix `R`. -`Polar()` produces a Hermitian and positive semidefinite `R`. `QRpos()` corrects the -standard QR decomposition such that the diagonal elements of `R` are positive. Only -`QRpos()` and `Polar()` are unique (no residual freedom) so that they always return the same -result for the same input tensor `t`. - -Orthogonality requires `InnerProductStyle(t) <: HasInnerProduct`, and -`leftorth(!)` is currently only implemented for - `InnerProductStyle(t) === EuclideanInnerProduct()`. -""" -function leftorth(t::AbstractTensorMap, p::Index2Tuple; kwargs...) - tcopy = permutedcopy_oftype(t, factorisation_scalartype(leftorth, t), p) - return leftorth!(tcopy; kwargs...) -end - -""" - rightorth(t::AbstractTensorMap, (leftind, rightind)::Index2Tuple; - alg::OrthogonalFactorizationAlgorithm = LQpos()) -> L, Q - -Create orthonormal basis `Q` for indices in `rightind`, and remainder `L` such that -`permute(t, (leftind, rightind)) = L*Q`. - -If `leftind` and `rightind` are not specified, the current partition of left and right -indices of `t` is used. In that case, less memory is allocated if one allows the data in `t` -to be destroyed/overwritten, by using `rightorth!(t, alg = LQpos())`. - -Different algorithms are available, namely `LQ()`, `LQpos()`, `RQ()`, `RQpos()`, `SVD()` and -`Polar()`. `LQ()` and `LQpos()` produce a lower triangular matrix `L` and are computed using -a QR decomposition of the transpose. `RQ()` and `RQpos()` produce an upper triangular -remainder `L` and only works if the total left dimension is smaller than or equal to the -total right dimension. `LQpos()` and `RQpos()` add an additional correction such that the -diagonal elements of `L` are positive. `Polar()` produces a Hermitian and positive -semidefinite `L`. Only `LQpos()`, `RQpos()` and `Polar()` are unique (no residual freedom) -so that they always return the same result for the same input tensor `t`. - -Orthogonality requires `InnerProductStyle(t) <: HasInnerProduct`, and -`rightorth(!)` is currently only implemented for -`InnerProductStyle(t) === EuclideanInnerProduct()`. -""" -function rightorth(t::AbstractTensorMap, p::Index2Tuple; kwargs...) - tcopy = permutedcopy_oftype(t, factorisation_scalartype(rightorth, t), p) - return rightorth!(tcopy; kwargs...) -end - -""" - leftnull(t::AbstractTensor, (leftind, rightind)::Index2Tuple; - alg::OrthogonalFactorizationAlgorithm = QRpos()) -> N - -Create orthonormal basis for the orthogonal complement of the support of the indices in -`leftind`, such that `N' * permute(t, (leftind, rightind)) = 0`. - -If `leftind` and `rightind` are not specified, the current partition of left and right -indices of `t` is used. In that case, less memory is allocated if one allows the data in `t` -to be destroyed/overwritten, by using `leftnull!(t, alg = QRpos())`. - -Different algorithms are available, namely `QR()` (or equivalently, `QRpos()`), `SVD()` and -`SDD()`. The first assumes that the matrix is full rank and requires `iszero(atol)` and -`iszero(rtol)`. With `SVD()` and `SDD()`, `rightnull` will use the corresponding singular -value decomposition, and one can specify an absolute or relative tolerance for which -singular values are to be considered zero, where `max(atol, norm(t)*rtol)` is used as upper -bound. - -Orthogonality requires `InnerProductStyle(t) <: HasInnerProduct`, and -`leftnull(!)` is currently only implemented for -`InnerProductStyle(t) === EuclideanInnerProduct()`. -""" -function leftnull(t::AbstractTensorMap, p::Index2Tuple; kwargs...) - tcopy = permutedcopy_oftype(t, factorisation_scalartype(leftnull, t), p) - return leftnull!(tcopy; kwargs...) -end - -""" - rightnull(t::AbstractTensor, (leftind, rightind)::Index2Tuple; - alg::OrthogonalFactorizationAlgorithm = LQ(), - atol::Real = 0.0, - rtol::Real = eps(real(float(one(scalartype(t)))))*iszero(atol)) -> N - -Create orthonormal basis for the orthogonal complement of the support of the indices in -`rightind`, such that `permute(t, (leftind, rightind))*N' = 0`. - -If `leftind` and `rightind` are not specified, the current partition of left and right -indices of `t` is used. In that case, less memory is allocated if one allows the data in `t` -to be destroyed/overwritten, by using `rightnull!(t, alg = LQpos())`. - -Different algorithms are available, namely `LQ()` (or equivalently, `LQpos`), `SVD()` and -`SDD()`. The first assumes that the matrix is full rank and requires `iszero(atol)` and -`iszero(rtol)`. With `SVD()` and `SDD()`, `rightnull` will use the corresponding singular -value decomposition, and one can specify an absolute or relative tolerance for which -singular values are to be considered zero, where `max(atol, norm(t)*rtol)` is used as upper -bound. - -Orthogonality requires `InnerProductStyle(t) <: HasInnerProduct`, and -`rightnull(!)` is currently only implemented for -`InnerProductStyle(t) === EuclideanInnerProduct()`. -""" -function rightnull(t::AbstractTensorMap, p::Index2Tuple; kwargs...) - tcopy = permutedcopy_oftype(t, factorisation_scalartype(rightnull, t), p) - return rightnull!(tcopy; kwargs...) -end - -""" - eigen(t::AbstractTensor, (leftind, rightind)::Index2Tuple; kwargs...) -> D, V - -Compute eigenvalue factorization of tensor `t` as linear map from `rightind` to `leftind`. - -If `leftind` and `rightind` are not specified, the current partition of left and right -indices of `t` is used. In that case, less memory is allocated if one allows the data in `t` -to be destroyed/overwritten, by using `eigen!(t)`. Note that the permuted tensor on which -`eigen!` is called should have equal domain and codomain, as otherwise the eigenvalue -decomposition is meaningless and cannot satisfy -``` -permute(t, (leftind, rightind)) * V = V * D -``` - -Accepts the same keyword arguments `scale` and `permute` as `eigen` of dense -matrices. See the corresponding documentation for more information. - -See also `eig` and `eigh` -""" -function LinearAlgebra.eigen(t::AbstractTensorMap, p::Index2Tuple; kwargs...) - tcopy = permutedcopy_oftype(t, factorisation_scalartype(eigen, t), p) - return eigen!(tcopy; kwargs...) -end - -function LinearAlgebra.eigvals(t::AbstractTensorMap; kwargs...) - tcopy = copy_oftype(t, factorisation_scalartype(eigen, t)) - return LinearAlgebra.eigvals!(tcopy; kwargs...) -end - -""" - eig(t::AbstractTensor, (leftind, rightind)::Index2Tuple; kwargs...) -> D, V - -Compute eigenvalue factorization of tensor `t` as linear map from `rightind` to `leftind`. -The function `eig` assumes that the linear map is not hermitian and returns type stable -complex valued `D` and `V` tensors for both real and complex valued `t`. See `eigh` for -hermitian linear maps - -If `leftind` and `rightind` are not specified, the current partition of left and right -indices of `t` is used. In that case, less memory is allocated if one allows the data in -`t` to be destroyed/overwritten, by using `eig!(t)`. Note that the permuted tensor on -which `eig!` is called should have equal domain and codomain, as otherwise the eigenvalue -decomposition is meaningless and cannot satisfy -``` -permute(t, (leftind, rightind)) * V = V * D -``` - -Accepts the same keyword arguments `scale` and `permute` as `eigen` of dense -matrices. See the corresponding documentation for more information. - -See also `eigen` and `eigh`. -""" -function eig(t::AbstractTensorMap, p::Index2Tuple; kwargs...) - tcopy = permutedcopy_oftype(t, factorisation_scalartype(eig, t), p) - return eig!(tcopy; kwargs...) -end - -""" - eigh(t::AbstractTensorMap, (leftind, rightind)::Index2Tuple) -> D, V - -Compute eigenvalue factorization of tensor `t` as linear map from `rightind` to `leftind`. -The function `eigh` assumes that the linear map is hermitian and `D` and `V` tensors with -the same `scalartype` as `t`. See `eig` and `eigen` for non-hermitian tensors. Hermiticity -requires that the tensor acts on inner product spaces, and the current implementation -requires `InnerProductStyle(t) === EuclideanInnerProduct()`. - -If `leftind` and `rightind` are not specified, the current partition of left and right -indices of `t` is used. In that case, less memory is allocated if one allows the data in -`t` to be destroyed/overwritten, by using `eigh!(t)`. Note that the permuted tensor on -which `eigh!` is called should have equal domain and codomain, as otherwise the eigenvalue -decomposition is meaningless and cannot satisfy -``` -permute(t, (leftind, rightind)) * V = V * D -``` - -See also `eigen` and `eig`. -""" -function eigh(t::AbstractTensorMap, p::Index2Tuple; kwargs...) - tcopy = permutedcopy_oftype(t, factorisation_scalartype(eigh, t), p) - return eigh!(tcopy; kwargs...) -end - -""" - isposdef(t::AbstractTensor, (leftind, rightind)::Index2Tuple) -> ::Bool - -Test whether a tensor `t` is positive definite as linear map from `rightind` to `leftind`. - -If `leftind` and `rightind` are not specified, the current partition of left and right -indices of `t` is used. In that case, less memory is allocated if one allows the data in -`t` to be destroyed/overwritten, by using `isposdef!(t)`. Note that the permuted tensor on -which `isposdef!` is called should have equal domain and codomain, as otherwise it is -meaningless. -""" -function LinearAlgebra.isposdef(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple) - tcopy = permutedcopy_oftype(t, factorisation_scalartype(isposdef, t), p) - return isposdef!(tcopy) -end - -function tsvd(t::AbstractTensorMap; kwargs...) - tcopy = copy_oftype(t, float(scalartype(t))) - return tsvd!(tcopy; kwargs...) -end -function leftorth(t::AbstractTensorMap; alg::OFA=QRpos(), kwargs...) - tcopy = copy_oftype(t, float(scalartype(t))) - return leftorth!(tcopy; alg=alg, kwargs...) -end -function rightorth(t::AbstractTensorMap; alg::OFA=LQpos(), kwargs...) - tcopy = copy_oftype(t, float(scalartype(t))) - return rightorth!(tcopy; alg=alg, kwargs...) -end -function leftnull(t::AbstractTensorMap; alg::OFA=QR(), kwargs...) - tcopy = copy_oftype(t, float(scalartype(t))) - return leftnull!(tcopy; alg=alg, kwargs...) -end -function rightnull(t::AbstractTensorMap; alg::OFA=LQ(), kwargs...) - tcopy = copy_oftype(t, float(scalartype(t))) - return rightnull!(tcopy; alg=alg, kwargs...) -end -function LinearAlgebra.eigen(t::AbstractTensorMap; kwargs...) - tcopy = copy_oftype(t, float(scalartype(t))) - return eigen!(tcopy; kwargs...) -end -function eig(t::AbstractTensorMap; kwargs...) - tcopy = copy_oftype(t, float(scalartype(t))) - return eig!(tcopy; kwargs...) -end -function eigh(t::AbstractTensorMap; kwargs...) - tcopy = copy_oftype(t, float(scalartype(t))) - return eigh!(tcopy; kwargs...) -end -function LinearAlgebra.isposdef(t::AbstractTensorMap) - tcopy = copy_oftype(t, float(scalartype(t))) - return isposdef!(tcopy) -end - -# Orthogonal factorizations (mutation for recycling memory): -# only possible if scalar type is floating point -# only correct if Euclidean inner product -#------------------------------------------------------------------------------------------ -const RealOrComplexFloat = Union{AbstractFloat,Complex{<:AbstractFloat}} - -function leftorth!(t::TensorMap{<:RealOrComplexFloat}; - alg::Union{QR,QRpos,QL,QLpos,SVD,SDD,Polar}=QRpos(), - atol::Real=zero(float(real(scalartype(t)))), - rtol::Real=(alg ∉ (SVD(), SDD())) ? zero(float(real(scalartype(t)))) : - eps(real(float(one(scalartype(t))))) * iszero(atol)) - InnerProductStyle(t) === EuclideanInnerProduct() || - throw_invalid_innerproduct(:leftorth!) - if !iszero(rtol) - atol = max(atol, rtol * norm(t)) - end - I = sectortype(t) - dims = SectorDict{I,Int}() - - # compute QR factorization for each block - if !isempty(blocks(t)) - generator = Base.Iterators.map(blocks(t)) do (c, b) - Qc, Rc = MatrixAlgebra.leftorth!(b, alg, atol) - dims[c] = size(Qc, 2) - return c => (Qc, Rc) - end - QRdata = SectorDict(generator) - end - - # construct new space - S = spacetype(t) - V = S(dims) - if alg isa Polar - @assert V ≅ domain(t) - W = domain(t) - elseif length(domain(t)) == 1 && domain(t) ≅ V - W = domain(t) - elseif length(codomain(t)) == 1 && codomain(t) ≅ V - W = codomain(t) - else - W = ProductSpace(V) - end - - # construct output tensors - T = float(scalartype(t)) - Q = similar(t, T, codomain(t) ← W) - R = similar(t, T, W ← domain(t)) - if !isempty(blocks(t)) - for (c, (Qc, Rc)) in QRdata - copy!(block(Q, c), Qc) - copy!(block(R, c), Rc) - end - end - return Q, R -end - -function leftnull!(t::TensorMap{<:RealOrComplexFloat}; - alg::Union{QR,QRpos,SVD,SDD}=QRpos(), - atol::Real=zero(float(real(scalartype(t)))), - rtol::Real=(alg ∉ (SVD(), SDD())) ? zero(float(real(scalartype(t)))) : - eps(real(float(one(scalartype(t))))) * iszero(atol)) - InnerProductStyle(t) === EuclideanInnerProduct() || - throw_invalid_innerproduct(:leftnull!) - if !iszero(rtol) - atol = max(atol, rtol * norm(t)) - end - I = sectortype(t) - dims = SectorDict{I,Int}() - - # compute QR factorization for each block - V = codomain(t) - if !isempty(blocksectors(V)) - generator = Base.Iterators.map(blocksectors(V)) do c - Nc = MatrixAlgebra.leftnull!(block(t, c), alg, atol) - dims[c] = size(Nc, 2) - return c => Nc - end - Ndata = SectorDict(generator) - end - - # construct new space - S = spacetype(t) - W = S(dims) - - # construct output tensor - T = float(scalartype(t)) - N = similar(t, T, V ← W) - if !isempty(blocksectors(V)) - for (c, Nc) in Ndata - copy!(block(N, c), Nc) - end - end - return N -end - -function rightorth!(t::TensorMap{<:RealOrComplexFloat}; - alg::Union{LQ,LQpos,RQ,RQpos,SVD,SDD,Polar}=LQpos(), - atol::Real=zero(float(real(scalartype(t)))), - rtol::Real=(alg ∉ (SVD(), SDD())) ? zero(float(real(scalartype(t)))) : - eps(real(float(one(scalartype(t))))) * iszero(atol)) - InnerProductStyle(t) === EuclideanInnerProduct() || - throw_invalid_innerproduct(:rightorth!) - if !iszero(rtol) - atol = max(atol, rtol * norm(t)) - end - I = sectortype(t) - dims = SectorDict{I,Int}() - - # compute LQ factorization for each block - if !isempty(blocks(t)) - generator = Base.Iterators.map(blocks(t)) do (c, b) - Lc, Qc = MatrixAlgebra.rightorth!(b, alg, atol) - dims[c] = size(Qc, 1) - return c => (Lc, Qc) - end - LQdata = SectorDict(generator) - end - - # construct new space - S = spacetype(t) - V = S(dims) - if alg isa Polar - @assert V ≅ codomain(t) - W = codomain(t) - elseif length(codomain(t)) == 1 && codomain(t) ≅ V - W = codomain(t) - elseif length(domain(t)) == 1 && domain(t) ≅ V - W = domain(t) - else - W = ProductSpace(V) - end - - # construct output tensors - T = float(scalartype(t)) - L = similar(t, T, codomain(t) ← W) - Q = similar(t, T, W ← domain(t)) - if !isempty(blocks(t)) - for (c, (Lc, Qc)) in LQdata - copy!(block(L, c), Lc) - copy!(block(Q, c), Qc) - end - end - return L, Q -end - -function rightnull!(t::TensorMap{<:RealOrComplexFloat}; - alg::Union{LQ,LQpos,SVD,SDD}=LQpos(), - atol::Real=zero(float(real(scalartype(t)))), - rtol::Real=(alg ∉ (SVD(), SDD())) ? zero(float(real(scalartype(t)))) : - eps(real(float(one(scalartype(t))))) * iszero(atol)) - InnerProductStyle(t) === EuclideanInnerProduct() || - throw_invalid_innerproduct(:rightnull!) - if !iszero(rtol) - atol = max(atol, rtol * norm(t)) - end - I = sectortype(t) - dims = SectorDict{I,Int}() - - # compute LQ factorization for each block - V = domain(t) - if !isempty(blocksectors(V)) - generator = Base.Iterators.map(blocksectors(V)) do c - Nc = MatrixAlgebra.rightnull!(block(t, c), alg, atol) - dims[c] = size(Nc, 1) - return c => Nc - end - Ndata = SectorDict(generator) - end - - # construct new space - S = spacetype(t) - W = S(dims) - - # construct output tensor - T = float(scalartype(t)) - N = similar(t, T, W ← V) - if !isempty(blocksectors(V)) - for (c, Nc) in Ndata - copy!(block(N, c), Nc) - end - end - return N -end - -function leftorth!(t::AdjointTensorMap; alg::OFA=QRpos()) - InnerProductStyle(t) === EuclideanInnerProduct() || - throw_invalid_innerproduct(:leftorth!) - return map(adjoint, reverse(rightorth!(adjoint(t); alg=alg'))) -end - -function rightorth!(t::AdjointTensorMap; alg::OFA=LQpos()) - InnerProductStyle(t) === EuclideanInnerProduct() || - throw_invalid_innerproduct(:rightorth!) - return map(adjoint, reverse(leftorth!(adjoint(t); alg=alg'))) -end - -function leftnull!(t::AdjointTensorMap; alg::OFA=QR(), kwargs...) - InnerProductStyle(t) === EuclideanInnerProduct() || - throw_invalid_innerproduct(:leftnull!) - return adjoint(rightnull!(adjoint(t); alg=alg', kwargs...)) -end - -function rightnull!(t::AdjointTensorMap; alg::OFA=LQ(), kwargs...) - InnerProductStyle(t) === EuclideanInnerProduct() || - throw_invalid_innerproduct(:rightnull!) - return adjoint(leftnull!(adjoint(t); alg=alg', kwargs...)) -end - -#------------------------------# -# Singular value decomposition # -#------------------------------# -function LinearAlgebra.svdvals!(t::TensorMap{<:RealOrComplexFloat}) - return SectorDict(c => LinearAlgebra.svdvals!(b) for (c, b) in blocks(t)) -end -LinearAlgebra.svdvals!(t::AdjointTensorMap) = svdvals!(adjoint(t)) - -function tsvd!(t::TensorMap{<:RealOrComplexFloat}; - trunc=NoTruncation(), p::Real=2, alg=SDD()) - return _tsvd!(t, alg, trunc, p) -end -function tsvd!(t::AdjointTensorMap; trunc=NoTruncation(), p::Real=2, alg=SDD()) - u, s, vt, err = tsvd!(adjoint(t); trunc=trunc, p=p, alg=alg) - return adjoint(vt), adjoint(s), adjoint(u), err -end - -# implementation dispatches on algorithm -function _tsvd!(t::TensorMap{<:RealOrComplexFloat}, alg::Union{SVD,SDD}, - trunc::TruncationScheme, p::Real=2) - # early return - if isempty(blocksectors(t)) - truncerr = zero(real(scalartype(t))) - return _empty_svdtensors(t)..., truncerr - end - - # compute SVD factorization for each block - S = spacetype(t) - SVDdata, dims = _compute_svddata!(t, alg) - Σdata = SectorDict(c => Σ for (c, (U, Σ, V)) in SVDdata) - truncdim = _compute_truncdim(Σdata, trunc, p) - truncerr = _compute_truncerr(Σdata, truncdim, p) - - # construct output tensors - U, Σ, V⁺ = _create_svdtensors(t, SVDdata, truncdim) - return U, Σ, V⁺, truncerr -end - -# helper functions -function _compute_svddata!(t::TensorMap, alg::Union{SVD,SDD}) - InnerProductStyle(t) === EuclideanInnerProduct() || throw_invalid_innerproduct(:tsvd!) - I = sectortype(t) - dims = SectorDict{I,Int}() - generator = Base.Iterators.map(blocks(t)) do (c, b) - U, Σ, V = MatrixAlgebra.svd!(b, alg) - dims[c] = length(Σ) - return c => (U, Σ, V) - end - SVDdata = SectorDict(generator) - return SVDdata, dims -end - -function _create_svdtensors(t::TensorMap{<:RealOrComplexFloat}, SVDdata, dims) - T = scalartype(t) - S = spacetype(t) - W = S(dims) - - Tr = real(T) - A = similarstoragetype(t, Tr) - Σ = DiagonalTensorMap{Tr,S,A}(undef, W) - - U = similar(t, codomain(t) ← W) - V⁺ = similar(t, W ← domain(t)) - for (c, (Uc, Σc, V⁺c)) in SVDdata - r = Base.OneTo(dims[c]) - copy!(block(U, c), view(Uc, :, r)) - copy!(block(Σ, c), Diagonal(view(Σc, r))) - copy!(block(V⁺, c), view(V⁺c, r, :)) - end - return U, Σ, V⁺ -end - -function _empty_svdtensors(t::TensorMap{<:RealOrComplexFloat}) - T = scalartype(t) - S = spacetype(t) - I = sectortype(t) - dims = SectorDict{I,Int}() - W = S(dims) - - Tr = real(T) - A = similarstoragetype(t, Tr) - Σ = DiagonalTensorMap{Tr,S,A}(undef, W) - - U = similar(t, codomain(t) ← W) - V⁺ = similar(t, W ← domain(t)) - return U, Σ, V⁺ -end - -#--------------------------# -# Eigenvalue decomposition # -#--------------------------# -function LinearAlgebra.eigen!(t::TensorMap{<:RealOrComplexFloat}) - return ishermitian(t) ? eigh!(t) : eig!(t) -end - -function LinearAlgebra.eigvals!(t::TensorMap{<:RealOrComplexFloat}; kwargs...) - return SectorDict(c => complex(LinearAlgebra.eigvals!(b; kwargs...)) - for (c, b) in blocks(t)) -end -function LinearAlgebra.eigvals!(t::AdjointTensorMap{<:RealOrComplexFloat}; kwargs...) - return SectorDict(c => conj!(complex(LinearAlgebra.eigvals!(b; kwargs...))) - for (c, b) in blocks(t)) -end - -function eigh!(t::TensorMap{<:RealOrComplexFloat}) - InnerProductStyle(t) === EuclideanInnerProduct() || throw_invalid_innerproduct(:eigh!) - domain(t) == codomain(t) || - throw(SpaceMismatch("`eigh!` requires domain and codomain to be the same")) - - T = scalartype(t) - I = sectortype(t) - S = spacetype(t) - dims = SectorDict{I,Int}(c => size(b, 1) for (c, b) in blocks(t)) - W = S(dims) - - Tr = real(T) - A = similarstoragetype(t, Tr) - D = DiagonalTensorMap{Tr,S,A}(undef, W) - V = similar(t, domain(t) ← W) - for (c, b) in blocks(t) - values, vectors = MatrixAlgebra.eigh!(b) - copy!(block(D, c), Diagonal(values)) - copy!(block(V, c), vectors) - end - return D, V -end - -function eig!(t::TensorMap{<:RealOrComplexFloat}; kwargs...) - domain(t) == codomain(t) || - throw(SpaceMismatch("`eig!` requires domain and codomain to be the same")) - - T = scalartype(t) - I = sectortype(t) - S = spacetype(t) - dims = SectorDict{I,Int}(c => size(b, 1) for (c, b) in blocks(t)) - W = S(dims) - - Tc = complex(T) - A = similarstoragetype(t, Tc) - D = DiagonalTensorMap{Tc,S,A}(undef, W) - V = similar(t, Tc, domain(t) ← W) - for (c, b) in blocks(t) - values, vectors = MatrixAlgebra.eig!(b; kwargs...) - copy!(block(D, c), Diagonal(values)) - copy!(block(V, c), vectors) - end - return D, V -end - -#--------------------------------------------------# -# Checks for hermiticity and positive definiteness # -#--------------------------------------------------# -function LinearAlgebra.ishermitian(t::TensorMap) - 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 - end - return true -end - -function LinearAlgebra.isposdef!(t::TensorMap) - 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 -end diff --git a/src/tensors/factorizations/adjoint.jl b/src/tensors/factorizations/adjoint.jl new file mode 100644 index 000000000..5f01d971b --- /dev/null +++ b/src/tensors/factorizations/adjoint.jl @@ -0,0 +1,78 @@ +# AdjointTensorMap +# ---------------- +# 1-arg functions +function initialize_output(::typeof(left_null!), t::AdjointTensorMap, + alg::AbstractAlgorithm) + return adjoint(initialize_output(right_null!, adjoint(t), alg)) +end +function initialize_output(::typeof(right_null!), t::AdjointTensorMap, + alg::AbstractAlgorithm) + return adjoint(initialize_output(left_null!, adjoint(t), alg)) +end + +function left_null!(t::AdjointTensorMap, N::AdjointTensorMap, alg::AbstractAlgorithm) + right_null!(adjoint(t), adjoint(N), alg) + return N +end +function right_null!(t::AdjointTensorMap, N::AdjointTensorMap, alg::AbstractAlgorithm) + left_null!(adjoint(t), adjoint(N), alg) + return N +end + +function MatrixAlgebraKit.is_left_isometry(t::AdjointTensorMap; kwargs...) + return is_right_isometry(adjoint(t); kwargs...) +end +function MatrixAlgebraKit.is_right_isometry(t::AdjointTensorMap; kwargs...) + return is_left_isometry(adjoint(t); kwargs...) +end + +# 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!)) + @eval function initialize_output(::typeof($left_f!), t::AdjointTensorMap, + alg::AbstractAlgorithm) + return reverse(adjoint.(initialize_output($right_f!, adjoint(t), alg))) + end + @eval function initialize_output(::typeof($right_f!), t::AdjointTensorMap, + alg::AbstractAlgorithm) + return reverse(adjoint.(initialize_output($left_f!, adjoint(t), alg))) + end + + @eval function $left_f!(t::AdjointTensorMap, + F::Tuple{AdjointTensorMap,AdjointTensorMap}, + alg::AbstractAlgorithm) + $right_f!(adjoint(t), reverse(adjoint.(F)), alg) + return F + end + @eval function $right_f!(t::AdjointTensorMap, + F::Tuple{AdjointTensorMap,AdjointTensorMap}, + alg::AbstractAlgorithm) + $left_f!(adjoint(t), reverse(adjoint.(F)), alg) + return F + end +end + +# 3-arg functions +for f! in (:svd_full!, :svd_compact!, :svd_trunc!) + @eval function initialize_output(::typeof($f!), t::AdjointTensorMap, + alg::AbstractAlgorithm) + return reverse(adjoint.(initialize_output($f!, adjoint(t), alg))) + end + _TS = f! === :svd_full! ? :AdjointTensorMap : DiagonalTensorMap + @eval function $f!(t::AdjointTensorMap, + F::Tuple{AdjointTensorMap,$_TS,AdjointTensorMap}, + alg::AbstractAlgorithm) + $f!(adjoint(t), reverse(adjoint.(F)), alg) + return F + end +end +# avoid amgiguity +function initialize_output(::typeof(svd_trunc!), t::AdjointTensorMap, + alg::TruncatedAlgorithm) + return initialize_output(svd_compact!, t, alg.alg) +end +# to fix ambiguity +function svd_trunc!(t::AdjointTensorMap, USVᴴ::Tuple{AdjointTensorMap,DiagonalTensorMap,AdjointTensorMap}, alg::TruncatedAlgorithm) + USVᴴ′ = svd_compact!(t, USVᴴ, alg.alg) + return truncate!(svd_trunc!, USVᴴ′, alg.trunc) +end diff --git a/src/tensors/factorizations/deprecations.jl b/src/tensors/factorizations/deprecations.jl new file mode 100644 index 000000000..8b1378917 --- /dev/null +++ b/src/tensors/factorizations/deprecations.jl @@ -0,0 +1 @@ + diff --git a/src/tensors/factorizations/factorizations.jl b/src/tensors/factorizations/factorizations.jl new file mode 100644 index 000000000..d6d3bef28 --- /dev/null +++ b/src/tensors/factorizations/factorizations.jl @@ -0,0 +1,169 @@ +# Tensor factorization +#---------------------- +# using submodule here to import MatrixAlgebraKit functions without polluting namespace +module Factorizations + +export eig, eig!, eigh, eigh! +export tsvd, tsvd!, svdvals, svdvals! +export leftorth, leftorth!, rightorth, rightorth! +export leftnull, leftnull!, rightnull, rightnull! +export copy_oftype, permutedcopy_oftype, one! +export TruncationScheme, notrunc, truncbelow, truncerr, truncdim, truncspace, PolarViaSVD +#export LAPACK_HouseholderQR, LAPACK_HouseholderLQ + +using ..TensorKit +using ..TensorKit: AdjointTensorMap, SectorDict, OFA, blocktype, foreachblock, one! + +using LinearAlgebra: LinearAlgebra, BlasFloat, Diagonal, svdvals, svdvals! +import LinearAlgebra: eigen, eigen!, isposdef, isposdef!, ishermitian + +using TensorOperations: Index2Tuple + +using MatrixAlgebraKit +using MatrixAlgebraKit: AbstractAlgorithm, TruncatedAlgorithm, TruncationStrategy, + NoTruncation, TruncationKeepAbove, TruncationKeepBelow, + TruncationIntersection, TruncationKeepFiltered, PolarViaSVD, + LAPACK_SVDAlgorithm, LAPACK_QRIteration, LAPACK_HouseholderQR, + LAPACK_HouseholderLQ +import MatrixAlgebraKit: default_algorithm, + copy_input, check_input, initialize_output, + qr_compact!, qr_full!, qr_null!, lq_compact!, lq_full!, lq_null!, + svd_compact!, svd_full!, svd_trunc!, svd_vals!, + eigh_full!, eigh_trunc!, eigh_vals!, + eig_full!, eig_trunc!, eig_vals!, + left_polar!, left_orth_polar!, right_polar!, right_orth_polar!, + left_null_svd!, right_null_svd!, + left_orth!, right_orth!, left_null!, right_null!, + truncate!, findtruncated, findtruncated_sorted, + diagview, isisometry + +include("utility.jl") +include("interface.jl") +include("implementations.jl") +include("matrixalgebrakit.jl") +include("truncation.jl") +include("deprecations.jl") +include("adjoint.jl") + +TensorKit.one!(A::AbstractMatrix) = MatrixAlgebraKit.one!(A) + +function isisometry(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple) + t = permute(t, (p₁, p₂); copy=false) + return isisometry(t) +end + +# Orthogonal factorizations (mutation for recycling memory): +# only possible if scalar type is floating point +# only correct if Euclidean inner product +#------------------------------------------------------------------------------------------ +const RealOrComplexFloat = Union{AbstractFloat,Complex{<:AbstractFloat}} + +# DiagonalTensorMap +# ----------------- +function leftorth!(d::DiagonalTensorMap; alg=QR(), kwargs...) + @assert alg isa Union{QR,QL} + return one(d), d # TODO: this is only correct for `alg = QR()` or `alg = QL()` +end +function rightorth!(d::DiagonalTensorMap; alg=LQ(), kwargs...) + @assert alg isa Union{LQ,RQ} + return d, one(d) # TODO: this is only correct for `alg = LQ()` or `alg = RQ()` +end +leftnull!(d::DiagonalTensorMap; kwargs...) = leftnull!(TensorMap(d); kwargs...) +rightnull!(d::DiagonalTensorMap; kwargs...) = rightnull!(TensorMap(d); kwargs...) + +function tsvd!(d::DiagonalTensorMap; trunc=NoTruncation(), p::Real=2, alg=SDD()) + return _tsvd!(d, alg, trunc, p) +end + +# helper function +function _compute_svddata!(d::DiagonalTensorMap, alg::Union{SVD,SDD}) + InnerProductStyle(d) === EuclideanInnerProduct() || throw_invalid_innerproduct(:tsvd!) + I = sectortype(d) + dims = SectorDict{I,Int}() + generator = Base.Iterators.map(blocks(d)) do (c, b) + lb = length(b.diag) + U = zerovector!(similar(b.diag, lb, lb)) + V = zerovector!(similar(b.diag, lb, lb)) + p = sortperm(b.diag; by=abs, rev=true) + for (i, pi) in enumerate(p) + U[pi, i] = safesign(b.diag[pi]) + V[i, pi] = 1 + end + Σ = abs.(view(b.diag, p)) + dims[c] = lb + return c => (U, Σ, V) + end + SVDdata = SectorDict(generator) + return SVDdata, dims +end + +eig!(d::DiagonalTensorMap) = d, one(d) +eigh!(d::DiagonalTensorMap{<:Real}) = d, one(d) +eigh!(d::DiagonalTensorMap{<:Complex}) = DiagonalTensorMap(real(d.data), d.domain), one(d) + +function LinearAlgebra.svdvals(d::DiagonalTensorMap) + return SectorDict(c => LinearAlgebra.svdvals(b) for (c, b) in blocks(d)) +end +function LinearAlgebra.eigvals(d::DiagonalTensorMap) + return SectorDict(c => LinearAlgebra.eigvals(b) for (c, b) in blocks(d)) +end + +function LinearAlgebra.cond(d::DiagonalTensorMap, p::Real=2) + return LinearAlgebra.cond(Diagonal(d.data), p) +end +#------------------------------# +# Singular value decomposition # +#------------------------------# +function LinearAlgebra.svdvals!(t::TensorMap{<:RealOrComplexFloat}) + return SectorDict(c => LinearAlgebra.svdvals!(b) for (c, b) in blocks(t)) +end +LinearAlgebra.svdvals!(t::AdjointTensorMap) = svdvals!(adjoint(t)) + +#--------------------------# +# Eigenvalue decomposition # +#--------------------------# + +function LinearAlgebra.eigvals!(t::TensorMap{<:RealOrComplexFloat}; kwargs...) + return SectorDict(c => complex(LinearAlgebra.eigvals!(b; kwargs...)) + for (c, b) in blocks(t)) +end +function LinearAlgebra.eigvals!(t::AdjointTensorMap{<:RealOrComplexFloat}; kwargs...) + return SectorDict(c => conj!(complex(LinearAlgebra.eigvals!(b; kwargs...))) + for (c, b) in blocks(t)) +end + +#--------------------------------------------------# +# Checks for hermiticity and positive definiteness # +#--------------------------------------------------# +function LinearAlgebra.ishermitian(t::TensorMap) + 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 + end + return true +end + +function LinearAlgebra.isposdef!(t::TensorMap) + 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 +end + +# TODO: tolerances are per-block, not global or weighted - does that matter? +function MatrixAlgebraKit.is_left_isometry(t::AbstractTensorMap; kwargs...) + domain(t) ≾ codomain(t) || return false + f((c, b)) = MatrixAlgebraKit.is_left_isometry(b; kwargs...) + return all(f, blocks(t)) +end +function MatrixAlgebraKit.is_right_isometry(t::AbstractTensorMap; kwargs...) + domain(t) ≿ codomain(t) || return false + f((c, b)) = MatrixAlgebraKit.is_right_isometry(b; kwargs...) + return all(f, blocks(t)) +end + +end diff --git a/src/tensors/factorizations/implementations.jl b/src/tensors/factorizations/implementations.jl new file mode 100644 index 000000000..49a63b9c1 --- /dev/null +++ b/src/tensors/factorizations/implementations.jl @@ -0,0 +1,182 @@ +_kindof(::Union{SVD,SDD}) = :svd +_kindof(::Union{QR,QRpos}) = :qr +_kindof(::Union{LQ,LQpos}) = :lq +_kindof(::Polar) = :polar + +_kindof(::LAPACK_HouseholderQR) = :qr +_kindof(::LAPACK_HouseholderLQ) = :lq +_kindof(::LAPACK_SVDAlgorithm) = :svd +_kindof(::PolarViaSVD) = :polar + +leftorth!(t; alg=nothing, kwargs...) = _leftorth!(t, alg; kwargs...) + +function _leftorth!(t::AbstractTensorMap, alg::Nothing, ; kwargs...) + return isempty(kwargs) ? left_orth!(t) : left_orth!(t; trunc=(; kwargs...)) +end +function _leftorth!(t::AbstractTensorMap, alg::Union{QL,QLpos}; kwargs...) + trunc = isempty(kwargs) ? nothing : (; kwargs...) + + if alg == QL() || alg == QLpos() + _reverse!(t; dims=2) + Q, R = left_orth!(t; kind=:qr, alg_qr=(; positive=alg == QLpos()), trunc) + _reverse!(Q; dims=2) + _reverse!(R) + return Q, R + end +end +function _leftorth!(t, alg::Union{OFA,AbstractAlgorithm}; kwargs...) + trunc = isempty(kwargs) ? nothing : (; kwargs...) + + alg isa OFA && Base.depwarn(lazy"$alg is deprecated", :leftorth!) + + kind = _kindof(alg) + if kind == :svd + alg_svd = alg === LAPACK_QRIteration() ? alg : + alg === LAPACK_DivideAndConquer() ? alg : + alg === SVD() ? LAPACK_QRIteration() : + alg === SDD() ? LAPACK_DivideAndConquer() : + throw(ArgumentError(lazy"Unknown algorithm $alg")) + return left_orth!(t; kind, alg_svd, trunc) + elseif kind == :qr + alg_qr = (; positive=(alg == QRpos())) + return left_orth!(t; kind, alg_qr, trunc) + elseif kind == :polar + return left_orth!(t; kind, trunc) + else + throw(ArgumentError(lazy"Invalid algorithm: $alg")) + end +end +# fallback to MatrixAlgebraKit version +_leftorth!(t, alg; kwargs...) = left_orth!(t, alg; kwargs...) + +function leftnull!(t::AbstractTensorMap; + alg::Union{LAPACK_HouseholderQR,LAPACK_QRIteration, LAPACK_DivideAndConquer,PolarViaSVD,QR,QRpos,SVD,SDD,Nothing}=nothing, kwargs...) + InnerProductStyle(t) === EuclideanInnerProduct() || + throw_invalid_innerproduct(:leftnull!) + trunc = isempty(kwargs) ? nothing : (; kwargs...) + alg isa OFA && Base.depwarn(lazy"$alg is deprecated", :leftnull!) + + isnothing(alg) && return left_null!(t; trunc) + + kind = _kindof(alg) + if kind == :svd + alg_svd = alg === LAPACK_QRIteration() ? alg : + alg === LAPACK_DivideAndConquer() ? alg : + alg === SVD() ? LAPACK_QRIteration() : + alg === SDD() ? LAPACK_DivideAndConquer() : + throw(ArgumentError(lazy"Unknown algorithm $alg")) + return left_null!(t; kind, alg_svd, trunc) + elseif kind == :qr + alg_qr = (; positive=(alg == QRpos())) + return left_null!(t; kind, alg_qr, trunc) + else + throw(ArgumentError(lazy"Invalid `leftnull!` algorithm: $alg")) + end +end + +function rightorth!(t::AbstractTensorMap; + alg::Union{LAPACK_HouseholderLQ,LAPACK_QRIteration, LAPACK_DivideAndConquer,PolarViaSVD,LQ,LQpos,RQ,RQpos,SVD,SDD,Polar,Nothing}=nothing, kwargs...) + InnerProductStyle(t) === EuclideanInnerProduct() || + throw_invalid_innerproduct(:rightorth!) + trunc = isempty(kwargs) ? nothing : (; kwargs...) + + alg isa OFA && Base.depwarn(lazy"$alg is deprecated", :rightorth!) + + isnothing(alg) && return right_orth!(t; trunc) + + if alg == RQ() || alg == RQpos() + _reverse!(t; dims=1) + L, Q = right_orth!(t; kind=:lq, alg_lq=(; positive=alg == RQpos()), trunc) + _reverse!(Q; dims=1) + _reverse!(L) + return L, Q + end + + kind = _kindof(alg) + if kind == :svd + alg_svd = alg === LAPACK_QRIteration() ? alg : + alg === LAPACK_DivideAndConquer() ? alg : + alg === SVD() ? LAPACK_QRIteration() : + alg === SDD() ? LAPACK_DivideAndConquer() : + throw(ArgumentError(lazy"Unknown algorithm $alg")) + return right_orth!(t; kind, alg_svd, trunc) + elseif kind == :lq + alg_lq = (; positive=(alg == LQpos())) + return right_orth!(t; kind, alg_lq, trunc) + elseif kind == :polar + return right_orth!(t; kind, trunc) + else + throw(ArgumentError(lazy"Invalid `rightorth!` algorithm: $alg")) + end +end + +function rightnull!(t::AbstractTensorMap; + alg::Union{LAPACK_HouseholderLQ, LAPACK_QRIteration, LAPACK_DivideAndConquer,PolarViaSVD,LQ,LQpos,SVD,SDD,Nothing}=nothing, kwargs...) + InnerProductStyle(t) === EuclideanInnerProduct() || + throw_invalid_innerproduct(:rightnull!) + trunc = isempty(kwargs) ? nothing : (; kwargs...) + + alg isa OFA && Base.depwarn(lazy"$alg is deprecated", :rightnull!) + + isnothing(alg) && return right_null!(t; trunc) + + kind = _kindof(alg) + if kind == :svd + alg_svd = alg === LAPACK_QRIteration() ? alg : + alg === LAPACK_DivideAndConquer() ? alg : + alg === SVD() ? LAPACK_QRIteration() : + alg === SDD() ? LAPACK_DivideAndConquer() : + throw(ArgumentError(lazy"Unknown algorithm $alg")) + return right_null!(t; kind, alg_svd, trunc) + elseif kind == :lq + alg_lq = (; positive=(alg == LQpos())) + return right_null!(t; kind, alg_lq, trunc) + else + throw(ArgumentError(lazy"Invalid `rightnull!` algorithm: $alg")) + end +end + +# Eigenvalue decomposition +# ------------------------ +function eigh!(t::AbstractTensorMap; trunc=notrunc(), kwargs...) + InnerProductStyle(t) === EuclideanInnerProduct() || throw_invalid_innerproduct(:eigh!) + if trunc == notrunc() + return eigh_full!(t; kwargs...) + else + return eigh_trunc!(t; trunc, kwargs...) + end +end + +function eig!(t::AbstractTensorMap; trunc=notrunc(), kwargs...) + InnerProductStyle(t) === EuclideanInnerProduct() || throw_invalid_innerproduct(:eig!) + + if trunc == notrunc() + return eig_full!(t; kwargs...) + else + return eig_trunc!(t; trunc, kwargs...) + end +end + +function eigen!(t::AbstractTensorMap; kwargs...) + return ishermitian(t) ? eigh!(t; kwargs...) : eig!(t; kwargs...) +end + +# Singular value decomposition +# ---------------------------- +function tsvd!(t::AbstractTensorMap; trunc=notrunc(), p=nothing, alg=nothing, kwargs...) + InnerProductStyle(t) === EuclideanInnerProduct() || throw_invalid_innerproduct(:tsvd!) + isnothing(p) || Base.depwarn("p is no longer supported", :tsvd!) + + if alg isa OFA + Base.depwarn(lazy"$alg is deprecated", :tsvd!) + alg = alg === SVD() ? LAPACK_QRIteration() : + alg === SDD() ? LAPACK_DivideAndConquer() : + throw(ArgumentError(lazy"Unknown algorithm $alg")) + end + + if trunc == notrunc() + return svd_compact!(t; alg, kwargs...) + else + return svd_trunc!(t; trunc, alg, kwargs...) + end +end diff --git a/src/tensors/factorizations/interface.jl b/src/tensors/factorizations/interface.jl new file mode 100644 index 000000000..22c3c75d5 --- /dev/null +++ b/src/tensors/factorizations/interface.jl @@ -0,0 +1,241 @@ +@doc """ + tsvd(t::AbstractTensorMap, [(leftind, rightind)::Index2Tuple]; + trunc::TruncationScheme = notrunc(), p::Real = 2, alg::Union{SVD, SDD} = SDD()) + -> U, S, V, ϵ + tsvd!(t::AbstractTensorMap, trunc::TruncationScheme = notrunc(), p::Real = 2, alg::Union{SVD, SDD} = SDD()) + -> U, S, V, ϵ + +Compute the (possibly truncated) singular value decomposition such that +`norm(permute(t, (leftind, rightind)) - U * S * V) ≈ ϵ`, where `ϵ` thus represents the truncation error. + +If `leftind` and `rightind` are not specified, the current partition of left and right +indices of `t` is used. In that case, less memory is allocated if one allows the data in +`t` to be destroyed/overwritten, by using `tsvd!(t, trunc = notrunc(), p = 2)`. + +A truncation parameter `trunc` can be specified for the new internal dimension, in which +case a truncated singular value decomposition will be computed. Choices are: +* `notrunc()`: no truncation (default); +* `truncerr(η::Real)`: truncates such that the p-norm of the truncated singular values is + smaller than `η`; +* `truncdim(χ::Int)`: truncates such that the equivalent total dimension of the internal + vector space is no larger than `χ`; +* `truncspace(V)`: truncates such that the dimension of the internal vector space is + smaller than that of `V` in any sector. +* `truncbelow(η::Real)`: truncates such that every singular value is larger then `η` ; + +Truncation options can also be combined using `&`, i.e. `truncbelow(η) & truncdim(χ)` will +choose the truncation space such that every singular value is larger than `η`, and the +equivalent total dimension of the internal vector space is no larger than `χ`. + +The method `tsvd` also returns the truncation error `ϵ`, computed as the `p` norm of the +singular values that were truncated. + +The keyword `alg` can be equal to `SVD()` or `SDD()`, corresponding to the underlying LAPACK +algorithm that computes the decomposition (`_gesvd` or `_gesdd`). + +Orthogonality requires `InnerProductStyle(t) <: HasInnerProduct`, and `tsvd(!)` +is currently only implemented for `InnerProductStyle(t) === EuclideanInnerProduct()`. +""" tsvd, tsvd! + +@doc """ + eig(t::AbstractTensorMap, [(leftind, rightind)::Index2Tuple]; kwargs...) -> D, V + eig!(t::AbstractTensorMap; kwargs...) -> D, V + +Compute eigenvalue factorization of tensor `t` as linear map from `rightind` to `leftind`. +The function `eig` assumes that the linear map is not hermitian and returns type stable +complex valued `D` and `V` tensors for both real and complex valued `t`. See `eigh` for +hermitian linear maps + +If `leftind` and `rightind` are not specified, the current partition of left and right +indices of `t` is used. In that case, less memory is allocated if one allows the data in +`t` to be destroyed/overwritten, by using `eig!(t)`. Note that the permuted tensor on +which `eig!` is called should have equal domain and codomain, as otherwise the eigenvalue +decomposition is meaningless and cannot satisfy +``` +permute(t, (leftind, rightind)) * V = V * D +``` + +Accepts the same keyword arguments `scale` and `permute` as `eigen` of dense +matrices. See the corresponding documentation for more information. + +See also `eigen` and `eigh`. +""" eig + +@doc """ + eigh(t::AbstractTensorMap, [(leftind, rightind)::Index2Tuple]; kwargs...) -> D, V + eigh!(t::AbstractTensorMap; kwargs...) -> D, V + +Compute eigenvalue factorization of tensor `t` as linear map from `rightind` to `leftind`. +The function `eigh` assumes that the linear map is hermitian and `D` and `V` tensors with +the same `scalartype` as `t`. See `eig` and `eigen` for non-hermitian tensors. Hermiticity +requires that the tensor acts on inner product spaces, and the current implementation +requires `InnerProductStyle(t) === EuclideanInnerProduct()`. + +If `leftind` and `rightind` are not specified, the current partition of left and right +indices of `t` is used. In that case, less memory is allocated if one allows the data in +`t` to be destroyed/overwritten, by using `eigh!(t)`. Note that the permuted tensor on +which `eigh!` is called should have equal domain and codomain, as otherwise the eigenvalue +decomposition is meaningless and cannot satisfy +``` +permute(t, (leftind, rightind)) * V = V * D +``` + +See also `eigen` and `eig`. +""" eigh, eigh! + +@doc """ + leftorth(t::AbstractTensorMap, (leftind, rightind)::Index2Tuple; + alg::OrthogonalFactorizationAlgorithm = QRpos()) -> Q, R + +Create orthonormal basis `Q` for indices in `leftind`, and remainder `R` such that +`permute(t, (leftind, rightind)) = Q*R`. + +If `leftind` and `rightind` are not specified, the current partition of left and right +indices of `t` is used. In that case, less memory is allocated if one allows the data in `t` +to be destroyed/overwritten, by using `leftorth!(t, alg = QRpos())`. + +Different algorithms are available, namely `QR()`, `QRpos()`, `SVD()` and `Polar()`. `QR()` +and `QRpos()` use a standard QR decomposition, producing an upper triangular matrix `R`. +`Polar()` produces a Hermitian and positive semidefinite `R`. `QRpos()` corrects the +standard QR decomposition such that the diagonal elements of `R` are positive. Only +`QRpos()` and `Polar()` are unique (no residual freedom) so that they always return the same +result for the same input tensor `t`. + +Orthogonality requires `InnerProductStyle(t) <: HasInnerProduct`, and +`leftorth(!)` is currently only implemented for + `InnerProductStyle(t) === EuclideanInnerProduct()`. +""" leftorth, leftorth! + +@doc """ + rightorth(t::AbstractTensorMap, [(leftind, rightind)::Index2Tuple]; + alg::OrthogonalFactorizationAlgorithm = LQpos()) -> L, Q + rightorth!(t::AbstractTensorMap; alg) -> L, Q + +Create orthonormal basis `Q` for indices in `rightind`, and remainder `L` such that +`permute(t, (leftind, rightind)) = L*Q`. + +If `leftind` and `rightind` are not specified, the current partition of left and right +indices of `t` is used. In that case, less memory is allocated if one allows the data in `t` +to be destroyed/overwritten, by using `rightorth!(t, alg = LQpos())`. + +Different algorithms are available, namely `LQ()`, `LQpos()`, `RQ()`, `RQpos()`, `SVD()` and +`Polar()`. `LQ()` and `LQpos()` produce a lower triangular matrix `L` and are computed using +a QR decomposition of the transpose. `RQ()` and `RQpos()` produce an upper triangular +remainder `L` and only works if the total left dimension is smaller than or equal to the +total right dimension. `LQpos()` and `RQpos()` add an additional correction such that the +diagonal elements of `L` are positive. `Polar()` produces a Hermitian and positive +semidefinite `L`. Only `LQpos()`, `RQpos()` and `Polar()` are unique (no residual freedom) +so that they always return the same result for the same input tensor `t`. + +Orthogonality requires `InnerProductStyle(t) <: HasInnerProduct`, and +`rightorth(!)` is currently only implemented for +`InnerProductStyle(t) === EuclideanInnerProduct()`. +""" rightorth, rightorth! + +@doc """ + leftnull(t::AbstractTensorMap, [(leftind, rightind)::Index2Tuple]; + alg::OrthogonalFactorizationAlgorithm = QRpos()) -> N + leftnull!(t::AbstractTensorMap; alg) -> N + +Create orthonormal basis for the orthogonal complement of the support of the indices in +`leftind`, such that `N' * permute(t, (leftind, rightind)) = 0`. + +If `leftind` and `rightind` are not specified, the current partition of left and right +indices of `t` is used. In that case, less memory is allocated if one allows the data in `t` +to be destroyed/overwritten, by using `leftnull!(t, alg = QRpos())`. + +Different algorithms are available, namely `QR()` (or equivalently, `QRpos()`), `SVD()` and +`SDD()`. The first assumes that the matrix is full rank and requires `iszero(atol)` and +`iszero(rtol)`. With `SVD()` and `SDD()`, `rightnull` will use the corresponding singular +value decomposition, and one can specify an absolute or relative tolerance for which +singular values are to be considered zero, where `max(atol, norm(t)*rtol)` is used as upper +bound. + +Orthogonality requires `InnerProductStyle(t) <: HasInnerProduct`, and +`leftnull(!)` is currently only implemented for +`InnerProductStyle(t) === EuclideanInnerProduct()`. +""" leftnull, leftnull! + +@doc """ + rightnull(t::AbstractTensorMap, [(leftind, rightind)::Index2Tuple]; + alg::OrthogonalFactorizationAlgorithm = LQ(), + atol::Real = 0.0, + rtol::Real = eps(real(float(one(scalartype(t)))))*iszero(atol)) -> N + rightnull!(t::AbstractTensorMap; alg, atol, rtol) + +Create orthonormal basis for the orthogonal complement of the support of the indices in +`rightind`, such that `permute(t, (leftind, rightind))*N' = 0`. + +If `leftind` and `rightind` are not specified, the current partition of left and right +indices of `t` is used. In that case, less memory is allocated if one allows the data in `t` +to be destroyed/overwritten, by using `rightnull!(t, alg = LQpos())`. + +Different algorithms are available, namely `LQ()` (or equivalently, `LQpos`), `SVD()` and +`SDD()`. The first assumes that the matrix is full rank and requires `iszero(atol)` and +`iszero(rtol)`. With `SVD()` and `SDD()`, `rightnull` will use the corresponding singular +value decomposition, and one can specify an absolute or relative tolerance for which +singular values are to be considered zero, where `max(atol, norm(t)*rtol)` is used as upper +bound. + +Orthogonality requires `InnerProductStyle(t) <: HasInnerProduct`, and +`rightnull(!)` is currently only implemented for +`InnerProductStyle(t) === EuclideanInnerProduct()`. +""" rightnull, rightnull! + +@doc """ + eigen(t::AbstractTensorMap, [(leftind, rightind)::Index2Tuple]; kwargs...) -> D, V + eigen!(t::AbstractTensorMap; kwargs...) -> D, V + +Compute eigenvalue factorization of tensor `t` as linear map from `rightind` to `leftind`. + +If `leftind` and `rightind` are not specified, the current partition of left and right +indices of `t` is used. In that case, less memory is allocated if one allows the data in `t` +to be destroyed/overwritten, by using `eigen!(t)`. Note that the permuted tensor on which +`eigen!` is called should have equal domain and codomain, as otherwise the eigenvalue +decomposition is meaningless and cannot satisfy +``` +permute(t, (leftind, rightind)) * V = V * D +``` + +Accepts the same keyword arguments `scale` and `permute` as `eigen` of dense +matrices. See the corresponding documentation for more information. + +See also [`eig(!)`](@ref eig) and [`eigh(!)`](@ref) +""" eigen(::AbstractTensorMap), eigen!(::AbstractTensorMap) + +@doc """ + isposdef(t::AbstractTensor, [(leftind, rightind)::Index2Tuple]) -> ::Bool + +Test whether a tensor `t` is positive definite as linear map from `rightind` to `leftind`. + +If `leftind` and `rightind` are not specified, the current partition of left and right +indices of `t` is used. In that case, less memory is allocated if one allows the data in +`t` to be destroyed/overwritten, by using `isposdef!(t)`. Note that the permuted tensor on +which `isposdef!` is called should have equal domain and codomain, as otherwise it is +meaningless. +""" isposdef(::AbstractTensorMap), isposdef!(::AbstractTensorMap) + +for f in + (:tsvd, :eig, :eigh, :eigen, :leftorth, :rightorth, :left_polar, :right_polar, + :leftnull, + :rightnull, :isposdef) + f! = Symbol(f, :!) + @eval function $f(t::AbstractTensorMap, p::Index2Tuple; kwargs...) + tcopy = permutedcopy_oftype(t, factorisation_scalartype($f, t), p) + return $f!(tcopy; kwargs...) + end + @eval function $f(t::AbstractTensorMap; kwargs...) + tcopy = copy_oftype(t, factorisation_scalartype($f, t)) + return $f!(tcopy; kwargs...) + end +end + +function LinearAlgebra.eigvals(t::AbstractTensorMap; kwargs...) + tcopy = copy_oftype(t, factorisation_scalartype(eigen, t)) + return LinearAlgebra.eigvals!(tcopy; kwargs...) +end + +function LinearAlgebra.svdvals(t::AbstractTensorMap) + tcopy = copy_oftype(t, factorisation_scalartype(tsvd, t)) + return LinearAlgebra.svdvals!(tcopy) +end diff --git a/src/tensors/factorizations/matrixalgebrakit.jl b/src/tensors/factorizations/matrixalgebrakit.jl new file mode 100644 index 000000000..b4fb41e9b --- /dev/null +++ b/src/tensors/factorizations/matrixalgebrakit.jl @@ -0,0 +1,562 @@ +# Algorithm selection +# ------------------- +for f! in + [:svd_compact!, :svd_full!, :svd_trunc!, :svd_vals!, :qr_compact!, :qr_full!, :qr_null!, + :lq_compact!, :lq_full!, :lq_null!, :eig_full!, :eig_trunc!, :eig_vals!, :eigh_full!, + :eigh_trunc!, :eigh_vals!, :left_polar!, :right_polar!] + @eval function default_algorithm(::typeof($f!), ::Type{T}; + kwargs...) where {T<:AbstractTensorMap} + return default_algorithm($f!, blocktype(T); kwargs...) + end +end + +function _select_truncation(f, ::AbstractTensorMap, + trunc::MatrixAlgebraKit.TruncationStrategy) + return trunc +end +function _select_truncation(::typeof(left_null!), ::AbstractTensorMap, trunc::NamedTuple) + return MatrixAlgebraKit.null_truncation_strategy(; trunc...) +end + +# Generic Implementations +# ----------------------_ +for f! in (:qr_compact!, :qr_full!, + :lq_compact!, :lq_full!, + :eig_full!, :eigh_full!, + :svd_compact!, :svd_full!, + :left_polar!, :left_orth_polar!, + :right_polar!, :right_orth_polar!, + :left_orth!, :right_orth!) + @eval function $f!(t::AbstractTensorMap, F, alg::AbstractAlgorithm) + check_input($f!, t, F, alg) + + foreachblock(t, F...) do _, bs + factors = Base.tail(bs) + factors′ = $f!(first(bs), factors, alg) + # deal with the case where the output is not in-place + for (f′, f) in zip(factors′, factors) + f′ === f || copyto!(f, f′) + end + return nothing + end + + return F + end +end + +# Handle these separately because single output instead of tuple +for f! in (:qr_null!, :lq_null!, :svd_vals!, :eig_vals!, :eigh_vals!) + @eval function $f!(t::AbstractTensorMap, N, alg::AbstractAlgorithm) + check_input($f!, t, N, alg) + + foreachblock(t, N) do _, (b, n) + n′ = $f!(b, n, alg) + # deal with the case where the output is not the same as the input + n === n′ || copyto!(n, n′) + return nothing + end + + return N + end +end + +# Singular value decomposition +# ---------------------------- +const _T_USVᴴ = Tuple{<:AbstractTensorMap,<:AbstractTensorMap,<:AbstractTensorMap} +const _T_USVᴴ_diag = Tuple{<:AbstractTensorMap,<:DiagonalTensorMap,<:AbstractTensorMap} + +function check_input(::typeof(svd_full!), t::AbstractTensorMap, (U, S, Vᴴ)::_T_USVᴴ, + ::AbstractAlgorithm) + # scalartype checks + @check_scalar U t + @check_scalar S t real + @check_scalar Vᴴ t + + # space checks + V_cod = fuse(codomain(t)) + V_dom = fuse(domain(t)) + @check_space(U, codomain(t) ← V_cod) + @check_space(S, V_cod ← V_dom) + @check_space(Vᴴ, V_dom ← domain(t)) + + return nothing +end + +function check_input(::typeof(svd_compact!), t::AbstractTensorMap, (U, S, Vᴴ)::_T_USVᴴ_diag, + ::AbstractAlgorithm) + # scalartype checks + @check_scalar U t + @check_scalar S t real + @check_scalar Vᴴ t + + # space checks + V_cod = V_dom = infimum(fuse(codomain(t)), fuse(domain(t))) + @check_space(U, codomain(t) ← V_cod) + @check_space(S, V_cod ← V_dom) + @check_space(Vᴴ, V_dom ← domain(t)) + + return nothing +end + +function check_input(::typeof(svd_vals!), t::AbstractTensorMap, S::SectorDict, + ::AbstractAlgorithm) + @check_scalar S t real + V_cod = infimum(fuse(codomain(t)), fuse(domain(t))) + @check_space(S, V_cod ← V_dom) + return nothing +end + +function initialize_output(::typeof(svd_full!), t::AbstractTensorMap, ::AbstractAlgorithm) + V_cod = fuse(codomain(t)) + V_dom = fuse(domain(t)) + U = similar(t, codomain(t) ← V_cod) + S = similar(t, real(scalartype(t)), V_cod ← V_dom) + Vᴴ = similar(t, V_dom ← domain(t)) + return U, S, Vᴴ +end + +function initialize_output(::typeof(svd_compact!), t::AbstractTensorMap, + ::AbstractAlgorithm) + V_cod = V_dom = infimum(fuse(codomain(t)), fuse(domain(t))) + U = similar(t, codomain(t) ← V_cod) + S = DiagonalTensorMap{real(scalartype(t))}(undef, V_cod) + Vᴴ = similar(t, V_dom ← domain(t)) + return U, S, Vᴴ +end + +function initialize_output(::typeof(svd_trunc!), t::AbstractTensorMap, + alg::TruncatedAlgorithm) + return initialize_output(svd_compact!, t, alg.alg) +end + +function initialize_output(::typeof(svd_vals!), t::AbstractTensorMap, + alg::AbstractAlgorithm) + V_cod = infimum(fuse(codomain(t)), fuse(domain(t))) + return DiagonalTensorMap{real(scalartype(t))}(undef, V_cod) +end + +function svd_trunc!(t::AbstractTensorMap, USVᴴ, alg::TruncatedAlgorithm) + USVᴴ′ = svd_compact!(t, USVᴴ, alg.alg) + return truncate!(svd_trunc!, USVᴴ′, alg.trunc) +end + +# Eigenvalue decomposition +# ------------------------ +const _T_DV = Tuple{<:DiagonalTensorMap,<:AbstractTensorMap} + +function check_input(::typeof(eigh_full!), t::AbstractTensorMap, (D, V)::_T_DV, + ::AbstractAlgorithm) + domain(t) == codomain(t) || + throw(ArgumentError("Eigenvalue decomposition requires square input tensor")) + + # scalartype checks + @check_scalar D t real + @check_scalar V t + + # space checks + V_D = fuse(domain(t)) + @check_space(D, V_D ← V_D) + @check_space(V, codomain(t) ← V_D) + + return nothing +end + +function check_input(::typeof(eig_full!), t::AbstractTensorMap, (D, V)::_T_DV, + ::AbstractAlgorithm) + domain(t) == codomain(t) || + throw(ArgumentError("Eigenvalue decomposition requires square input tensor")) + + # scalartype checks + @check_scalar D t complex + @check_scalar V t complex + + # space checks + V_D = fuse(domain(t)) + @check_space(D, V_D ← V_D) + @check_space(V, codomain(t) ← V_D) + + return nothing +end + +function check_input(::typeof(eigh_vals!), t::AbstractTensorMap, D::DiagonalTensorMap, + ::AbstractAlgorithm) + @check_scalar D t real + V_D = fuse(domain(t)) + @check_space(D, V_D ← V_D) + return nothing +end + +function check_input(::typeof(eig_vals!), t::AbstractTensorMap, D::DiagonalTensorMap, + ::AbstractAlgorithm) + @check_scalar D t complex + V_D = fuse(domain(t)) + @check_space(D, V_D ← V_D) + return nothing +end + +function initialize_output(::typeof(eigh_full!), t::AbstractTensorMap, ::AbstractAlgorithm) + V_D = fuse(domain(t)) + T = real(scalartype(t)) + D = DiagonalTensorMap{T}(undef, V_D) + V = similar(t, codomain(t) ← V_D) + return D, V +end + +function initialize_output(::typeof(eig_full!), t::AbstractTensorMap, ::AbstractAlgorithm) + V_D = fuse(domain(t)) + Tc = complex(scalartype(t)) + D = DiagonalTensorMap{Tc}(undef, V_D) + V = similar(t, Tc, codomain(t) ← V_D) + return D, V +end + +function initialize_output(::typeof(eigh_vals!), t::AbstractTensorMap, + alg::AbstractAlgorithm) + V_D = fuse(domain(t)) + T = real(scalartype(t)) + return D = DiagonalTensorMap{Tc}(undef, V_D) +end + +function initialize_output(::typeof(eig_vals!), t::AbstractTensorMap, + alg::AbstractAlgorithm) + V_D = fuse(domain(t)) + Tc = complex(scalartype(t)) + return D = DiagonalTensorMap{Tc}(undef, V_D) +end + +function initialize_output(::typeof(eigh_trunc!), t::AbstractTensorMap, + alg::TruncatedAlgorithm) + return initialize_output(eigh_full!, t, alg.alg) +end + +function initialize_output(::typeof(eig_trunc!), t::AbstractTensorMap, + alg::TruncatedAlgorithm) + return initialize_output(eig_full!, t, alg.alg) +end + +function eigh_trunc!(t::AbstractTensorMap, DV, alg::TruncatedAlgorithm) + DV′ = eigh_full!(t, DV, alg.alg) + return truncate!(eigh_trunc!, DV′, alg.trunc) +end + +function eig_trunc!(t::AbstractTensorMap, DV, alg::TruncatedAlgorithm) + DV′ = eig_full!(t, DV, alg.alg) + return truncate!(eig_trunc!, DV′, alg.trunc) +end + +# QR decomposition +# ---------------- +const _T_QR = Tuple{<:AbstractTensorMap,<:AbstractTensorMap} + +function check_input(::typeof(qr_full!), t::AbstractTensorMap, (Q, R)::_T_QR, + ::AbstractAlgorithm) + # scalartype checks + @check_scalar Q t + @check_scalar R t + + # space checks + V_Q = fuse(codomain(t)) + @check_space(Q, codomain(t) ← V_Q) + @check_space(R, V_Q ← domain(t)) + + return nothing +end + +function check_input(::typeof(qr_compact!), t::AbstractTensorMap, (Q, R)::_T_QR, + ::AbstractAlgorithm) + # scalartype checks + @check_scalar Q t + @check_scalar R t + + # space checks + V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) + @check_space(Q, codomain(t) ← V_Q) + @check_space(R, V_Q ← domain(t)) + + return nothing +end + +function check_input(::typeof(qr_null!), t::AbstractTensorMap, N::AbstractTensorMap, + ::AbstractAlgorithm) + # scalartype checks + @check_scalar N t + + # space checks + V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) + V_N = ⊖(fuse(codomain(t)), V_Q) + @check_space(N, codomain(t) ← V_N) + + return nothing +end + +function initialize_output(::typeof(qr_full!), t::AbstractTensorMap, ::AbstractAlgorithm) + V_Q = fuse(codomain(t)) + Q = similar(t, codomain(t) ← V_Q) + R = similar(t, V_Q ← domain(t)) + return Q, R +end + +function initialize_output(::typeof(qr_compact!), t::AbstractTensorMap, ::AbstractAlgorithm) + V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) + Q = similar(t, codomain(t) ← V_Q) + R = similar(t, V_Q ← domain(t)) + return Q, R +end + +function initialize_output(::typeof(qr_null!), t::AbstractTensorMap, ::AbstractAlgorithm) + V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) + V_N = ⊖(fuse(codomain(t)), V_Q) + N = similar(t, codomain(t) ← V_N) + return N +end + +# LQ decomposition +# ---------------- +const _T_LQ = Tuple{<:AbstractTensorMap,<:AbstractTensorMap} + +function check_input(::typeof(lq_full!), t::AbstractTensorMap, (L, Q)::_T_LQ, + ::AbstractAlgorithm) + # scalartype checks + @check_scalar L t + @check_scalar Q t + + # space checks + V_Q = fuse(domain(t)) + @check_space(L, codomain(t) ← V_Q) + @check_space(Q, V_Q ← domain(t)) + + return nothing +end + +function check_input(::typeof(lq_compact!), t::AbstractTensorMap, (L, Q)::_T_LQ, + ::AbstractAlgorithm) + # scalartype checks + @check_scalar L t + @check_scalar Q t + + # space checks + V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) + @check_space(L, codomain(t) ← V_Q) + @check_space(Q, V_Q ← domain(t)) + + return nothing +end + +function check_input(::typeof(lq_null!), t::AbstractTensorMap, N, ::AbstractAlgorithm) + # scalartype checks + @check_scalar N t + + # space checks + V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) + V_N = ⊖(fuse(domain(t)), V_Q) + @check_space(N, V_N ← domain(t)) + + return nothing +end + +function initialize_output(::typeof(lq_full!), t::AbstractTensorMap, ::AbstractAlgorithm) + V_Q = fuse(domain(t)) + L = similar(t, codomain(t) ← V_Q) + Q = similar(t, V_Q ← domain(t)) + return L, Q +end + +function initialize_output(::typeof(lq_compact!), t::AbstractTensorMap, ::AbstractAlgorithm) + V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) + L = similar(t, codomain(t) ← V_Q) + Q = similar(t, V_Q ← domain(t)) + return L, Q +end + +function initialize_output(::typeof(lq_null!), t::AbstractTensorMap, ::AbstractAlgorithm) + V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) + V_N = ⊖(fuse(domain(t)), V_Q) + N = similar(t, V_N ← domain(t)) + return N +end + +# Polar decomposition +# ------------------- +const _T_WP = Tuple{<:AbstractTensorMap,<:AbstractTensorMap} +const _T_PWᴴ = Tuple{<:AbstractTensorMap,<:AbstractTensorMap} + +function check_input(::typeof(left_polar!), t::AbstractTensorMap, (W, P)::_T_WP, + ::AbstractAlgorithm) + codomain(t) ≿ domain(t) || + throw(ArgumentError("Polar decomposition requires `codomain(t) ≿ domain(t)`")) + + # scalartype checks + @check_scalar W t + @check_scalar P t + + # space checks + @check_space(W, space(t)) + @check_space(P, domain(t) ← domain(t)) + + return nothing +end + +function check_input(::typeof(left_orth_polar!), t::AbstractTensorMap, (W, P)::_T_WP, + ::AbstractAlgorithm) + codomain(t) ≿ domain(t) || + throw(ArgumentError("Polar decomposition requires `codomain(t) ≿ domain(t)`")) + + # scalartype checks + @check_scalar W t + @check_scalar P t + + # space checks + VW = fuse(domain(t)) + @check_space(W, codomain(t) ← VW) + @check_space(P, VW ← domain(t)) + + return nothing +end + +function initialize_output(::typeof(left_polar!), t::AbstractTensorMap, ::AbstractAlgorithm) + W = similar(t, space(t)) + P = similar(t, domain(t) ← domain(t)) + return W, P +end + +function check_input(::typeof(right_polar!), t::AbstractTensorMap, (P, Wᴴ)::_T_PWᴴ, + ::AbstractAlgorithm) + codomain(t) ≾ domain(t) || + throw(ArgumentError("Polar decomposition requires `domain(t) ≿ codomain(t)`")) + + # scalartype checks + @check_scalar P t + @check_scalar Wᴴ t + + # space checks + @check_space(P, codomain(t) ← codomain(t)) + @check_space(Wᴴ, space(t)) + + return nothing +end + +function check_input(::typeof(right_orth_polar!), t::AbstractTensorMap, (P, Wᴴ)::_T_PWᴴ, + ::AbstractAlgorithm) + codomain(t) ≾ domain(t) || + throw(ArgumentError("Polar decomposition requires `domain(t) ≿ codomain(t)`")) + + # scalartype checks + @check_scalar P t + @check_scalar Wᴴ t + + # space checks + VW = fuse(codomain(t)) + @check_space(P, codomain(t) ← VW) + @check_space(Wᴴ, VW ← domain(t)) + + return nothing +end + +function initialize_output(::typeof(right_polar!), t::AbstractTensorMap, + ::AbstractAlgorithm) + P = similar(t, codomain(t) ← codomain(t)) + Wᴴ = similar(t, space(t)) + return Wᴴ, P +end + +# Needed to get algorithm selection to behave +function left_orth_polar!(t::AbstractTensorMap, VC, alg) + alg′ = MatrixAlgebraKit.select_algorithm(left_polar!, t, alg) + return left_orth_polar!(t, VC, alg′) +end +function right_orth_polar!(t::AbstractTensorMap, CVᴴ, alg) + alg′ = MatrixAlgebraKit.select_algorithm(right_polar!, t, alg) + return right_orth_polar!(t, CVᴴ, alg′) +end + +# Orthogonalization +# ----------------- +const _T_VC = Tuple{<:AbstractTensorMap,<:AbstractTensorMap} +const _T_CVᴴ = Tuple{<:AbstractTensorMap,<:AbstractTensorMap} + +function check_input(::typeof(left_orth!), t::AbstractTensorMap, (V, C)::_T_VC, + ::AbstractAlgorithm) + # scalartype checks + @check_scalar V t + isnothing(C) || @check_scalar C t + + # space checks + V_C = infimum(fuse(codomain(t)), fuse(domain(t))) + @check_space(V, codomain(t) ← V_C) + isnothing(C) || @check_space(C, V_C ← domain(t)) + + return nothing +end + +function check_input(::typeof(right_orth!), t::AbstractTensorMap, (C, Vᴴ)::_T_CVᴴ, + ::AbstractAlgorithm) + # scalartype checks + isnothing(C) || @check_scalar C t + @check_scalar Vᴴ t + + # space checks + V_C = infimum(fuse(codomain(t)), fuse(domain(t))) + isnothing(C) || @check_space(C, codomain(t) ← V_C) + @check_space(Vᴴ, V_C ← domain(t)) + + return nothing +end + +function initialize_output(::typeof(left_orth!), t::AbstractTensorMap) + V_C = infimum(fuse(codomain(t)), fuse(domain(t))) + V = similar(t, codomain(t) ← V_C) + C = similar(t, V_C ← domain(t)) + return V, C +end + +function initialize_output(::typeof(right_orth!), t::AbstractTensorMap) + V_C = infimum(fuse(codomain(t)), fuse(domain(t))) + C = similar(t, codomain(t) ← V_C) + Vᴴ = similar(t, V_C ← domain(t)) + return C, Vᴴ +end + +# Nullspace +# --------- +function check_input(::typeof(left_null!), t::AbstractTensorMap, N, ::AbstractAlgorithm) + # scalartype checks + @check_scalar N t + + # space checks + V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) + V_N = ⊖(fuse(codomain(t)), V_Q) + @check_space(N, codomain(t) ← V_N) + + return nothing +end + +function check_input(::typeof(right_null!), t::AbstractTensorMap, N, ::AbstractAlgorithm) + @check_scalar N t + + # space checks + V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) + V_N = ⊖(fuse(domain(t)), V_Q) + @check_space(N, V_N ← domain(t)) + + return nothing +end + +function initialize_output(::typeof(left_null!), t::AbstractTensorMap) + V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) + V_N = ⊖(fuse(codomain(t)), V_Q) + N = similar(t, codomain(t) ← V_N) + return N +end + +function initialize_output(::typeof(right_null!), t::AbstractTensorMap) + V_Q = infimum(fuse(codomain(t)), fuse(domain(t))) + V_N = ⊖(fuse(domain(t)), V_Q) + N = similar(t, V_N ← domain(t)) + return N +end + +for (f!, f_svd!) in zip((:left_null!, :right_null!), (:left_null_svd!, :right_null_svd!)) + @eval function $f_svd!(t::AbstractTensorMap, N, alg, ::Nothing=nothing) + return $f!(t, N; alg_svd=alg) + end +end diff --git a/src/tensors/factorizations/truncation.jl b/src/tensors/factorizations/truncation.jl new file mode 100644 index 000000000..d553c5c93 --- /dev/null +++ b/src/tensors/factorizations/truncation.jl @@ -0,0 +1,235 @@ +# Strategies +# ---------- +notrunc() = NoTruncation() + +# deprecate +const TruncationScheme = TruncationStrategy +@deprecate truncdim(d::Int) truncrank(d) +@deprecate truncbelow(ϵ::Real, add_back::Int=0) trunctol(ϵ) + +# TODO: add this to MatrixAlgebraKit +struct TruncationError{T<:Real} <: TruncationStrategy + ϵ::T + p::Real +end +truncerr(epsilon::Real, p::Real=2) = TruncationError(epsilon, p) + +struct TruncationSpace{S<:ElementarySpace} <: TruncationStrategy + space::S +end +truncspace(space::ElementarySpace) = TruncationSpace(space) + +# Truncation +# ---------- +function truncate!(::typeof(svd_trunc!), (U, S, Vᴴ)::_T_USVᴴ, strategy::TruncationStrategy) + ind = findtruncated_sorted(diagview(S), strategy) + V_truncated = spacetype(S)(c => length(I) for (c, I) in ind) + + Ũ = similar(U, codomain(U) ← V_truncated) + for (c, b) in blocks(Ũ) + I = get(ind, c, nothing) + @assert !isnothing(I) + copy!(b, @view(block(U, c)[:, I])) + end + + S̃ = DiagonalTensorMap{scalartype(S)}(undef, V_truncated) + for (c, b) in blocks(S̃) + I = get(ind, c, nothing) + @assert !isnothing(I) + copy!(b.diag, @view(block(S, c).diag[I])) + end + + Ṽᴴ = similar(Vᴴ, V_truncated ← domain(Vᴴ)) + for (c, b) in blocks(Ṽᴴ) + I = get(ind, c, nothing) + @assert !isnothing(I) + copy!(b, @view(block(Vᴴ, c)[I, :])) + end + + return Ũ, S̃, Ṽᴴ +end + +function truncate!(::typeof(left_null!), + (U, S)::Tuple{<:AbstractTensorMap, + <:AbstractTensorMap}, + strategy::MatrixAlgebraKit.TruncationStrategy) + extended_S = SectorDict(c => vcat(diagview(b), + zeros(eltype(b), max(0, size(b, 2) - size(b, 1)))) + for (c, b) in blocks(S)) + ind = findtruncated(extended_S, strategy) + V_truncated = spacetype(S)(c => length(axes(b, 1)[ind[c]]) for (c, b) in blocks(S)) + Ũ = similar(U, codomain(U) ← V_truncated) + for (c, b) in blocks(Ũ) + copy!(b, @view(block(U, c)[:, ind[c]])) + end + return Ũ +end + +function truncate!(::typeof(eigh_trunc!), (D, V)::_T_DV, strategy::TruncationStrategy) + ind = findtruncated(diagview(D), strategy) + V_truncated = spacetype(D)(c => length(I) for (c, I) in ind) + + D̃ = DiagonalTensorMap{scalartype(D)}(undef, V_truncated) + for (c, b) in blocks(D̃) + I = get(ind, c, nothing) + @assert !isnothing(I) + copy!(b.diag, @view(block(D, c).diag[I])) + end + + Ṽ = similar(V, V_truncated ← domain(V)) + for (c, b) in blocks(Ṽ) + I = get(ind, c, nothing) + @assert !isnothing(I) + copy!(b, @view(block(V, c)[I, :])) + end + + return D̃, Ṽ +end +function truncate!(::typeof(eig_trunc!), (D, V)::_T_DV, strategy::TruncationStrategy) + ind = findtruncated(diagview(D), strategy) + V_truncated = spacetype(D)(c => length(I) for (c, I) in ind) + + D̃ = DiagonalTensorMap{scalartype(D)}(undef, V_truncated) + for (c, b) in blocks(D̃) + I = get(ind, c, nothing) + @assert !isnothing(I) + copy!(b.diag, @view(block(D, c).diag[I])) + end + + Ṽ = similar(V, V_truncated ← domain(V)) + for (c, b) in blocks(Ṽ) + I = get(ind, c, nothing) + @assert !isnothing(I) + copy!(b, @view(block(V, c)[I, :])) + end + + return D̃, Ṽ +end + +# Find truncation +# --------------- +# auxiliary functions +rtol_to_atol(S, p, atol, rtol) = rtol > 0 ? max(atol, _norm(S, p) * rtol) : atol + +function _compute_truncerr(Σdata, truncdim, p=2) + I = keytype(Σdata) + S = scalartype(valtype(Σdata)) + return TensorKit._norm((c => @view(v[(get(truncdim, c, 0) + 1):end]) + for (c, v) in Σdata), + p, zero(S)) +end + +function _findnexttruncvalue(S, truncdim::SectorDict{I,Int}; by=identity, + rev::Bool=true) where {I<:Sector} + # early return + (isempty(S) || all(iszero, values(truncdim))) && return nothing + if rev + σmin, imin = findmin(keys(truncdim)) do c + d = truncdim[c] + return by(S[c][d]) + end + return σmin, keys(truncdim)[imin] + else + σmax, imax = findmax(keys(truncdim)) do c + d = truncdim[c] + return by(S[c][d]) + end + return σmax, keys(truncdim)[imax] + end +end + +# implementations +function findtruncated_sorted(S::SectorDict, strategy::TruncationKeepAbove) + atol = rtol_to_atol(S, strategy.p, strategy.atol, strategy.rtol) + findtrunc = Base.Fix2(findtruncated_sorted, truncbelow(atol)) + return SectorDict(c => findtrunc(d) for (c, d) in Sd) +end +function findtruncated(S::SectorDict, strategy::TruncationKeepAbove) + atol = rtol_to_atol(S, strategy.p, strategy.atol, strategy.rtol) + findtrunc = Base.Fix2(findtruncated, truncbelow(atol)) + return SectorDict(c => findtrunc(d) for (c, d) in Sd) +end + +function findtruncated_sorted(S::SectorDict, strategy::TruncationKeepBelow) + atol = rtol_to_atol(S, strategy.p, strategy.atol, strategy.rtol) + findtrunc = Base.Fix2(findtruncated_sorted, truncabove(atol)) + return SectorDict(c => findtrunc(d) for (c, d) in Sd) +end +function findtruncated(S::SectorDict, strategy::TruncationKeepBelow) + atol = rtol_to_atol(S, strategy.p, strategy.atol, strategy.rtol) + findtrunc = Base.Fix2(findtruncated, truncabove(atol)) + return SectorDict(c => findtrunc(d) for (c, d) in Sd) +end + +function findtruncated_sorted(Sd::SectorDict, strategy::TruncationError) + I = keytype(Sd) + truncdim = SectorDict{I,Int}(c => length(d) for (c, d) in Sd) + while true + next = _findnexttruncvalue(Sd, truncdim) + isnothing(next) && break + σmin, cmin = next + truncdim[cmin] -= 1 + err = _compute_truncerr(Sd, truncdim, strategy.p) + if err > strategy.ϵ + truncdim[cmin] += 1 + break + end + if truncdim[cmin] == 0 + delete!(truncdim, cmin) + end + end + return SectorDict{I,Base.OneTo{Int}}(c => Base.OneTo(d) for (c, d) in truncdim) +end + +function findtruncated_sorted(Sd::SectorDict, strategy::TruncationKeepSorted) + return findtruncated(Sd, strategy) +end +function findtruncated(Sd::SectorDict, strategy::TruncationKeepSorted) + permutations = SectorDict(c => (sortperm(d; strategy.by, strategy.rev)) + for (c, d) in Sd) + Sd = SectorDict(c => sort(d; strategy.by, strategy.rev) for (c, d) in Sd) + + I = keytype(Sd) + truncdim = SectorDict{I,Int}(c => length(d) for (c, d) in Sd) + totaldim = sum(dim(c) * d for (c, d) in truncdim; init=0) + while true + next = _findnexttruncvalue(Sd, truncdim; strategy.by, strategy.rev) + isnothing(next) && break + _, cmin = next + truncdim[cmin] -= 1 + totaldim -= dim(cmin) + if totaldim < strategy.howmany + truncdim[cmin] += 1 + break + end + if truncdim[cmin] == 0 + delete!(truncdim, cmin) + end + end + return SectorDict(c => permutations[c][Base.OneTo(d)] for (c, d) in truncdim) +end + +function findtruncated_sorted(Sd::SectorDict, strategy::TruncationSpace) + I = keytype(Sd) + return SectorDict{I,Base.OneTo{Int}}(c => Base.OneTo(min(length(d), + dim(strategy.space, c))) + for (c, d) in Sd) +end + +function findtruncated_sorted(Sd::SectorDict, strategy::TruncationKeepFiltered) + return SectorDict(c => findtruncated_sorted(d, strategy) for (c, d) in Sd) +end +function findtruncated(Sd::SectorDict, strategy::TruncationKeepFiltered) + return SectorDict(c => findtruncated(d, strategy) for (c, d) in Sd) +end + +function findtruncated_sorted(Sd::SectorDict, strategy::TruncationIntersection) + inds = map(Base.Fix1(findtruncated_sorted, Sd), strategy) + return SectorDict(c => intersect(map(Base.Fix2(getindex, c), inds)...) + for c in intersect(map(keys, inds)...)) +end +function findtruncated(Sd::SectorDict, strategy::TruncationIntersection) + inds = map(Base.Fix1(findtruncated, Sd), strategy) + return SectorDict(c => intersect(map(Base.Fix2(getindex, c), inds)...) + for c in intersect(map(keys, inds)...)) +end diff --git a/src/tensors/factorizations/utility.jl b/src/tensors/factorizations/utility.jl new file mode 100644 index 000000000..e23fb8b73 --- /dev/null +++ b/src/tensors/factorizations/utility.jl @@ -0,0 +1,29 @@ +# convenience to set default +macro check_space(x, V) + return esc(:($MatrixAlgebraKit.@check_size($x, $V, $space))) +end +macro check_scalar(x, y, op=:identity, eltype=:scalartype) + return esc(:($MatrixAlgebraKit.@check_scalar($x, $y, $op, $eltype))) +end + +function factorisation_scalartype(t::AbstractTensorMap) + T = scalartype(t) + return promote_type(Float32, typeof(zero(T) / sqrt(abs2(one(T))))) +end +factorisation_scalartype(f, t) = factorisation_scalartype(t) + +function permutedcopy_oftype(t::AbstractTensorMap, T::Type{<:Number}, p::Index2Tuple) + return permute!(similar(t, T, permute(space(t), p)), t, p) +end +function copy_oftype(t::AbstractTensorMap, T::Type{<:Number}) + return copy!(similar(t, T), t) +end + +function _reverse!(t::AbstractTensorMap; dims=:) + for (c, b) in blocks(t) + reverse!(b; dims) + end + return t +end + +diagview(t::AbstractTensorMap) = SectorDict(c => diagview(b) for (c, b) in blocks(t)) diff --git a/src/tensors/linalg.jl b/src/tensors/linalg.jl index f29bdf809..51567dec7 100644 --- a/src/tensors/linalg.jl +++ b/src/tensors/linalg.jl @@ -59,7 +59,7 @@ function one!(t::AbstractTensorMap) domain(t) == codomain(t) || throw(SectorMismatch("no identity if domain and codomain are different")) for (c, b) in blocks(t) - MatrixAlgebra.one!(b) + one!(b) end return t end @@ -106,7 +106,7 @@ function isomorphism!(t::AbstractTensorMap) domain(t) ≅ codomain(t) || throw(SpaceMismatch(lazy"domain and codomain are not isomorphic: $(space(t))")) for (_, b) in blocks(t) - MatrixAlgebra.one!(b) + one!(b) end return t end @@ -155,7 +155,7 @@ function isometry!(t::AbstractTensorMap) domain(t) ≾ codomain(t) || throw(SpaceMismatch(lazy"domain and codomain are not isometrically embeddable: $(space(t))")) for (_, b) in blocks(t) - MatrixAlgebra.one!(b) + one!(b) end return t end @@ -377,7 +377,7 @@ function Base.inv(t::AbstractTensorMap) T = float(scalartype(t)) tinv = similar(t, T, dom ← cod) for (c, b) in blocks(t) - binv = MatrixAlgebra.one!(block(tinv, c)) + binv = one!(block(tinv, c)) ldiv!(lu(b), binv) end return tinv @@ -449,11 +449,11 @@ for f in (:cos, :sin, :tan, :cot, :cosh, :sinh, :tanh, :coth, :atan, :acot, :asi tf = similar(t, T) if T <: Real for (c, b) in blocks(t) - copy!(block(tf, c), real(MatrixAlgebra.$f(b))) + copy!(block(tf, c), real(MatrixAlgebraKit.$f(b))) end else for (c, b) in blocks(t) - copy!(block(tf, c), MatrixAlgebra.$f(b)) + copy!(block(tf, c), MatrixAlgebraKit.$f(b)) end end return tf diff --git a/src/tensors/truncation.jl b/src/tensors/truncation.jl deleted file mode 100644 index e49cdc94d..000000000 --- a/src/tensors/truncation.jl +++ /dev/null @@ -1,163 +0,0 @@ -# truncation.jl -# -# Implements truncation schemes for truncating a tensor with svd, leftorth or rightorth -abstract type TruncationScheme end - -struct NoTruncation <: TruncationScheme -end -notrunc() = NoTruncation() - -struct TruncationError{T<:Real} <: TruncationScheme - ϵ::T -end -truncerr(epsilon::Real) = TruncationError(epsilon) - -struct TruncationDimension <: TruncationScheme - dim::Int -end -truncdim(d::Int) = TruncationDimension(d) - -struct TruncationSpace{S<:ElementarySpace} <: TruncationScheme - space::S -end -truncspace(space::ElementarySpace) = TruncationSpace(space) - -struct TruncationCutoff{T<:Real} <: TruncationScheme - ϵ::T - add_back::Int -end -truncbelow(epsilon::Real, add_back::Int=0) = TruncationCutoff(epsilon, add_back) - -# Compute the total truncation error given truncation dimensions -function _compute_truncerr(Σdata, truncdim, p=2) - I = keytype(Σdata) - S = scalartype(valtype(Σdata)) - return _norm((c => view(v, (truncdim[c] + 1):length(v)) for (c, v) in Σdata), p, - zero(S)) -end - -# Compute truncation dimensions -function _compute_truncdim(Σdata, ::NoTruncation, p=2) - I = keytype(Σdata) - truncdim = SectorDict{I,Int}(c => length(v) for (c, v) in Σdata) - return truncdim -end -function _compute_truncdim(Σdata, trunc::TruncationError, p=2) - I = keytype(Σdata) - S = real(eltype(valtype(Σdata))) - truncdim = SectorDict{I,Int}(c => length(Σc) for (c, Σc) in Σdata) - truncerr = zero(S) - while true - cmin = _findnexttruncvalue(Σdata, truncdim, p) - isnothing(cmin) && break - truncdim[cmin] -= 1 - truncerr = _compute_truncerr(Σdata, truncdim, p) - if truncerr > trunc.ϵ - truncdim[cmin] += 1 - break - end - end - return truncdim -end -function _compute_truncdim(Σdata, trunc::TruncationDimension, p=2) - I = keytype(Σdata) - truncdim = SectorDict{I,Int}(c => length(v) for (c, v) in Σdata) - while sum(dim(c) * d for (c, d) in truncdim) > trunc.dim - cmin = _findnexttruncvalue(Σdata, truncdim, p) - isnothing(cmin) && break - truncdim[cmin] -= 1 - end - return truncdim -end -function _compute_truncdim(Σdata, trunc::TruncationSpace, p=2) - I = keytype(Σdata) - truncdim = SectorDict{I,Int}(c => min(length(v), dim(trunc.space, c)) - for (c, v) in Σdata) - return truncdim -end - -function _compute_truncdim(Σdata, trunc::TruncationCutoff, p=2) - I = keytype(Σdata) - truncdim = SectorDict{I,Int}(c => length(v) for (c, v) in Σdata) - for (c, v) in Σdata - newdim = findlast(Base.Fix2(>, trunc.ϵ), v) - if newdim === nothing - truncdim[c] = 0 - else - truncdim[c] = newdim - end - end - for i in 1:(trunc.add_back) - cmax = _findnextgrowvalue(Σdata, truncdim, p) - isnothing(cmax) && break - truncdim[cmax] += 1 - end - return truncdim -end - -# Combine truncations -struct MultipleTruncation{T<:Tuple{Vararg{TruncationScheme}}} <: TruncationScheme - truncations::T -end -function Base.:&(a::MultipleTruncation, b::MultipleTruncation) - return MultipleTruncation((a.truncations..., b.truncations...)) -end -function Base.:&(a::MultipleTruncation, b::TruncationScheme) - return MultipleTruncation((a.truncations..., b)) -end -function Base.:&(a::TruncationScheme, b::MultipleTruncation) - return MultipleTruncation((a, b.truncations...)) -end -Base.:&(a::TruncationScheme, b::TruncationScheme) = MultipleTruncation((a, b)) - -function _compute_truncdim(Σdata, trunc::MultipleTruncation, p::Real=2) - truncdim = _compute_truncdim(Σdata, trunc.truncations[1], p) - for k in 2:length(trunc.truncations) - truncdimₖ = _compute_truncdim(Σdata, trunc.truncations[k], p) - for (c, d) in truncdim - truncdim[c] = min(d, truncdimₖ[c]) - end - end - return truncdim -end - -# auxiliary function -function _findnexttruncvalue(Σdata, truncdim::SectorDict{I,Int}, p::Real) where {I<:Sector} - # early return - (isempty(Σdata) || all(iszero, values(truncdim))) && return nothing - - # find some suitable starting candidate - cmin = findfirst(>(0), truncdim) - σmin = Σdata[cmin][truncdim[cmin]] - - # find the actual minimum singular value - for (c, σs) in Σdata - if truncdim[c] > 0 - σ = σs[truncdim[c]] - if σ < σmin - cmin, σmin = c, σ - end - end - end - return cmin -end -function _findnextgrowvalue(Σdata, truncdim::SectorDict{I,Int}, p::Real) where {I<:Sector} - istruncated = SectorDict{I,Bool}(c => (d < length(Σdata[c])) for (c, d) in truncdim) - # early return - (isempty(Σdata) || !any(values(istruncated))) && return nothing - - # find some suitable starting candidate - cmax = findfirst(istruncated) - σmax = Σdata[cmax][truncdim[cmax] + 1] - - # find the actual maximal singular value - for (c, σs) in Σdata - if istruncated[c] - σ = σs[truncdim[c] + 1] - if σ > σmax - cmax, σmax = c, σ - end - end - end - return cmax -end diff --git a/test/ad.jl b/test/ad.jl index e5e2d884d..4013846d0 100644 --- a/test/ad.jl +++ b/test/ad.jl @@ -398,7 +398,7 @@ Vlist = ((ℂ^2, (ℂ^3)', ℂ^3, ℂ^2, (ℂ^2)'), test_rrule(eigh′, H; atol, output_tangent=(ΔD, ΔU)) end - let (U, S, V, ϵ) = tsvd(A) + let (U, S, V) = tsvd(A) ΔU = randn(scalartype(U), space(U)) ΔS = randn(scalartype(S), space(S)) ΔV = randn(scalartype(V), space(V)) @@ -408,54 +408,53 @@ Vlist = ((ℂ^2, (ℂ^3)', ℂ^3, ℂ^2, (ℂ^2)'), mul!(block(ΔU, c), block(U, c), Diagonal(imag(diag(b))), -im, 1) end end - test_rrule(tsvd, A; atol, output_tangent=(ΔU, ΔS, ΔV, 0.0)) + test_rrule(tsvd, A; atol, output_tangent=(ΔU, ΔS, ΔV)) allS = mapreduce(x -> diag(x[2]), vcat, blocks(S)) truncval = (maximum(allS) + minimum(allS)) / 2 - U, S, V, ϵ = tsvd(A; trunc=truncerr(truncval)) + U, S, V = tsvd(A; trunc=truncerr(truncval)) ΔU = randn(scalartype(U), space(U)) ΔS = randn(scalartype(S), space(S)) ΔV = randn(scalartype(V), space(V)) T <: Complex && remove_svdgauge_depence!(ΔU, ΔV, U, S, V) - test_rrule(tsvd, A; atol, output_tangent=(ΔU, ΔS, ΔV, 0.0), + test_rrule(tsvd, A; atol, output_tangent=(ΔU, ΔS, ΔV), fkwargs=(; trunc=truncerr(truncval))) end - let (U, S, V, ϵ) = tsvd(B) + let (U, S, V) = tsvd(B) ΔU = randn(scalartype(U), space(U)) ΔS = randn(scalartype(S), space(S)) ΔV = randn(scalartype(V), space(V)) T <: Complex && remove_svdgauge_depence!(ΔU, ΔV, U, S, V) - test_rrule(tsvd, B; atol, output_tangent=(ΔU, ΔS, ΔV, 0.0)) + test_rrule(tsvd, B; atol, output_tangent=(ΔU, ΔS, ΔV)) Vtrunc = spacetype(S)(TensorKit.SectorDict(c => ceil(Int, size(b, 1) / 2) for (c, b) in blocks(S))) - U, S, V, ϵ = tsvd(B; trunc=truncspace(Vtrunc)) + U, S, V = tsvd(B; trunc=truncspace(Vtrunc)) ΔU = randn(scalartype(U), space(U)) ΔS = randn(scalartype(S), space(S)) ΔV = randn(scalartype(V), space(V)) T <: Complex && remove_svdgauge_depence!(ΔU, ΔV, U, S, V) - test_rrule(tsvd, B; atol, output_tangent=(ΔU, ΔS, ΔV, 0.0), + test_rrule(tsvd, B; atol, output_tangent=(ΔU, ΔS, ΔV), fkwargs=(; trunc=truncspace(Vtrunc))) end - let (U, S, V, ϵ) = tsvd(C) + let (U, S, V) = tsvd(C) ΔU = randn(scalartype(U), space(U)) ΔS = randn(scalartype(S), space(S)) ΔV = randn(scalartype(V), space(V)) T <: Complex && remove_svdgauge_depence!(ΔU, ΔV, U, S, V) - test_rrule(tsvd, C; atol, output_tangent=(ΔU, ΔS, ΔV, 0.0)) + test_rrule(tsvd, C; atol, output_tangent=(ΔU, ΔS, ΔV)) - c, = TensorKit.MatrixAlgebra._argmax(x -> sqrt(dim(x[1])) * maximum(diag(x[2])), - blocks(S)) + c, = argmax(x -> sqrt(dim(x[1])) * maximum(diag(x[2])), blocks(S)) trunc = truncdim(round(Int, 2 * dim(c))) - U, S, V, ϵ = tsvd(C; trunc) + U, S, V = tsvd(C; trunc) ΔU = randn(scalartype(U), space(U)) ΔS = randn(scalartype(S), space(S)) ΔV = randn(scalartype(V), space(V)) T <: Complex && remove_svdgauge_depence!(ΔU, ΔV, U, S, V) - test_rrule(tsvd, C; atol, output_tangent=(ΔU, ΔS, ΔV, 0.0), fkwargs=(; trunc)) + test_rrule(tsvd, C; atol, output_tangent=(ΔU, ΔS, ΔV), fkwargs=(; trunc)) end let D = LinearAlgebra.eigvals(C) diff --git a/test/diagonal.jl b/test/diagonal.jl index e5fcaa430..0f67836dd 100644 --- a/test/diagonal.jl +++ b/test/diagonal.jl @@ -200,6 +200,16 @@ diagspacelist = ((ℂ^4)', ℂ[Z2Irrep](0 => 2, 1 => 3), @test all(((s, t),) -> isapprox(s, t), zip(values(LinearAlgebra.eigvals(D)), values(LinearAlgebra.eigvals(t)))) + D, W = @constinferred eig!(t) + @test D === t + @test W == one(t) + @test t * W ≈ W * D + D2, V2 = @constinferred eigh!(t2) + if T <: Real + @test D2 === t2 + end + @test V2 == one(t) + @test t2 * V2 ≈ V2 * D2 end @testset "leftorth with $alg" for alg in (TensorKit.QR(), TensorKit.QL()) Q, R = @constinferred leftorth(t; alg=alg) diff --git a/test/factorizations.jl b/test/factorizations.jl new file mode 100644 index 000000000..c76511bd8 --- /dev/null +++ b/test/factorizations.jl @@ -0,0 +1,202 @@ +spacelist = try + if ENV["CI"] == "true" + println("Detected running on CI") + if Sys.iswindows() + (Vtr, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂) + elseif Sys.isapple() + (Vtr, Vℤ₃, VfU₁, VfSU₂) + else + (Vtr, VU₁, VCU₁, VSU₂, VfSU₂) + end + else + (Vtr, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂) + end +catch + (Vtr, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂) +end + +for V in spacelist + I = sectortype(first(V)) + Istr = TensorKit.type_repr(I) + println("---------------------------------------") + println("Tensors with symmetry: $Istr") + println("---------------------------------------") + @timedtestset "Tensors with symmetry: $Istr" verbose = true begin + V1, V2, V3, V4, V5 = V + @timedtestset "Factorization" begin + W = V1 ⊗ V2 + @testset for T in (Float32, ComplexF64) + # Test both a normal tensor and an adjoint one. + ts = (rand(T, W, W'), rand(T, W, W')') + @testset for t in ts + # test squares and rectangles here + @testset "leftorth with $alg" for alg in + (TensorKit.LAPACK_HouseholderQR(), + TensorKit.LAPACK_HouseholderQR(positive=true), + #TensorKit.QL(), + #TensorKit.QLpos(), + TensorKit.PolarViaSVD(TensorKit.LAPACK_QRIteration()), + TensorKit.PolarViaSVD(TensorKit.LAPACK_DivideAndConquer()), + TensorKit.LAPACK_QRIteration(), + TensorKit.LAPACK_DivideAndConquer()) + Q, R = @constinferred leftorth(t; alg=alg) + @test isisometry(Q) + tQR = Q * R + @test tQR ≈ t + end + @testset "leftnull with $alg" for alg in + (TensorKit.LAPACK_HouseholderQR(), + TensorKit.LAPACK_QRIteration(), + TensorKit.LAPACK_DivideAndConquer()) + N = @constinferred leftnull(t; alg=alg) + @test isisometry(N) + @test norm(N' * t) < 100 * eps(norm(t)) + end + @testset "rightorth with $alg" for alg in + (#TensorKit.RQ(), TensorKit.RQpos(), + TensorKit.LAPACK_HouseholderLQ(), + TensorKit.LAPACK_HouseholderLQ(positive=true), + TensorKit.PolarViaSVD(TensorKit.LAPACK_QRIteration()), + TensorKit.PolarViaSVD(TensorKit.LAPACK_DivideAndConquer()), + TensorKit.LAPACK_QRIteration(), + TensorKit.LAPACK_DivideAndConquer()) + L, Q = @constinferred rightorth(t; alg=alg) + @test isisometry(Q; side=:right) + @test L * Q ≈ t + end + @testset "rightnull with $alg" for alg in + (TensorKit.LAPACK_HouseholderLQ(), + TensorKit.LAPACK_QRIteration(), + TensorKit.LAPACK_DivideAndConquer()) + M = @constinferred rightnull(t; alg=alg) + @test isisometry(M; side=:right) + @test norm(t * M') < 100 * eps(norm(t)) + end + @testset "tsvd with $alg" for alg in (TensorKit.LAPACK_QRIteration(), + TensorKit.LAPACK_DivideAndConquer()) + U, S, V = @constinferred tsvd(t; alg=alg) + @test isisometry(U) + @test isisometry(V; side=:right) + @test U * S * V ≈ t + + s = LinearAlgebra.svdvals(t) + s′ = LinearAlgebra.diag(S) + for (c, b) in s + @test b ≈ s′[c] + end + s = LinearAlgebra.svdvals(t') + s′ = LinearAlgebra.diag(S') + for (c, b) in s + @test b ≈ s′[c] + end + end + @testset "cond and rank" begin + d1 = dim(codomain(t)) + d2 = dim(domain(t)) + @test rank(t) == min(d1, d2) + M = leftnull(t) + @test rank(M) == max(d1, d2) - min(d1, d2) + t3 = unitary(T, V1 ⊗ V2, V1 ⊗ V2) + @test cond(t3) ≈ one(real(T)) + @test rank(t3) == dim(V1 ⊗ V2) + t4 = randn(T, V1 ⊗ V2, V1 ⊗ V2) + t4 = (t4 + t4') / 2 + vals = LinearAlgebra.eigvals(t4) + λmax = maximum(s -> maximum(abs, s), values(vals)) + λmin = minimum(s -> minimum(abs, s), values(vals)) + @test cond(t4) ≈ λmax / λmin + vals = LinearAlgebra.eigvals(t4') + λmax = maximum(s -> maximum(abs, s), values(vals)) + λmin = minimum(s -> minimum(abs, s), values(vals)) + @test cond(t4') ≈ λmax / λmin + end + end + @testset "empty tensor" begin + t = randn(T, V1 ⊗ V2, zero(V1)) + @testset "leftorth with $alg" for alg in + (TensorKit.LAPACK_HouseholderQR(), + TensorKit.LAPACK_HouseholderQR(positive=true), + #TensorKit.QL(), TensorKit.QLpos(), + TensorKit.PolarViaSVD(TensorKit.LAPACK_QRIteration()), + TensorKit.PolarViaSVD(TensorKit.LAPACK_DivideAndConquer()), + TensorKit.LAPACK_QRIteration(), + TensorKit.LAPACK_DivideAndConquer()) + Q, R = @constinferred leftorth(t; alg=alg) + @test Q == t + @test dim(Q) == dim(R) == 0 + end + @testset "leftnull with $alg" for alg in + (TensorKit.LAPACK_HouseholderQR(), + TensorKit.LAPACK_QRIteration(), + TensorKit.LAPACK_DivideAndConquer()) + N = @constinferred leftnull(t; alg=alg) + @test isunitary(N) + end + @testset "rightorth with $alg" for alg in + (#TensorKit.RQ(), TensorKit.RQpos(), + TensorKit.LAPACK_HouseholderLQ(), + TensorKit.LAPACK_HouseholderLQ(positive=true), + TensorKit.PolarViaSVD(TensorKit.LAPACK_QRIteration()), + TensorKit.PolarViaSVD(TensorKit.LAPACK_DivideAndConquer()), + TensorKit.LAPACK_QRIteration(), + TensorKit.LAPACK_DivideAndConquer()) + L, Q = @constinferred rightorth(copy(t'); alg=alg) + @test Q == t' + @test dim(Q) == dim(L) == 0 + end + @testset "rightnull with $alg" for alg in + (TensorKit.LAPACK_HouseholderLQ(), + TensorKit.LAPACK_QRIteration(), + TensorKit.LAPACK_DivideAndConquer()) + M = @constinferred rightnull(copy(t'); alg=alg) + @test isunitary(M) + end + @testset "tsvd with $alg" for alg in (TensorKit.LAPACK_QRIteration(), + TensorKit.LAPACK_DivideAndConquer()) + U, S, V = @constinferred tsvd(t; alg=alg) + @test U == t + @test dim(U) == dim(S) == dim(V) + end + @testset "cond and rank" begin + @test rank(t) == 0 + W2 = zero(V1) * zero(V2) + t2 = rand(W2, W2) + @test rank(t2) == 0 + @test cond(t2) == 0.0 + end + end + @testset "eig and isposdef" begin + t = rand(T, V1, V1) + D, V = eigen(t) + @test t * V ≈ V * D + + d = LinearAlgebra.eigvals(t; sortby=nothing) + d′ = LinearAlgebra.diag(D) + for (c, b) in d + @test b ≈ d′[c] + end + + # Somehow moving these test before the previous one gives rise to errors + # with T=Float32 on x86 platforms. Is this an OpenBLAS issue? + VdV = V' * V + VdV = (VdV + VdV') / 2 + @test isposdef(VdV) + + @test !isposdef(t) # unlikely for non-hermitian map + t2 = (t + t') + D, V = eigen(t2) + @test isisometry(V) + D̃, Ṽ = @constinferred eigh(t2) + @test D ≈ D̃ + @test V ≈ Ṽ + λ = minimum(minimum(real(LinearAlgebra.diag(b))) + for (c, b) in blocks(D)) + @test cond(Ṽ) ≈ one(real(T)) + @test isposdef(t2) == isposdef(λ) + @test isposdef(t2 - λ * one(t2) + 0.1 * one(t2)) + @test !isposdef(t2 - λ * one(t2) - 0.1 * one(t2)) + end + end + end + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 9e44b8a3a..a1b9963b1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -113,15 +113,16 @@ VSU₂U₁ = (Vect[SU2Irrep ⊠ U1Irrep]((0, 0) => 1, (1 // 2, -1) => 1), # ℂ[SU3Irrep]((0, 0, 0) => 1, (1, 0, 0) => 1, (1, 1, 0) => 1)') Ti = time() -include("fusiontrees.jl") -include("spaces.jl") -include("tensors.jl") -include("diagonal.jl") -include("planar.jl") +#include("fusiontrees.jl") +#include("spaces.jl") +#include("tensors.jl") +include("factorizations.jl") +#include("diagonal.jl") +#include("planar.jl") # TODO: remove once we know AD is slow on macOS CI -if !(Sys.isapple() && get(ENV, "CI", "false") == "true") && isempty(VERSION.prerelease) +#=if !(Sys.isapple() && get(ENV, "CI", "false") == "true") && isempty(VERSION.prerelease) include("ad.jl") -end +end=# include("bugfixes.jl") Tf = time() printstyled("Finished all tests in ", diff --git a/test/tensors.jl b/test/tensors.jl index 25bc157b9..d5ca54139 100644 --- a/test/tensors.jl +++ b/test/tensors.jl @@ -1,9 +1,3 @@ -for V in (Vtr, Vℤ₂, Vfℤ₂, Vℤ₃, VU₁, VfU₁, VCU₁, VSU₂, VfSU₂, VSU₂U₁)#, VSU₃) - V1, V2, V3, V4, V5 = V - @assert V3 * V4 * V2 ≿ V1' * V5' # necessary for leftorth tests - @assert V3 * V4 ≾ V1' * V2' * V5' # necessary for rightorth tests -end - spacelist = try if ENV["CI"] == "true" println("Detected running on CI") @@ -369,9 +363,8 @@ for V in spacelist for T in (Float64, ComplexF64) t1 = randisometry(T, W1, W2) t2 = randisometry(T, W2 ← W2) - @test t1' * t1 ≈ one(t2) - @test t2' * t2 ≈ one(t2) - @test t2 * t2' ≈ one(t2) + @test isisometry(t1) + @test isunitary(t2) P = t1 * t1' @test P * P ≈ P end @@ -439,173 +432,6 @@ for V in spacelist @test LinearAlgebra.isdiag(D) @test LinearAlgebra.diag(D) == d end - @timedtestset "Factorization" begin - W = V1 ⊗ V2 ⊗ V3 ⊗ V4 ⊗ V5 - for T in (Float32, ComplexF64) - # Test both a normal tensor and an adjoint one. - ts = (rand(T, W), rand(T, W)') - for t in ts - @testset "leftorth with $alg" for alg in - (TensorKit.QR(), TensorKit.QRpos(), - TensorKit.QL(), TensorKit.QLpos(), - TensorKit.Polar(), TensorKit.SVD(), - TensorKit.SDD()) - Q, R = @constinferred leftorth(t, ((3, 4, 2), (1, 5)); alg=alg) - QdQ = Q' * Q - @test QdQ ≈ one(QdQ) - @test Q * R ≈ permute(t, ((3, 4, 2), (1, 5))) - if alg isa Polar - @test isposdef(R) - @test domain(R) == codomain(R) == space(t, 1)' ⊗ space(t, 5)' - end - end - @testset "leftnull with $alg" for alg in - (TensorKit.QR(), TensorKit.SVD(), - TensorKit.SDD()) - N = @constinferred leftnull(t, ((3, 4, 2), (1, 5)); alg=alg) - NdN = N' * N - @test NdN ≈ one(NdN) - @test norm(N' * permute(t, ((3, 4, 2), (1, 5)))) < - 100 * eps(norm(t)) - end - @testset "rightorth with $alg" for alg in - (TensorKit.RQ(), TensorKit.RQpos(), - TensorKit.LQ(), TensorKit.LQpos(), - TensorKit.Polar(), TensorKit.SVD(), - TensorKit.SDD()) - L, Q = @constinferred rightorth(t, ((3, 4), (2, 1, 5)); alg=alg) - QQd = Q * Q' - @test QQd ≈ one(QQd) - @test L * Q ≈ permute(t, ((3, 4), (2, 1, 5))) - if alg isa Polar - @test isposdef(L) - @test domain(L) == codomain(L) == space(t, 3) ⊗ space(t, 4) - end - end - @testset "rightnull with $alg" for alg in - (TensorKit.LQ(), TensorKit.SVD(), - TensorKit.SDD()) - M = @constinferred rightnull(t, ((3, 4), (2, 1, 5)); alg=alg) - MMd = M * M' - @test MMd ≈ one(MMd) - @test norm(permute(t, ((3, 4), (2, 1, 5))) * M') < - 100 * eps(norm(t)) - end - @testset "tsvd with $alg" for alg in (TensorKit.SVD(), TensorKit.SDD()) - U, S, V = @constinferred tsvd(t, ((3, 4, 2), (1, 5)); alg=alg) - UdU = U' * U - @test UdU ≈ one(UdU) - VVd = V * V' - @test VVd ≈ one(VVd) - t2 = permute(t, ((3, 4, 2), (1, 5))) - @test U * S * V ≈ t2 - - s = LinearAlgebra.svdvals(t2) - s′ = LinearAlgebra.diag(S) - for (c, b) in s - @test b ≈ s′[c] - end - end - @testset "cond and rank" begin - t2 = permute(t, ((3, 4, 2), (1, 5))) - d1 = dim(codomain(t2)) - d2 = dim(domain(t2)) - @test rank(t2) == min(d1, d2) - M = leftnull(t2) - @test rank(M) == max(d1, d2) - min(d1, d2) - t3 = unitary(T, V1 ⊗ V2, V1 ⊗ V2) - @test cond(t3) ≈ one(real(T)) - @test rank(t3) == dim(V1 ⊗ V2) - t4 = randn(T, V1 ⊗ V2, V1 ⊗ V2) - t4 = (t4 + t4') / 2 - vals = LinearAlgebra.eigvals(t4) - λmax = maximum(s -> maximum(abs, s), values(vals)) - λmin = minimum(s -> minimum(abs, s), values(vals)) - @test cond(t4) ≈ λmax / λmin - end - end - @testset "empty tensor" begin - t = randn(T, V1 ⊗ V2, zero(V1)) - @testset "leftorth with $alg" for alg in - (TensorKit.QR(), TensorKit.QRpos(), - TensorKit.QL(), TensorKit.QLpos(), - TensorKit.Polar(), TensorKit.SVD(), - TensorKit.SDD()) - Q, R = @constinferred leftorth(t; alg=alg) - @test Q == t - @test dim(Q) == dim(R) == 0 - end - @testset "leftnull with $alg" for alg in - (TensorKit.QR(), TensorKit.SVD(), - TensorKit.SDD()) - N = @constinferred leftnull(t; alg=alg) - @test N' * N ≈ id(domain(N)) - @test N * N' ≈ id(codomain(N)) - end - @testset "rightorth with $alg" for alg in - (TensorKit.RQ(), TensorKit.RQpos(), - TensorKit.LQ(), TensorKit.LQpos(), - TensorKit.Polar(), TensorKit.SVD(), - TensorKit.SDD()) - L, Q = @constinferred rightorth(copy(t'); alg=alg) - @test Q == t' - @test dim(Q) == dim(L) == 0 - end - @testset "rightnull with $alg" for alg in - (TensorKit.LQ(), TensorKit.SVD(), - TensorKit.SDD()) - M = @constinferred rightnull(copy(t'); alg=alg) - @test M * M' ≈ id(codomain(M)) - @test M' * M ≈ id(domain(M)) - end - @testset "tsvd with $alg" for alg in (TensorKit.SVD(), TensorKit.SDD()) - U, S, V = @constinferred tsvd(t; alg=alg) - @test U == t - @test dim(U) == dim(S) == dim(V) - end - @testset "cond and rank" begin - @test rank(t) == 0 - W2 = zero(V1) * zero(V2) - t2 = rand(W2, W2) - @test rank(t2) == 0 - @test cond(t2) == 0.0 - end - end - t = rand(T, V1 ⊗ V1' ⊗ V2 ⊗ V2') - @testset "eig and isposdef" begin - D, V = eigen(t, ((1, 3), (2, 4))) - t2 = permute(t, ((1, 3), (2, 4))) - @test t2 * V ≈ V * D - - d = LinearAlgebra.eigvals(t2; sortby=nothing) - d′ = LinearAlgebra.diag(D) - for (c, b) in d - @test b ≈ d′[c] - end - - # Somehow moving these test before the previous one gives rise to errors - # with T=Float32 on x86 platforms. Is this an OpenBLAS issue? - VdV = V' * V - VdV = (VdV + VdV') / 2 - @test isposdef(VdV) - - @test !isposdef(t2) # unlikely for non-hermitian map - t2 = (t2 + t2') - D, V = eigen(t2) - VdV = V' * V - @test VdV ≈ one(VdV) - D̃, Ṽ = @constinferred eigh(t2) - @test D ≈ D̃ - @test V ≈ Ṽ - λ = minimum(minimum(real(LinearAlgebra.diag(b))) - for (c, b) in blocks(D)) - @test cond(Ṽ) ≈ one(real(T)) - @test isposdef(t2) == isposdef(λ) - @test isposdef(t2 - λ * one(t2) + 0.1 * one(t2)) - @test !isposdef(t2 - λ * one(t2) - 0.1 * one(t2)) - end - end - end @timedtestset "Tensor truncation" begin for T in (Float32, ComplexF64) for p in (1, 2, 3, Inf) @@ -615,23 +441,36 @@ for V in spacelist for t in ts U₀, S₀, V₀, = tsvd(t) t = rmul!(t, 1 / norm(S₀, p)) - U, S, V, ϵ = @constinferred tsvd(t; trunc=truncerr(5e-1), p=p) + U, S, V = @constinferred tsvd(t; trunc=truncerr(5e-1, p)) + ϵ = TensorKit._norm(LinearAlgebra.svdvals(U * S * V - t), p, + zero(scalartype(S))) + p == 2 && @test ϵ < 5e-1 # @show p, ϵ # @show domain(S) # @test min(space(S,1), space(S₀,1)) != space(S₀,1) - U′, S′, V′, ϵ′ = tsvd(t; trunc=truncerr(nextfloat(ϵ)), p=p) + U′, S′, V′ = tsvd(t; trunc=truncerr(ϵ + 10eps(ϵ), p)) + ϵ′ = TensorKit._norm(LinearAlgebra.svdvals(U′ * S′ * V′ - t), p, + zero(scalartype(S))) + @test (U, S, V, ϵ) == (U′, S′, V′, ϵ′) - U′, S′, V′, ϵ′ = tsvd(t; trunc=truncdim(ceil(Int, dim(domain(S)))), - p=p) + U′, S′, V′ = tsvd(t; trunc=truncdim(ceil(Int, dim(domain(S))))) + ϵ′ = TensorKit._norm(LinearAlgebra.svdvals(U′ * S′ * V′ - t), p, + zero(scalartype(S))) @test (U, S, V, ϵ) == (U′, S′, V′, ϵ′) - U′, S′, V′, ϵ′ = tsvd(t; trunc=truncspace(space(S, 1)), p=p) + U′, S′, V′ = tsvd(t; trunc=truncspace(space(S, 1))) + ϵ′ = TensorKit._norm(LinearAlgebra.svdvals(U′ * S′ * V′ - t), p, + zero(scalartype(S))) @test (U, S, V, ϵ) == (U′, S′, V′, ϵ′) # results with truncationcutoff cannot be compared because they don't take degeneracy into account, and thus truncate differently - U, S, V, ϵ = tsvd(t; trunc=truncbelow(1 / dim(domain(S₀))), p=p) + U, S, V = tsvd(t; trunc=truncbelow(1 / dim(domain(S₀)))) + ϵ = TensorKit._norm(LinearAlgebra.svdvals(U * S * V - t), p, + zero(scalartype(S))) # @show p, ϵ # @show domain(S) # @test min(space(S,1), space(S₀,1)) != space(S₀,1) - U′, S′, V′, ϵ′ = tsvd(t; trunc=truncspace(space(S, 1)), p=p) + U′, S′, V′ = tsvd(t; trunc=truncspace(space(S, 1))) + ϵ′ = TensorKit._norm(LinearAlgebra.svdvals(U′ * S′ * V′ - t), p, + zero(scalartype(S))) @test (U, S, V, ϵ) == (U′, S′, V′, ϵ′) end end @@ -691,8 +530,8 @@ for V in spacelist for T in (Float32, ComplexF64) tA = rand(T, V1 ⊗ V3, V1 ⊗ V3) tB = rand(T, V2 ⊗ V4, V2 ⊗ V4) - tA = 3 // 2 * leftorth(tA; alg=Polar())[1] - tB = 1 // 5 * leftorth(tB; alg=Polar())[1] + tA = 3 // 2 * left_polar(tA)[1] + tB = 1 // 5 * left_polar(tB)[1] tC = rand(T, V1 ⊗ V3, V2 ⊗ V4) t = @constinferred sylvester(tA, tB, tC) @test codomain(t) == V1 ⊗ V3