From 836836d273bfacfcd449a794245f59360fe7695b Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 12 Mar 2025 14:27:37 -0400 Subject: [PATCH 01/26] Add `qr` and `svd` --- Project.toml | 2 + src/factorizations.jl | 100 +++++++++++++++++++++++------------- test/Project.toml | 1 + test/test_basics.jl | 23 --------- test/test_factorizations.jl | 90 ++++++++++++++++++++++++++++++++ 5 files changed, 156 insertions(+), 60 deletions(-) create mode 100644 test/test_factorizations.jl diff --git a/Project.toml b/Project.toml index 9fcc3ac..1c2ade8 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4" TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6" TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" @@ -23,6 +24,7 @@ BlockArrays = "1.2.0" EllipsisNotation = "1.8.0" GradedUnitRanges = "0.1.0" LinearAlgebra = "1.10" +MatrixAlgebraKit = "0.1.0" TupleTools = "1.6.0" TypeParameterAccessors = "0.2.1, 0.3" julia = "1.10" diff --git a/src/factorizations.jl b/src/factorizations.jl index a017ca1..f482cb3 100644 --- a/src/factorizations.jl +++ b/src/factorizations.jl @@ -1,45 +1,71 @@ -using ArrayLayouts: LayoutMatrix -using LinearAlgebra: LinearAlgebra, Diagonal +using MatrixAlgebraKit: qr_full, qr_compact, svd_full, svd_compact, svd_trunc -function qr(a::AbstractArray, biperm::BlockedPermutation{2}) - a_matricized = fusedims(a, biperm) - # TODO: Make this more generic, allow choosing thin or full, - # make sure this works on GPU. - q_fact, r_matricized = LinearAlgebra.qr(a_matricized) - q_matricized = typeof(a_matricized)(q_fact) - axes_codomain, axes_domain = blockpermute(axes(a), biperm) - axes_q = (axes_codomain..., axes(q_matricized, 2)) - axes_r = (axes(r_matricized, 1), axes_domain...) - q = splitdims(q_matricized, axes_q) - r = splitdims(r_matricized, axes_r) - return q, r +# TODO: consider in-place version +# TODO: figure out kwargs and document +# +""" + qr(A::AbstractArray, labels_A, labels_codomain, labels_domain; full=true, kwargs...) -> Q, R + qr(A::AbstractArray, biperm::BlockedPermutation{2}; full=true, kwargs...) -> Q, R +Compute the QR decomposition of a generic N-dimensional array, by interpreting it as +a linear map from the domain to the codomain indices. These can be specified either via +their labels, or directly through a `biperm`. +""" +function qr(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) + biperm = blockedperm_indexin(Tuple.((labels_A, labels_codomain, labels_domain))...) + return qr(A, biperm) end +function qr(A::AbstractArray, biperm::BlockedPermutation{2}; full::Bool=true, kwargs...) + # tensor to matrix + A_mat = fusedims(A, biperm) -function qr(a::AbstractArray, labels_a, labels_codomain, labels_domain) - # TODO: Generalize to conversion to `Tuple` isn't needed. - return qr( - a, blockedperm_indexin(Tuple(labels_a), Tuple(labels_codomain), Tuple(labels_domain)) - ) + # factorization + Q, R = full ? qr_full(A_mat; kwargs...) : qr_compact(A_mat; kwargs...) + + # matrix to tensor + axes_codomain, axes_domain = blockpermute(axes(A), biperm) + axes_Q = (axes_codomain..., axes(Q, 2)) + axes_R = (axes(R, 1), axes_domain...) + return splitdims(Q, axes_Q), splitdims(R, axes_R) end -function svd(a::AbstractArray, biperm::BlockedPermutation{2}) - a_matricized = fusedims(a, biperm) - usv_matricized = LinearAlgebra.svd(a_matricized) - u_matricized = usv_matricized.U - s_diag = usv_matricized.S - v_matricized = usv_matricized.Vt - axes_codomain, axes_domain = blockpermute(axes(a), biperm) - axes_u = (axes_codomain..., axes(u_matricized, 2)) - axes_v = (axes(v_matricized, 1), axes_domain...) - u = splitdims(u_matricized, axes_u) - # TODO: Use `DiagonalArrays.diagonal` to make it more general. - s = Diagonal(s_diag) - v = splitdims(v_matricized, axes_v) - return u, s, v +# TODO: separate out the algorithm selection step from the implementation +""" + svd(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) -> U, S, Vᴴ + svd(A::AbstractArray, biperm::BlockedPermutation{2}; kwargs...) -> U, S, Vᴴ +Compute the SVD decomposition of a generic N-dimensional array, by interpreting it as +a linear map from the domain to the codomain indices. These can be specified either via +their labels, or directly through a `biperm`. +## Keyword arguments +- `full::Bool=false`: select between a "thick" or a "thin" decomposition, where both `U` and `Vᴴ` + are unitary or isometric. +- `trunc`: Truncation keywords for `svd_trunc`. Not compatible with `full=true`. +- Other keywords are passed on directly to MatrixAlgebraKit +""" +function svd(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) + biperm = blockedperm_indexin(Tuple.((labels_A, labels_codomain, labels_domain))...) + return svd(A, biperm; kwargs...) end +function svd( + A::AbstractArray, + biperm::BlockedPermutation{2}; + full::Bool=false, + trunc=nothing, + kwargs..., +) + # tensor to matrix + A_mat = fusedims(A, biperm) + + # factorization + if !isnothing(trunc) + @assert !full "Specified both full and truncation, currently not supported" + U, S, Vᴴ = svd_trunc(A_mat; trunc, kwargs...) + else + U, S, Vᴴ = full ? svd_full(A_mat; kwargs...) : svd_compact(A_mat; kwargs...) + end -function svd(a::AbstractArray, labels_a, labels_codomain, labels_domain) - return svd( - a, blockedperm_indexin(Tuple(labels_a), Tuple(labels_codomain), Tuple(labels_domain)) - ) + # matrix to tensor + axes_codomain, axes_domain = blockpermute(axes(A), biperm) + axes_U = (axes_codomain..., axes(U, 2)) + axes_Vᴴ = (axes(Vᴴ, 1), axes_domain...) + return splitdims(U, axes_U), S, splitdims(Vᴴ, axes_Vᴴ) end diff --git a/test/Project.toml b/test/Project.toml index 7258ca4..97ff043 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -6,6 +6,7 @@ GradedUnitRanges = "e2de450a-8a67-46c7-b59c-01d5a3d041c5" JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb" LabelledNumbers = "f856a3a6-4152-4ec4-b2a7-02c1a55d7993" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" diff --git a/test/test_basics.jl b/test/test_basics.jl index 6b75c69..3bce1b7 100644 --- a/test/test_basics.jl +++ b/test/test_basics.jl @@ -212,26 +212,3 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64}) @test a_dest[] ≈ s[] * t[] end end -@testset "qr (eltype=$elt)" for elt in elts - a = randn(elt, 5, 4, 3, 2) - labels_a = (:a, :b, :c, :d) - labels_q = (:b, :a) - labels_r = (:d, :c) - q, r = qr(a, labels_a, labels_q, labels_r) - label_qr = :qr - a′ = contract(labels_a, q, (labels_q..., label_qr), r, (label_qr, labels_r...)) - @test a ≈ a′ -end -@testset "svd (eltype=$elt)" for elt in elts - a = randn(elt, 5, 4, 3, 2) - labels_a = (:a, :b, :c, :d) - labels_u = (:b, :a) - labels_v = (:d, :c) - u, s, v = svd(a, labels_a, labels_u, labels_v) - label_u = :u - label_v = :v - # TODO: Define multi-arg `contract`? - us, labels_us = contract(u, (labels_u..., label_u), s, (label_u, label_v)) - a′ = contract(labels_a, us, labels_us, v, (label_v, labels_v...)) - @test a ≈ a′ -end diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl new file mode 100644 index 0000000..06cc18c --- /dev/null +++ b/test/test_factorizations.jl @@ -0,0 +1,90 @@ +using Test: @test, @testset, @inferred +using TestExtras: @constinferred +using TensorAlgebra: contract, qr, svd, TensorAlgebra +using MatrixAlgebraKit: truncrank + +elts = (Float64, ComplexF64) + +# QR Decomposition +# ---------------- +@testset "Full QR ($T)" for T in elts + A = randn(T, 5, 4, 3, 2) + labels_A = (:a, :b, :c, :d) + labels_Q = (:b, :a) + labels_R = (:d, :c) + + Acopy = deepcopy(A) + Q, R = @constinferred qr(A, labels_A, labels_Q, labels_R; full=true) + @test A == Acopy # should not have altered initial array + A′ = contract(labels_A, Q, (labels_Q..., :q), R, (:q, labels_R...)) + @test A ≈ A′ + @test size(Q, 1) * size(Q, 2) == size(Q, 3) # Q is unitary +end + +@testset "Compact QR ($T)" for T in elts + A = randn(T, 2, 3, 4, 5) # compact only makes a difference for less columns + labels_A = (:a, :b, :c, :d) + labels_Q = (:b, :a) + labels_R = (:d, :c) + + Acopy = deepcopy(A) + Q, R = @constinferred qr(A, labels_A, labels_Q, labels_R; full=false) + @test A == Acopy # should not have altered initial array + A′ = contract(labels_A, Q, (labels_Q..., :q), R, (:q, labels_R...)) + @test A ≈ A′ + @test size(Q, 3) == min(size(A, 1) * size(A, 2), size(A, 3) * size(A, 4)) +end + +# Singular Value Decomposition +# ---------------------------- +@testset "Full SVD ($T)" for T in elts + A = randn(T, 5, 4, 3, 2) + labels_A = (:a, :b, :c, :d) + labels_U = (:b, :a) + labels_Vᴴ = (:d, :c) + + Acopy = deepcopy(A) + U, S, Vᴴ = @constinferred svd(A, labels_A, labels_U, labels_Vᴴ; full=true) + @test A == Acopy # should not have altered initial array + US, labels_US = contract(U, (labels_U..., :u), S, (:u, :v)) + A′ = contract(labels_A, US, labels_US, Vᴴ, (:v, labels_Vᴴ...)) + @test A ≈ A′ + @test size(U, 1) * size(U, 2) == size(U, 3) # U is unitary + @test size(Vᴴ, 1) == size(Vᴴ, 2) * size(Vᴴ, 3) # V is unitary +end + +@testset "Compact SVD ($T)" for T in elts + A = randn(T, 5, 4, 3, 2) + labels_A = (:a, :b, :c, :d) + labels_U = (:b, :a) + labels_Vᴴ = (:d, :c) + + Acopy = deepcopy(A) + U, S, Vᴴ = @constinferred svd(A, labels_A, labels_U, labels_Vᴴ; full=false) + @test A == Acopy # should not have altered initial array + US, labels_US = contract(U, (labels_U..., :u), S, (:u, :v)) + A′ = contract(labels_A, US, labels_US, Vᴴ, (:v, labels_Vᴴ...)) + @test A ≈ A′ + k = min(size(S)...) + @test size(U, 3) == k == size(Vᴴ, 1) +end + +@testset "Truncated SVD ($T)" for T in elts + A = randn(T, 5, 4, 3, 2) + labels_A = (:a, :b, :c, :d) + labels_U = (:b, :a) + labels_Vᴴ = (:d, :c) + + # test truncated SVD + Acopy = deepcopy(A) + _, S_untrunc, _ = svd(A, labels_A, labels_U, labels_Vᴴ) + + trunc = truncrank(size(S_untrunc, 1) - 1) + U, S, Vᴴ = @constinferred svd(A, labels_A, labels_U, labels_Vᴴ; trunc) + + @test A == Acopy # should not have altered initial array + US, labels_US = contract(U, (labels_U..., :u), S, (:u, :v)) + A′ = contract(labels_A, US, labels_US, Vᴴ, (:v, labels_Vᴴ...)) + @test norm(A - A′) ≈ S_untrunc[end] + @test size(S, 1) == size(S_untrunc, 1) - 1 +end From 4a81a4b0742ab19b7e331849752d026181a1f6f4 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 12 Mar 2025 14:42:35 -0400 Subject: [PATCH 02/26] Fix test dependency --- test/test_factorizations.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl index 06cc18c..dc3fbc9 100644 --- a/test/test_factorizations.jl +++ b/test/test_factorizations.jl @@ -2,6 +2,7 @@ using Test: @test, @testset, @inferred using TestExtras: @constinferred using TensorAlgebra: contract, qr, svd, TensorAlgebra using MatrixAlgebraKit: truncrank +using LinearAlgebra: norm elts = (Float64, ComplexF64) From c8df746e51a4e37b034511edd19b3e4d51bf24fa Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 12 Mar 2025 14:42:40 -0400 Subject: [PATCH 03/26] Add `eig` --- src/factorizations.jl | 78 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 77 insertions(+), 1 deletion(-) diff --git a/src/factorizations.jl b/src/factorizations.jl index f482cb3..0cb7cfd 100644 --- a/src/factorizations.jl +++ b/src/factorizations.jl @@ -1,4 +1,5 @@ -using MatrixAlgebraKit: qr_full, qr_compact, svd_full, svd_compact, svd_trunc +using MatrixAlgebraKit: + lq_full, lq_compact, qr_full, qr_compact, svd_full, svd_compact, svd_trunc # TODO: consider in-place version # TODO: figure out kwargs and document @@ -6,6 +7,7 @@ using MatrixAlgebraKit: qr_full, qr_compact, svd_full, svd_compact, svd_trunc """ qr(A::AbstractArray, labels_A, labels_codomain, labels_domain; full=true, kwargs...) -> Q, R qr(A::AbstractArray, biperm::BlockedPermutation{2}; full=true, kwargs...) -> Q, R + Compute the QR decomposition of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain indices. These can be specified either via their labels, or directly through a `biperm`. @@ -28,14 +30,88 @@ function qr(A::AbstractArray, biperm::BlockedPermutation{2}; full::Bool=true, kw return splitdims(Q, axes_Q), splitdims(R, axes_R) end +""" + lq(A::AbstractArray, labels_A, labels_codomain, labels_domain; full=true, kwargs...) -> L, Q + lq(A::AbstractArray, biperm::BlockedPermutation{2}; full=true, kwargs...) -> L, Q + +Compute the LQ decomposition of a generic N-dimensional array, by interpreting it as +a linear map from the domain to the codomain indices. These can be specified either via +their labels, or directly through a `biperm`. +""" +function lq(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) + biperm = blockedperm_indexin(Tuple.((labels_A, labels_codomain, labels_domain))...) + return lq(A, biperm) +end +function lq(A::AbstractArray, biperm::BlockedPermutation{2}; full::Bool=true, kwargs...) + # tensor to matrix + A_mat = fusedims(A, biperm) + + # factorization + L, Q = full ? lq_full(A_mat; kwargs...) : lq_compact(A_mat; kwargs...) + + # matrix to tensor + axes_codomain, axes_domain = blockpermute(axes(A), biperm) + axes_L = (axes_codomain..., axes(L, 2)) + axes_Q = (axes(Q, 1), axes_domain...) + return splitdims(L, axes_L), splitdims(Q, axes_Q) +end + +# TODO: what name do we want? +""" + eig(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) -> D, V + eig(A::AbstractArray, biperm::BlockedPermutation{2}; kwargs...) -> D, V + +Compute the eigenvalue decomposition of a generic N-dimensional array, by interpreting it as +a linear map from the domain to the codomain indices. These can be specified either via +their labels, or directly through a `biperm`. + +## Keyword arguments + +- `ishermitian::Bool`: specify if the matrix is Hermitian, which can be used to speed up the + computation. If `false`, the output `eltype` will always be `<:Complex`. +- `trunc`: Truncation keywords for `eig(h)_trunc`. +- Other keywords are passed on directly to MatrixAlgebraKit +""" +function eig(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) + biperm = blockedperm_indexin(Tuple.((labels_A, labels_codomain, labels_domain))...) + return eig(A, biperm; kwargs...) +end +function eig( + A::AbstractArray, + biperm::BlockedPermutation{2}; + trunc=nothing, + ishermitian=nothing, + kwargs..., +) + # tensor to matrix + A_mat = fusedims(A, biperm) + + ishermitian = @something ishermitian LinearAlgebra.ishermitian(A_mat) + + # factorization + if !isnothing(trunc) + D, V = (ishermitian ? eigh_trunc : eig_trunc)(A_mat; trunc, kwargs...) + else + D, V = (ishermitian ? eigh_full : eig_full)(A_mat; kwargs...) + end + + # matrix to tensor + axes_codomain, = blockpermute(axes(A), biperm) + axes_V = (axes_codomain..., axes(V, ndims(V))) + return D, splitdims(V, axes_V) +end + # TODO: separate out the algorithm selection step from the implementation """ svd(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) -> U, S, Vᴴ svd(A::AbstractArray, biperm::BlockedPermutation{2}; kwargs...) -> U, S, Vᴴ + Compute the SVD decomposition of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain indices. These can be specified either via their labels, or directly through a `biperm`. + ## Keyword arguments + - `full::Bool=false`: select between a "thick" or a "thin" decomposition, where both `U` and `Vᴴ` are unitary or isometric. - `trunc`: Truncation keywords for `svd_trunc`. Not compatible with `full=true`. From 36d3c110f662b1bb15c8e6cba06c8231205c8f31 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 12 Mar 2025 15:06:42 -0400 Subject: [PATCH 04/26] Export factorizations --- src/TensorAlgebra.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TensorAlgebra.jl b/src/TensorAlgebra.jl index caa2cc5..ad8d3b6 100644 --- a/src/TensorAlgebra.jl +++ b/src/TensorAlgebra.jl @@ -1,6 +1,6 @@ module TensorAlgebra -export contract, contract! +export contract, contract!, eig, lq, qr, svd include("blockedtuple.jl") include("blockedpermutation.jl") From 39f96ea51075d116a254956fc5377bc440b1ec61 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 12 Mar 2025 15:07:07 -0400 Subject: [PATCH 05/26] import eig(h)_full/trunc etc --- src/factorizations.jl | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/factorizations.jl b/src/factorizations.jl index 0cb7cfd..e346ba5 100644 --- a/src/factorizations.jl +++ b/src/factorizations.jl @@ -1,5 +1,16 @@ using MatrixAlgebraKit: - lq_full, lq_compact, qr_full, qr_compact, svd_full, svd_compact, svd_trunc + eig_full, + eigh_full, + eig_trunc, + eigh_trunc, + lq_full, + lq_compact, + qr_full, + qr_compact, + svd_full, + svd_compact, + svd_trunc +using LinearAlgebra: LinearAlgebra # TODO: consider in-place version # TODO: figure out kwargs and document From 4ba57a7b90b930dda0cea27ad68225bac350f813 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 12 Mar 2025 15:07:18 -0400 Subject: [PATCH 06/26] Add tests eigenvalue decompositions --- test/test_factorizations.jl | 40 ++++++++++++++++++++++++++++++++++++- 1 file changed, 39 insertions(+), 1 deletion(-) diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl index dc3fbc9..fe9d708 100644 --- a/test/test_factorizations.jl +++ b/test/test_factorizations.jl @@ -1,6 +1,6 @@ using Test: @test, @testset, @inferred using TestExtras: @constinferred -using TensorAlgebra: contract, qr, svd, TensorAlgebra +using TensorAlgebra: contract, qr, svd, TensorAlgebra, eig using MatrixAlgebraKit: truncrank using LinearAlgebra: norm @@ -36,6 +36,44 @@ end @test size(Q, 3) == min(size(A, 1) * size(A, 2), size(A, 3) * size(A, 4)) end +# Eigenvalue Decomposition +# ------------------------ +@testset "Eigenvalue decomposition ($T)" for T in elts + A = randn(T, 4, 3, 4, 3) # needs to be square + labels_A = (:a, :b, :c, :d) + labels_V = (:b, :a) + labels_V′ = (:d, :c) + + Acopy = deepcopy(A) + # type-unstable if `ishermitian` not set + D, V = @constinferred eig(A, labels_A, labels_V, labels_V′; ishermitian=false) + @test A == Acopy # should not have altered initial array + @test eltype(D) == eltype(V) && eltype(D) <: Complex + + AV = contract((:a, :b, :D), A, labels_A, V, (labels_V′..., :D)) + VD = contract((:a, :b, :D), V, (labels_V..., :D′), D, (:D′, :D)) + @test AV ≈ VD +end + +@testset "Hermitian eigenvalue decomposition ($T)" for T in elts + A = randn(T, 12, 12) + A = reshape(A + A', 4, 3, 4, 3) + labels_A = (:a, :b, :c, :d) + labels_V = (:b, :a) + labels_V′ = (:d, :c) + + Acopy = deepcopy(A) + # type-unstable if `ishermitian` not set + D, V = @constinferred eig(A, labels_A, labels_V, labels_V′; ishermitian=true) + @test A == Acopy # should not have altered initial array + @test eltype(D) <: Real + @test eltype(V) == eltype(A) + + AV = contract((:a, :b, :D), A, labels_A, V, (labels_V′..., :D)) + VD = contract((:a, :b, :D), V, (labels_V..., :D′), D, (:D′, :D)) + @test AV ≈ VD +end + # Singular Value Decomposition # ---------------------------- @testset "Full SVD ($T)" for T in elts From e8e56e79342128628de2b3e2814cc0b582931f89 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 12 Mar 2025 15:17:28 -0400 Subject: [PATCH 07/26] Add `svdvals` and `eigvals` --- src/TensorAlgebra.jl | 2 +- src/factorizations.jl | 41 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 42 insertions(+), 1 deletion(-) diff --git a/src/TensorAlgebra.jl b/src/TensorAlgebra.jl index ad8d3b6..5f48988 100644 --- a/src/TensorAlgebra.jl +++ b/src/TensorAlgebra.jl @@ -1,6 +1,6 @@ module TensorAlgebra -export contract, contract!, eig, lq, qr, svd +export contract, contract!, eig, eigvals, lq, qr, svd, svdvals include("blockedtuple.jl") include("blockedpermutation.jl") diff --git a/src/factorizations.jl b/src/factorizations.jl index e346ba5..58a65b9 100644 --- a/src/factorizations.jl +++ b/src/factorizations.jl @@ -112,6 +112,30 @@ function eig( return D, splitdims(V, axes_V) end +""" + eigvals(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) -> D + eigvals(A::AbstractArray, biperm::BlockedPermutation{2}; kwargs...) -> D + +Compute the eigenvalues of a generic N-dimensional array, by interpreting it as +a linear map from the domain to the codomain indices. These can be specified either via +their labels, or directly through a `biperm`. The output is a vector of eigenvalues. + +## Keyword arguments + +- `ishermitian::Bool`: specify if the matrix is Hermitian, which can be used to speed up the + computation. If `false`, the output `eltype` will always be `<:Complex`. +- Other keywords are passed on directly to MatrixAlgebraKit +""" +function eigvals(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) + biperm = blockedperm_indexin(Tuple.((labels_A, labels_codomain, labels_domain))...) + return eigvals(A, biperm; kwargs...) +end +function eigvals(A::AbstractArray, biperm::BlockedPermutation{2}; kwargs...) + A_mat = fusedims(A, biperm) + ishermitian = @something ishermitian LinearAlgebra.ishermitian(A_mat) + return (ishermitian ? eigh_vals : eig_vals)(A_mat; kwargs...) +end + # TODO: separate out the algorithm selection step from the implementation """ svd(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) -> U, S, Vᴴ @@ -156,3 +180,20 @@ function svd( axes_Vᴴ = (axes(Vᴴ, 1), axes_domain...) return splitdims(U, axes_U), S, splitdims(Vᴴ, axes_Vᴴ) end + +""" + svdvals(A::AbstractArray, labels_A, labels_codomain, labels_domain) -> S + svdvals(A::AbstractArray, biperm::BlockedPermutation{2}) -> S + +Compute the singular values of a generic N-dimensional array, by interpreting it as +a linear map from the domain to the codomain indices. These can be specified either via +their labels, or directly through a `biperm`. The output is a vector of singular values. +""" +function svdvals(A::AbstractArray, labels_A, labels_codomain, labels_domain) + biperm = blockedperm_indexin(Tuple.((labels_A, labels_codomain, labels_domain))...) + return svdvals(A, biperm) +end +function svdvals(A::AbstractArray, biperm::BlockedPermutation{2}) + A_mat = fusedims(A, biperm) + return svd_vals(A_mat) +end From ee10d27841ee38d07f5cc04b34af47f8eed479c2 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 12 Mar 2025 15:18:37 -0400 Subject: [PATCH 08/26] implement skeleton --- .github/workflows/IntegrationTestRequest.yml | 14 ++++++++++++++ docs/make.jl | 6 ++++-- docs/src/reference.md | 5 +++++ test/runtests.jl | 6 ++++-- 4 files changed, 27 insertions(+), 4 deletions(-) create mode 100644 .github/workflows/IntegrationTestRequest.yml create mode 100644 docs/src/reference.md diff --git a/.github/workflows/IntegrationTestRequest.yml b/.github/workflows/IntegrationTestRequest.yml new file mode 100644 index 0000000..d42fcca --- /dev/null +++ b/.github/workflows/IntegrationTestRequest.yml @@ -0,0 +1,14 @@ +name: "Integration Test Request" + +on: + issue_comment: + types: [created] + +jobs: + integrationrequest: + if: | + github.event.issue.pull_request && + contains(fromJSON('["OWNER", "COLLABORATOR", "MEMBER"]'), github.event.comment.author_association) + uses: ITensor/ITensorActions/.github/workflows/IntegrationTestRequest.yml@main + with: + localregistry: https://github.com/ITensor/ITensorRegistry.git diff --git a/docs/make.jl b/docs/make.jl index 8d582b7..3c6aa14 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,7 +1,9 @@ using TensorAlgebra: TensorAlgebra using Documenter: Documenter, DocMeta, deploydocs, makedocs -DocMeta.setdocmeta!(TensorAlgebra, :DocTestSetup, :(using TensorAlgebra); recursive=true) +DocMeta.setdocmeta!( + TensorAlgebra, :DocTestSetup, :(using TensorAlgebra); recursive=true +) include("make_index.jl") @@ -14,7 +16,7 @@ makedocs(; edit_link="main", assets=String[], ), - pages=["Home" => "index.md"], + pages=["Home" => "index.md", "Reference" => "reference.md"], ) deploydocs(; diff --git a/docs/src/reference.md b/docs/src/reference.md new file mode 100644 index 0000000..bc6be34 --- /dev/null +++ b/docs/src/reference.md @@ -0,0 +1,5 @@ +# Reference + +```@autodocs +Modules = [TensorAlgebra] +``` diff --git a/test/runtests.jl b/test/runtests.jl index bd97441..1c52c3e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -24,9 +24,11 @@ isexamplefile(fn) = # tests in groups based on folder structure for testgroup in filter(isdir, readdir(@__DIR__)) if GROUP == "ALL" || GROUP == uppercase(testgroup) - for file in filter(istestfile, readdir(joinpath(@__DIR__, testgroup); join=true)) + groupdir = joinpath(@__DIR__, testgroup) + for file in filter(istestfile, readdir(groupdir)) + filename = joinpath(groupdir, file) @eval @safetestset $file begin - include($file) + include($filename) end end end From 721db2a7922ab8eb576f483e42c3e4541b259396 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 12 Mar 2025 15:31:22 -0400 Subject: [PATCH 09/26] Fix exports test --- test/test_exports.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/test_exports.jl b/test/test_exports.jl index f07654d..7187a0a 100644 --- a/test/test_exports.jl +++ b/test/test_exports.jl @@ -1,6 +1,8 @@ using TensorAlgebra: TensorAlgebra using Test: @test, @testset @testset "Test exports" begin - exports = [:TensorAlgebra, :contract, :contract!] + exports = [ + :TensorAlgebra, :contract, :contract!, :eig, :eigvals, :lq, :qr, :svd, :svdvals + ] @test issetequal(names(TensorAlgebra), exports) end From 5fa0807a35bf6697c527dd2e3222848c7f890f37 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 12 Mar 2025 15:31:34 -0400 Subject: [PATCH 10/26] Fix factorization imports --- src/factorizations.jl | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/factorizations.jl b/src/factorizations.jl index 58a65b9..c3c2caf 100644 --- a/src/factorizations.jl +++ b/src/factorizations.jl @@ -1,15 +1,18 @@ using MatrixAlgebraKit: eig_full, - eigh_full, eig_trunc, + eig_vals, + eigh_full, eigh_trunc, + eigh_vals, lq_full, lq_compact, qr_full, qr_compact, svd_full, svd_compact, - svd_trunc + svd_trunc, + svd_vals using LinearAlgebra: LinearAlgebra # TODO: consider in-place version @@ -130,7 +133,9 @@ function eigvals(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwa biperm = blockedperm_indexin(Tuple.((labels_A, labels_codomain, labels_domain))...) return eigvals(A, biperm; kwargs...) end -function eigvals(A::AbstractArray, biperm::BlockedPermutation{2}; kwargs...) +function eigvals( + A::AbstractArray, biperm::BlockedPermutation{2}; ishermitian=nothing, kwargs... +) A_mat = fusedims(A, biperm) ishermitian = @something ishermitian LinearAlgebra.ishermitian(A_mat) return (ishermitian ? eigh_vals : eig_vals)(A_mat; kwargs...) From c0bdf887593868753e60dab01a88a4ca8c55180e Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 12 Mar 2025 15:31:48 -0400 Subject: [PATCH 11/26] Add `svdvals` and `eigvals` tests --- test/test_factorizations.jl | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl index fe9d708..baf8b0d 100644 --- a/test/test_factorizations.jl +++ b/test/test_factorizations.jl @@ -1,8 +1,8 @@ using Test: @test, @testset, @inferred using TestExtras: @constinferred -using TensorAlgebra: contract, qr, svd, TensorAlgebra, eig +using TensorAlgebra: contract, qr, svd, svdvals, TensorAlgebra, eig, eigvals using MatrixAlgebraKit: truncrank -using LinearAlgebra: norm +using LinearAlgebra: norm, diag elts = (Float64, ComplexF64) @@ -53,6 +53,10 @@ end AV = contract((:a, :b, :D), A, labels_A, V, (labels_V′..., :D)) VD = contract((:a, :b, :D), V, (labels_V..., :D′), D, (:D′, :D)) @test AV ≈ VD + + Dvals = @constinferred eigvals(A, labels_A, labels_V, labels_V′; ishermitian=false) + @test Dvals ≈ diag(D) + @test eltype(Dvals) <: Complex end @testset "Hermitian eigenvalue decomposition ($T)" for T in elts @@ -72,6 +76,11 @@ end AV = contract((:a, :b, :D), A, labels_A, V, (labels_V′..., :D)) VD = contract((:a, :b, :D), V, (labels_V..., :D′), D, (:D′, :D)) @test AV ≈ VD + + # TODO: broken type stability + Dvals = eigvals(A, labels_A, labels_V, labels_V′; ishermitian=true) + @test Dvals ≈ diag(D) + @test eltype(Dvals) <: Real end # Singular Value Decomposition @@ -106,6 +115,9 @@ end @test A ≈ A′ k = min(size(S)...) @test size(U, 3) == k == size(Vᴴ, 1) + + Svals = @constinferred svdvals(A, labels_A, labels_U, labels_Vᴴ) + @test Svals ≈ diag(S) end @testset "Truncated SVD ($T)" for T in elts From b9e5049ad1b05ecd71423418250b46466e9c0a13 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 12 Mar 2025 15:36:35 -0400 Subject: [PATCH 12/26] Add `lq` tests --- test/test_factorizations.jl | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl index baf8b0d..68a91ce 100644 --- a/test/test_factorizations.jl +++ b/test/test_factorizations.jl @@ -36,6 +36,36 @@ end @test size(Q, 3) == min(size(A, 1) * size(A, 2), size(A, 3) * size(A, 4)) end +# LQ Decomposition +# ---------------- +@testset "Full LQ ($T)" for T in elts + A = randn(T, 2, 3, 4, 5) + labels_A = (:a, :b, :c, :d) + labels_Q = (:d, :c) + labels_L = (:b, :a) + + Acopy = deepcopy(A) + L, Q = @constinferred lq(A, labels_A, labels_L, labels_Q; full=true) + @test A == Acopy # should not have altered initial array + A′ = contract(labels_A, L, (labels_L..., :q), Q, (:q, labels_Q...)) + @test A ≈ A′ + @test size(Q, 1) == size(Q, 2) * size(Q, 3) # Q is unitary +end + +@testset "Compact LQ ($T)" for T in elts + A = randn(T, 2, 3, 4, 5) + labels_A = (:a, :b, :c, :d) + labels_Q = (:d, :c) + labels_L = (:b, :a) + + Acopy = deepcopy(A) + L, Q = @constinferred lq(A, labels_A, labels_L, labels_Q; full=true) + @test A == Acopy # should not have altered initial array + A′ = contract(labels_A, L, (labels_L..., :q), Q, (:q, labels_Q...)) + @test A ≈ A′ + @test size(Q, 1) == min(size(A, 1) * size(A, 2), size(A, 3) * size(A, 4)) # Q is unitary +end + # Eigenvalue Decomposition # ------------------------ @testset "Eigenvalue decomposition ($T)" for T in elts From 08906c6befc6b09a37fd1abf746d4f33df93e0af Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 12 Mar 2025 15:36:58 -0400 Subject: [PATCH 13/26] Formatter --- docs/make.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 3c6aa14..e89c24e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,9 +1,7 @@ using TensorAlgebra: TensorAlgebra using Documenter: Documenter, DocMeta, deploydocs, makedocs -DocMeta.setdocmeta!( - TensorAlgebra, :DocTestSetup, :(using TensorAlgebra); recursive=true -) +DocMeta.setdocmeta!(TensorAlgebra, :DocTestSetup, :(using TensorAlgebra); recursive=true) include("make_index.jl") From cbc3ec61a065e985c14157e25c15a9d3ef8c6be1 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 12 Mar 2025 15:47:45 -0400 Subject: [PATCH 14/26] Fix some tests but also segfault ? --- test/test_factorizations.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl index 68a91ce..c8159f0 100644 --- a/test/test_factorizations.jl +++ b/test/test_factorizations.jl @@ -1,6 +1,6 @@ using Test: @test, @testset, @inferred using TestExtras: @constinferred -using TensorAlgebra: contract, qr, svd, svdvals, TensorAlgebra, eig, eigvals +using TensorAlgebra: TensorAlgebra, contract, lq, qr, svd, svdvals, eig, eigvals using MatrixAlgebraKit: truncrank using LinearAlgebra: norm, diag @@ -53,13 +53,13 @@ end end @testset "Compact LQ ($T)" for T in elts - A = randn(T, 2, 3, 4, 5) + A = randn(T, 5, 4, 3, 2) # compact only makes a difference for less rows labels_A = (:a, :b, :c, :d) labels_Q = (:d, :c) labels_L = (:b, :a) Acopy = deepcopy(A) - L, Q = @constinferred lq(A, labels_A, labels_L, labels_Q; full=true) + L, Q = @constinferred lq(A, labels_A, labels_L, labels_Q; full=false) @test A == Acopy # should not have altered initial array A′ = contract(labels_A, L, (labels_L..., :q), Q, (:q, labels_Q...)) @test A ≈ A′ @@ -84,7 +84,8 @@ end VD = contract((:a, :b, :D), V, (labels_V..., :D′), D, (:D′, :D)) @test AV ≈ VD - Dvals = @constinferred eigvals(A, labels_A, labels_V, labels_V′; ishermitian=false) + # TODO: broken type stability + Dvals = eigvals(A, labels_A, labels_V, labels_V′; ishermitian=false) @test Dvals ≈ diag(D) @test eltype(Dvals) <: Complex end From 7694766ffe289bec50d75b7025bba7bdb55c4c87 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 12 Mar 2025 16:29:29 -0400 Subject: [PATCH 15/26] Fix typo --- src/factorizations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/factorizations.jl b/src/factorizations.jl index c3c2caf..ec5b705 100644 --- a/src/factorizations.jl +++ b/src/factorizations.jl @@ -65,7 +65,7 @@ function lq(A::AbstractArray, biperm::BlockedPermutation{2}; full::Bool=true, kw # matrix to tensor axes_codomain, axes_domain = blockpermute(axes(A), biperm) - axes_L = (axes_codomain..., axes(L, 2)) + axes_L = (axes_codomain..., axes(L, ndims(L))) axes_Q = (axes(Q, 1), axes_domain...) return splitdims(L, axes_L), splitdims(Q, axes_Q) end From 959c96987d192969a9217d86d70db72a46c84b72 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 12 Mar 2025 16:29:40 -0400 Subject: [PATCH 16/26] Temporarily disable `lq` tests --- test/test_factorizations.jl | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl index c8159f0..9168860 100644 --- a/test/test_factorizations.jl +++ b/test/test_factorizations.jl @@ -52,19 +52,20 @@ end @test size(Q, 1) == size(Q, 2) * size(Q, 3) # Q is unitary end -@testset "Compact LQ ($T)" for T in elts - A = randn(T, 5, 4, 3, 2) # compact only makes a difference for less rows - labels_A = (:a, :b, :c, :d) - labels_Q = (:d, :c) - labels_L = (:b, :a) - - Acopy = deepcopy(A) - L, Q = @constinferred lq(A, labels_A, labels_L, labels_Q; full=false) - @test A == Acopy # should not have altered initial array - A′ = contract(labels_A, L, (labels_L..., :q), Q, (:q, labels_Q...)) - @test A ≈ A′ - @test size(Q, 1) == min(size(A, 1) * size(A, 2), size(A, 3) * size(A, 4)) # Q is unitary -end +# TODO: broken and segfaults because of MatrixAlgebraKit issue +# @testset "Compact LQ ($T)" for T in elts +# A = randn(T, 5, 4, 3, 2) # compact only makes a difference for less rows +# labels_A = (:a, :b, :c, :d) +# labels_Q = (:d, :c) +# labels_L = (:b, :a) + +# Acopy = deepcopy(A) +# L, Q = @constinferred lq(A, labels_A, labels_L, labels_Q; full=false) +# @test A == Acopy # should not have altered initial array +# A′ = contract(labels_A, L, (labels_L..., :q), Q, (:q, labels_Q...)) +# @test A ≈ A′ +# @test size(Q, 1) == min(size(A, 1) * size(A, 2), size(A, 3) * size(A, 4)) # Q is unitary +# end # Eigenvalue Decomposition # ------------------------ From bdc3e5818aed3f4267af82868e873cc598455204 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 13 Mar 2025 11:32:38 -0400 Subject: [PATCH 17/26] re-enable LQ tests --- Project.toml | 2 +- test/test_factorizations.jl | 27 +++++++++++++-------------- 2 files changed, 14 insertions(+), 15 deletions(-) diff --git a/Project.toml b/Project.toml index 1c2ade8..00bc3e1 100644 --- a/Project.toml +++ b/Project.toml @@ -24,7 +24,7 @@ BlockArrays = "1.2.0" EllipsisNotation = "1.8.0" GradedUnitRanges = "0.1.0" LinearAlgebra = "1.10" -MatrixAlgebraKit = "0.1.0" +MatrixAlgebraKit = "0.1.1" TupleTools = "1.6.0" TypeParameterAccessors = "0.2.1, 0.3" julia = "1.10" diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl index 9168860..c8159f0 100644 --- a/test/test_factorizations.jl +++ b/test/test_factorizations.jl @@ -52,20 +52,19 @@ end @test size(Q, 1) == size(Q, 2) * size(Q, 3) # Q is unitary end -# TODO: broken and segfaults because of MatrixAlgebraKit issue -# @testset "Compact LQ ($T)" for T in elts -# A = randn(T, 5, 4, 3, 2) # compact only makes a difference for less rows -# labels_A = (:a, :b, :c, :d) -# labels_Q = (:d, :c) -# labels_L = (:b, :a) - -# Acopy = deepcopy(A) -# L, Q = @constinferred lq(A, labels_A, labels_L, labels_Q; full=false) -# @test A == Acopy # should not have altered initial array -# A′ = contract(labels_A, L, (labels_L..., :q), Q, (:q, labels_Q...)) -# @test A ≈ A′ -# @test size(Q, 1) == min(size(A, 1) * size(A, 2), size(A, 3) * size(A, 4)) # Q is unitary -# end +@testset "Compact LQ ($T)" for T in elts + A = randn(T, 5, 4, 3, 2) # compact only makes a difference for less rows + labels_A = (:a, :b, :c, :d) + labels_Q = (:d, :c) + labels_L = (:b, :a) + + Acopy = deepcopy(A) + L, Q = @constinferred lq(A, labels_A, labels_L, labels_Q; full=false) + @test A == Acopy # should not have altered initial array + A′ = contract(labels_A, L, (labels_L..., :q), Q, (:q, labels_Q...)) + @test A ≈ A′ + @test size(Q, 1) == min(size(A, 1) * size(A, 2), size(A, 3) * size(A, 4)) # Q is unitary +end # Eigenvalue Decomposition # ------------------------ From 512cfadcef5bf9bd55eceab2cf933a35a3147c27 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 13 Mar 2025 15:04:25 -0400 Subject: [PATCH 18/26] default to `full=false` for `qr` and `lq` --- src/factorizations.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/factorizations.jl b/src/factorizations.jl index ec5b705..8b751f6 100644 --- a/src/factorizations.jl +++ b/src/factorizations.jl @@ -19,8 +19,8 @@ using LinearAlgebra: LinearAlgebra # TODO: figure out kwargs and document # """ - qr(A::AbstractArray, labels_A, labels_codomain, labels_domain; full=true, kwargs...) -> Q, R - qr(A::AbstractArray, biperm::BlockedPermutation{2}; full=true, kwargs...) -> Q, R + qr(A::AbstractArray, labels_A, labels_codomain, labels_domain; full=false, kwargs...) -> Q, R + qr(A::AbstractArray, biperm::BlockedPermutation{2}; full=false, kwargs...) -> Q, R Compute the QR decomposition of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain indices. These can be specified either via @@ -30,7 +30,7 @@ function qr(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs.. biperm = blockedperm_indexin(Tuple.((labels_A, labels_codomain, labels_domain))...) return qr(A, biperm) end -function qr(A::AbstractArray, biperm::BlockedPermutation{2}; full::Bool=true, kwargs...) +function qr(A::AbstractArray, biperm::BlockedPermutation{2}; full::Bool=false, kwargs...) # tensor to matrix A_mat = fusedims(A, biperm) @@ -45,8 +45,8 @@ function qr(A::AbstractArray, biperm::BlockedPermutation{2}; full::Bool=true, kw end """ - lq(A::AbstractArray, labels_A, labels_codomain, labels_domain; full=true, kwargs...) -> L, Q - lq(A::AbstractArray, biperm::BlockedPermutation{2}; full=true, kwargs...) -> L, Q + lq(A::AbstractArray, labels_A, labels_codomain, labels_domain; full=false, kwargs...) -> L, Q + lq(A::AbstractArray, biperm::BlockedPermutation{2}; full=false, kwargs...) -> L, Q Compute the LQ decomposition of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain indices. These can be specified either via @@ -56,7 +56,7 @@ function lq(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs.. biperm = blockedperm_indexin(Tuple.((labels_A, labels_codomain, labels_domain))...) return lq(A, biperm) end -function lq(A::AbstractArray, biperm::BlockedPermutation{2}; full::Bool=true, kwargs...) +function lq(A::AbstractArray, biperm::BlockedPermutation{2}; full::Bool=false, kwargs...) # tensor to matrix A_mat = fusedims(A, biperm) From b1e3104276c52c6da43e1a145d3f0f95ab298920 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 13 Mar 2025 15:05:58 -0400 Subject: [PATCH 19/26] Replace `eig` with `eigen` --- src/TensorAlgebra.jl | 2 +- src/factorizations.jl | 10 +++++----- test/test_exports.jl | 2 +- test/test_factorizations.jl | 6 +++--- 4 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/TensorAlgebra.jl b/src/TensorAlgebra.jl index 5f48988..fc3c157 100644 --- a/src/TensorAlgebra.jl +++ b/src/TensorAlgebra.jl @@ -1,6 +1,6 @@ module TensorAlgebra -export contract, contract!, eig, eigvals, lq, qr, svd, svdvals +export contract, contract!, eigen, eigvals, lq, qr, svd, svdvals include("blockedtuple.jl") include("blockedpermutation.jl") diff --git a/src/factorizations.jl b/src/factorizations.jl index 8b751f6..8b1a1ba 100644 --- a/src/factorizations.jl +++ b/src/factorizations.jl @@ -72,8 +72,8 @@ end # TODO: what name do we want? """ - eig(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) -> D, V - eig(A::AbstractArray, biperm::BlockedPermutation{2}; kwargs...) -> D, V + eigen(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) -> D, V + eigen(A::AbstractArray, biperm::BlockedPermutation{2}; kwargs...) -> D, V Compute the eigenvalue decomposition of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain indices. These can be specified either via @@ -86,11 +86,11 @@ their labels, or directly through a `biperm`. - `trunc`: Truncation keywords for `eig(h)_trunc`. - Other keywords are passed on directly to MatrixAlgebraKit """ -function eig(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) +function eigen(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) biperm = blockedperm_indexin(Tuple.((labels_A, labels_codomain, labels_domain))...) - return eig(A, biperm; kwargs...) + return eigen(A, biperm; kwargs...) end -function eig( +function eigen( A::AbstractArray, biperm::BlockedPermutation{2}; trunc=nothing, diff --git a/test/test_exports.jl b/test/test_exports.jl index 7187a0a..6a59cb5 100644 --- a/test/test_exports.jl +++ b/test/test_exports.jl @@ -2,7 +2,7 @@ using TensorAlgebra: TensorAlgebra using Test: @test, @testset @testset "Test exports" begin exports = [ - :TensorAlgebra, :contract, :contract!, :eig, :eigvals, :lq, :qr, :svd, :svdvals + :TensorAlgebra, :contract, :contract!, :eigen, :eigvals, :lq, :qr, :svd, :svdvals ] @test issetequal(names(TensorAlgebra), exports) end diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl index c8159f0..4f5b677 100644 --- a/test/test_factorizations.jl +++ b/test/test_factorizations.jl @@ -1,6 +1,6 @@ using Test: @test, @testset, @inferred using TestExtras: @constinferred -using TensorAlgebra: TensorAlgebra, contract, lq, qr, svd, svdvals, eig, eigvals +using TensorAlgebra: TensorAlgebra, contract, lq, qr, svd, svdvals, eigen, eigvals using MatrixAlgebraKit: truncrank using LinearAlgebra: norm, diag @@ -76,7 +76,7 @@ end Acopy = deepcopy(A) # type-unstable if `ishermitian` not set - D, V = @constinferred eig(A, labels_A, labels_V, labels_V′; ishermitian=false) + D, V = @constinferred eigen(A, labels_A, labels_V, labels_V′; ishermitian=false) @test A == Acopy # should not have altered initial array @test eltype(D) == eltype(V) && eltype(D) <: Complex @@ -99,7 +99,7 @@ end Acopy = deepcopy(A) # type-unstable if `ishermitian` not set - D, V = @constinferred eig(A, labels_A, labels_V, labels_V′; ishermitian=true) + D, V = @constinferred eigen(A, labels_A, labels_V, labels_V′; ishermitian=true) @test A == Acopy # should not have altered initial array @test eltype(D) <: Real @test eltype(V) == eltype(A) From 232a2f0198c9a727a49dd60b2a0b5645a8b21ae8 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 13 Mar 2025 15:10:48 -0400 Subject: [PATCH 20/26] Change for in-place factorizations --- src/factorizations.jl | 44 +++++++++++++++++++++---------------------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/src/factorizations.jl b/src/factorizations.jl index 8b1a1ba..338e4b3 100644 --- a/src/factorizations.jl +++ b/src/factorizations.jl @@ -1,18 +1,18 @@ using MatrixAlgebraKit: - eig_full, - eig_trunc, - eig_vals, - eigh_full, - eigh_trunc, - eigh_vals, - lq_full, - lq_compact, - qr_full, - qr_compact, - svd_full, - svd_compact, - svd_trunc, - svd_vals + eig_full!, + eig_trunc!, + eig_vals!, + eigh_full!, + eigh_trunc!, + eigh_vals!, + lq_full!, + lq_compact!, + qr_full!, + qr_compact!, + svd_full!, + svd_compact!, + svd_trunc!, + svd_vals! using LinearAlgebra: LinearAlgebra # TODO: consider in-place version @@ -35,7 +35,7 @@ function qr(A::AbstractArray, biperm::BlockedPermutation{2}; full::Bool=false, k A_mat = fusedims(A, biperm) # factorization - Q, R = full ? qr_full(A_mat; kwargs...) : qr_compact(A_mat; kwargs...) + Q, R = full ? qr_full!(A_mat; kwargs...) : qr_compact!(A_mat; kwargs...) # matrix to tensor axes_codomain, axes_domain = blockpermute(axes(A), biperm) @@ -61,7 +61,7 @@ function lq(A::AbstractArray, biperm::BlockedPermutation{2}; full::Bool=false, k A_mat = fusedims(A, biperm) # factorization - L, Q = full ? lq_full(A_mat; kwargs...) : lq_compact(A_mat; kwargs...) + L, Q = full ? lq_full!(A_mat; kwargs...) : lq_compact!(A_mat; kwargs...) # matrix to tensor axes_codomain, axes_domain = blockpermute(axes(A), biperm) @@ -104,9 +104,9 @@ function eigen( # factorization if !isnothing(trunc) - D, V = (ishermitian ? eigh_trunc : eig_trunc)(A_mat; trunc, kwargs...) + D, V = (ishermitian ? eigh_trunc! : eig_trunc!)(A_mat; trunc, kwargs...) else - D, V = (ishermitian ? eigh_full : eig_full)(A_mat; kwargs...) + D, V = (ishermitian ? eigh_full! : eig_full!)(A_mat; kwargs...) end # matrix to tensor @@ -138,7 +138,7 @@ function eigvals( ) A_mat = fusedims(A, biperm) ishermitian = @something ishermitian LinearAlgebra.ishermitian(A_mat) - return (ishermitian ? eigh_vals : eig_vals)(A_mat; kwargs...) + return (ishermitian ? eigh_vals! : eig_vals!)(A_mat; kwargs...) end # TODO: separate out the algorithm selection step from the implementation @@ -174,9 +174,9 @@ function svd( # factorization if !isnothing(trunc) @assert !full "Specified both full and truncation, currently not supported" - U, S, Vᴴ = svd_trunc(A_mat; trunc, kwargs...) + U, S, Vᴴ = svd_trunc!(A_mat; trunc, kwargs...) else - U, S, Vᴴ = full ? svd_full(A_mat; kwargs...) : svd_compact(A_mat; kwargs...) + U, S, Vᴴ = full ? svd_full!(A_mat; kwargs...) : svd_compact!(A_mat; kwargs...) end # matrix to tensor @@ -200,5 +200,5 @@ function svdvals(A::AbstractArray, labels_A, labels_codomain, labels_domain) end function svdvals(A::AbstractArray, biperm::BlockedPermutation{2}) A_mat = fusedims(A, biperm) - return svd_vals(A_mat) + return svd_vals!(A_mat) end From 3714df3dfe805e1de57578e31d21e64a444907fd Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 13 Mar 2025 16:15:13 -0400 Subject: [PATCH 21/26] Bump project 0.2.3 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 00bc3e1..f10c911 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "TensorAlgebra" uuid = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a" authors = ["ITensor developers and contributors"] -version = "0.2.2" +version = "0.2.3" [deps] ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" From 7c94b258c8636feaef699ca179463e5dece6c807 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 13 Mar 2025 18:04:50 -0400 Subject: [PATCH 22/26] Fix some tests --- src/factorizations.jl | 4 ++-- test/test_factorizations.jl | 12 ++++++------ 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/factorizations.jl b/src/factorizations.jl index 338e4b3..f9da496 100644 --- a/src/factorizations.jl +++ b/src/factorizations.jl @@ -28,7 +28,7 @@ their labels, or directly through a `biperm`. """ function qr(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) biperm = blockedperm_indexin(Tuple.((labels_A, labels_codomain, labels_domain))...) - return qr(A, biperm) + return qr(A, biperm; kwargs...) end function qr(A::AbstractArray, biperm::BlockedPermutation{2}; full::Bool=false, kwargs...) # tensor to matrix @@ -54,7 +54,7 @@ their labels, or directly through a `biperm`. """ function lq(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) biperm = blockedperm_indexin(Tuple.((labels_A, labels_codomain, labels_domain))...) - return lq(A, biperm) + return lq(A, biperm; kwargs...) end function lq(A::AbstractArray, biperm::BlockedPermutation{2}; full::Bool=false, kwargs...) # tensor to matrix diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl index 4f5b677..bef1d97 100644 --- a/test/test_factorizations.jl +++ b/test/test_factorizations.jl @@ -75,8 +75,8 @@ end labels_V′ = (:d, :c) Acopy = deepcopy(A) - # type-unstable if `ishermitian` not set - D, V = @constinferred eigen(A, labels_A, labels_V, labels_V′; ishermitian=false) + # type-unstable because of `ishermitian` difference + D, V = eigen(A, labels_A, labels_V, labels_V′; ishermitian=false) @test A == Acopy # should not have altered initial array @test eltype(D) == eltype(V) && eltype(D) <: Complex @@ -84,7 +84,7 @@ end VD = contract((:a, :b, :D), V, (labels_V..., :D′), D, (:D′, :D)) @test AV ≈ VD - # TODO: broken type stability + # type-unstable because of `ishermitian` difference Dvals = eigvals(A, labels_A, labels_V, labels_V′; ishermitian=false) @test Dvals ≈ diag(D) @test eltype(Dvals) <: Complex @@ -98,8 +98,8 @@ end labels_V′ = (:d, :c) Acopy = deepcopy(A) - # type-unstable if `ishermitian` not set - D, V = @constinferred eigen(A, labels_A, labels_V, labels_V′; ishermitian=true) + # type-unstable because of `ishermitian` difference + D, V = eigen(A, labels_A, labels_V, labels_V′; ishermitian=true) @test A == Acopy # should not have altered initial array @test eltype(D) <: Real @test eltype(V) == eltype(A) @@ -108,7 +108,7 @@ end VD = contract((:a, :b, :D), V, (labels_V..., :D′), D, (:D′, :D)) @test AV ≈ VD - # TODO: broken type stability + # type-unstable because of `ishermitian` difference Dvals = eigvals(A, labels_A, labels_V, labels_V′; ishermitian=true) @test Dvals ≈ diag(D) @test eltype(Dvals) <: Real From 0168bd6db10beb5744f9d8198195617cd2c16389 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 13 Mar 2025 18:09:56 -0400 Subject: [PATCH 23/26] Update docstrings --- src/factorizations.jl | 42 ++++++++++++++++++++++++++++++++---------- 1 file changed, 32 insertions(+), 10 deletions(-) diff --git a/src/factorizations.jl b/src/factorizations.jl index f9da496..831def3 100644 --- a/src/factorizations.jl +++ b/src/factorizations.jl @@ -15,16 +15,21 @@ using MatrixAlgebraKit: svd_vals! using LinearAlgebra: LinearAlgebra -# TODO: consider in-place version -# TODO: figure out kwargs and document -# """ - qr(A::AbstractArray, labels_A, labels_codomain, labels_domain; full=false, kwargs...) -> Q, R - qr(A::AbstractArray, biperm::BlockedPermutation{2}; full=false, kwargs...) -> Q, R + qr(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) -> Q, R + qr(A::AbstractArray, biperm::BlockedPermutation{2}; kwargs...) -> Q, R Compute the QR decomposition of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain indices. These can be specified either via their labels, or directly through a `biperm`. + +## Keyword arguments + +- `full::Bool=false`: select between a "full" or a "compact" decomposition, where `Q` is unitary or `R` is square, respectively. +- `positive::Bool=false`: specify if the diagonal of `R` should be positive, leading to a unique decomposition. +- Other keywords are passed on directly to MatrixAlgebraKit. + +See also `MatrixAlgebraKit.qr_full!` and `MatrixAlgebraKit.qr_compact!`. """ function qr(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) biperm = blockedperm_indexin(Tuple.((labels_A, labels_codomain, labels_domain))...) @@ -45,12 +50,20 @@ function qr(A::AbstractArray, biperm::BlockedPermutation{2}; full::Bool=false, k end """ - lq(A::AbstractArray, labels_A, labels_codomain, labels_domain; full=false, kwargs...) -> L, Q - lq(A::AbstractArray, biperm::BlockedPermutation{2}; full=false, kwargs...) -> L, Q + lq(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) -> L, Q + lq(A::AbstractArray, biperm::BlockedPermutation{2}; kwargs...) -> L, Q Compute the LQ decomposition of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain indices. These can be specified either via their labels, or directly through a `biperm`. + +## Keyword arguments + +- `full::Bool=false`: select between a "full" or a "compact" decomposition, where `Q` is unitary or `L` is square, respectively. +- `positive::Bool=false`: specify if the diagonal of `L` should be positive, leading to a unique decomposition. +- Other keywords are passed on directly to MatrixAlgebraKit. + +See also `MatrixAlgebraKit.lq_full!` and `MatrixAlgebraKit.lq_compact!`. """ function lq(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) biperm = blockedperm_indexin(Tuple.((labels_A, labels_codomain, labels_domain))...) @@ -84,7 +97,10 @@ their labels, or directly through a `biperm`. - `ishermitian::Bool`: specify if the matrix is Hermitian, which can be used to speed up the computation. If `false`, the output `eltype` will always be `<:Complex`. - `trunc`: Truncation keywords for `eig(h)_trunc`. -- Other keywords are passed on directly to MatrixAlgebraKit +- Other keywords are passed on directly to MatrixAlgebraKit. + +See also `MatrixAlgebraKit.eig_full!`, `MatrixAlgebraKit.eig_trunc!`, `MatrixAlgebraKit.eig_vals!`, +`MatrixAlgebraKit.eigh_full!`, `MatrixAlgebraKit.eigh_trunc!`, and `MatrixAlgebraKit.eigh_vals!`. """ function eigen(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) biperm = blockedperm_indexin(Tuple.((labels_A, labels_codomain, labels_domain))...) @@ -127,7 +143,9 @@ their labels, or directly through a `biperm`. The output is a vector of eigenval - `ishermitian::Bool`: specify if the matrix is Hermitian, which can be used to speed up the computation. If `false`, the output `eltype` will always be `<:Complex`. -- Other keywords are passed on directly to MatrixAlgebraKit +- Other keywords are passed on directly to MatrixAlgebraKit. + +See also `MatrixAlgebraKit.eig_vals!` and `MatrixAlgebraKit.eigh_vals!`. """ function eigvals(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) biperm = blockedperm_indexin(Tuple.((labels_A, labels_codomain, labels_domain))...) @@ -155,7 +173,9 @@ their labels, or directly through a `biperm`. - `full::Bool=false`: select between a "thick" or a "thin" decomposition, where both `U` and `Vᴴ` are unitary or isometric. - `trunc`: Truncation keywords for `svd_trunc`. Not compatible with `full=true`. -- Other keywords are passed on directly to MatrixAlgebraKit +- Other keywords are passed on directly to MatrixAlgebraKit. + +See also `MatrixAlgebraKit.svd_full!`, `MatrixAlgebraKit.svd_compact!`, and `MatrixAlgebraKit.svd_trunc!`. """ function svd(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) biperm = blockedperm_indexin(Tuple.((labels_A, labels_codomain, labels_domain))...) @@ -193,6 +213,8 @@ end Compute the singular values of a generic N-dimensional array, by interpreting it as a linear map from the domain to the codomain indices. These can be specified either via their labels, or directly through a `biperm`. The output is a vector of singular values. + +See also `MatrixAlgebraKit.svd_vals!`. """ function svdvals(A::AbstractArray, labels_A, labels_codomain, labels_domain) biperm = blockedperm_indexin(Tuple.((labels_A, labels_codomain, labels_domain))...) From ab353a90d5ee54dc9c272bd5b23c725d23823bd2 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 17 Mar 2025 09:39:40 -0400 Subject: [PATCH 24/26] Add `left_null` and `right_null` --- src/TensorAlgebra.jl | 2 +- src/factorizations.jl | 61 +++++++++++++++++++++++++++++++++++++ test/test_exports.jl | 11 ++++++- test/test_factorizations.jl | 28 +++++++++++++++-- 4 files changed, 98 insertions(+), 4 deletions(-) diff --git a/src/TensorAlgebra.jl b/src/TensorAlgebra.jl index fc3c157..591133b 100644 --- a/src/TensorAlgebra.jl +++ b/src/TensorAlgebra.jl @@ -1,6 +1,6 @@ module TensorAlgebra -export contract, contract!, eigen, eigvals, lq, qr, svd, svdvals +export contract, contract!, eigen, eigvals, lq, left_null, qr, right_null, svd, svdvals include("blockedtuple.jl") include("blockedpermutation.jl") diff --git a/src/factorizations.jl b/src/factorizations.jl index 831def3..0f61b2f 100644 --- a/src/factorizations.jl +++ b/src/factorizations.jl @@ -5,10 +5,12 @@ using MatrixAlgebraKit: eigh_full!, eigh_trunc!, eigh_vals!, + left_null!, lq_full!, lq_compact!, qr_full!, qr_compact!, + right_null!, svd_full!, svd_compact!, svd_trunc!, @@ -224,3 +226,62 @@ function svdvals(A::AbstractArray, biperm::BlockedPermutation{2}) A_mat = fusedims(A, biperm) return svd_vals!(A_mat) end + +""" + left_null(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) -> N + left_null(A::AbstractArray, biperm::BlockedPermutation{2}; kwargs...) -> N + +Compute the left nullspace of a generic N-dimensional array, by interpreting it as +a linear map from the domain to the codomain indices. These can be specified either via +their labels, or directly through a `biperm`. +The output satisfies `N' * A ≈ 0` and `N' * N ≈ I`. + +## Keyword arguments + +- `atol::Real=0`: absolute tolerance for the nullspace computation. +- `rtol::Real=0`: relative tolerance for the nullspace computation. +- `kind::Symbol`: specify the kind of decomposition used to compute the nullspace. + The options are `:qr`, `:qrpos` and `:svd`. The former two require `0 == atol == rtol`. + The default is `:qrpos` if `atol == rtol == 0`, and `:svd` otherwise. +""" +function left_null(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) + biperm = blockedperm_indexin(Tuple.((labels_A, labels_codomain, labels_domain))...) + return left_null(A, biperm; kwargs...) +end +function left_null(A::AbstractArray, biperm::BlockedPermutation{2}; kwargs...) + A_mat = fusedims(A, biperm) + N = left_null!(A_mat; kwargs...) + axes_codomain, _ = blockpermute(axes(A), biperm) + axes_N = (axes_codomain..., axes(N, 2)) + N_tensor = splitdims(N, axes_N) + return N_tensor +end + +""" + right_null(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) -> Nᴴ + right_null(A::AbstractArray, biperm::BlockedPermutation{2}; kwargs...) -> Nᴴ + +Compute the right nullspace of a generic N-dimensional array, by interpreting it as +a linear map from the domain to the codomain indices. These can be specified either via +their labels, or directly through a `biperm`. +The output satisfies `A * Nᴴ' ≈ 0` and `Nᴴ * Nᴴ' ≈ I`. + +## Keyword arguments + +- `atol::Real=0`: absolute tolerance for the nullspace computation. +- `rtol::Real=0`: relative tolerance for the nullspace computation. +- `kind::Symbol`: specify the kind of decomposition used to compute the nullspace. + The options are `:lq`, `:lqpos` and `:svd`. The former two require `0 == atol == rtol`. + The default is `:lqpos` if `atol == rtol == 0`, and `:svd` otherwise. +""" +function right_null(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) + biperm = blockedperm_indexin(Tuple.((labels_A, labels_codomain, labels_domain))...) + return right_null(A, biperm; kwargs...) +end +function right_null(A::AbstractArray, biperm::BlockedPermutation{2}; kwargs...) + A_mat = fusedims(A, biperm) + Nᴴ = right_null!(A_mat; kwargs...) + _, axes_domain = blockpermute(axes(A), biperm) + axes_Nᴴ = (axes(Nᴴ, 1), axes_domain...) + return splitdims(Nᴴ, axes_Nᴴ) +end diff --git a/test/test_exports.jl b/test/test_exports.jl index 6a59cb5..1e8693d 100644 --- a/test/test_exports.jl +++ b/test/test_exports.jl @@ -2,7 +2,16 @@ using TensorAlgebra: TensorAlgebra using Test: @test, @testset @testset "Test exports" begin exports = [ - :TensorAlgebra, :contract, :contract!, :eigen, :eigvals, :lq, :qr, :svd, :svdvals + :TensorAlgebra, + :contract, + :contract!, + :eigen, + :eigvals, + :left_null, + :lq, + :qr, + (:right_null):svd, + :svdvals, ] @test issetequal(names(TensorAlgebra), exports) end diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl index bef1d97..bdd944a 100644 --- a/test/test_factorizations.jl +++ b/test/test_factorizations.jl @@ -1,8 +1,9 @@ using Test: @test, @testset, @inferred using TestExtras: @constinferred -using TensorAlgebra: TensorAlgebra, contract, lq, qr, svd, svdvals, eigen, eigvals +using TensorAlgebra: + TensorAlgebra, contract, lq, qr, svd, svdvals, eigen, eigvals, left_null, right_null using MatrixAlgebraKit: truncrank -using LinearAlgebra: norm, diag +using LinearAlgebra: LinearAlgebra, norm, diag elts = (Float64, ComplexF64) @@ -170,3 +171,26 @@ end @test norm(A - A′) ≈ S_untrunc[end] @test size(S, 1) == size(S_untrunc, 1) - 1 end + +@testset "Nullspace ($T)" for T in elts + A = randn(T, 5, 4, 3, 2) + labels_A = (:a, :b, :c, :d) + labels_codomain = (:b, :a) + labels_domain = (:d, :c) + + Acopy = deepcopy(A) + N = @constinferred left_null(A, labels_A, labels_codomain, labels_domain) + @test A == Acopy # should not have altered initial array + # N^ba_n' * A^ba_dc = 0 + NA = contract((:n, labels_domain...), conj(N), (labels_codomain..., :n), A, labels_A) + @test norm(NA) ≈ 0 atol = 1e-14 + NN = contract((:n, :n′), conj(N), (labels_codomain..., :n), N, (labels_codomain..., :n′)) + @test NN ≈ LinearAlgebra.I + + Nᴴ = @constinferred right_null(A, labels_A, labels_codomain, labels_domain) + @test A == Acopy # should not have altered initial array + # A^ba_dc * N^dc_n' = 0 + AN = contract((labels_codomain..., :n), A, labels_A, conj(Nᴴ), (:n, labels_domain...)) + @test norm(AN) ≈ 0 atol = 1e-14 + NN = contract((:n, :n′), Nᴴ, (:n, labels_domain...), Nᴴ, (:n′, labels_domain...)) +end From 3355baa1e0656035f9986effca4c97d5794371d3 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 17 Mar 2025 09:51:36 -0400 Subject: [PATCH 25/26] Fix typo --- test/test_exports.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_exports.jl b/test/test_exports.jl index 1e8693d..dbef64f 100644 --- a/test/test_exports.jl +++ b/test/test_exports.jl @@ -10,7 +10,8 @@ using Test: @test, @testset :left_null, :lq, :qr, - (:right_null):svd, + :right_null, + :svd, :svdvals, ] @test issetequal(names(TensorAlgebra), exports) From eb6276f40a5d2c67dc1c8586e6b8be56ede60844 Mon Sep 17 00:00:00 2001 From: Matt Fishman Date: Mon, 17 Mar 2025 17:06:19 -0400 Subject: [PATCH 26/26] Remove outdated comment --- src/factorizations.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/factorizations.jl b/src/factorizations.jl index 0f61b2f..4adcadb 100644 --- a/src/factorizations.jl +++ b/src/factorizations.jl @@ -85,7 +85,6 @@ function lq(A::AbstractArray, biperm::BlockedPermutation{2}; full::Bool=false, k return splitdims(L, axes_L), splitdims(Q, axes_Q) end -# TODO: what name do we want? """ eigen(A::AbstractArray, labels_A, labels_codomain, labels_domain; kwargs...) -> D, V eigen(A::AbstractArray, biperm::BlockedPermutation{2}; kwargs...) -> D, V