From 1f4f04440521df91ba8179575e0a263b2f5c05d7 Mon Sep 17 00:00:00 2001 From: Matthew Priddin Date: Wed, 2 Nov 2022 11:52:48 +0000 Subject: [PATCH 1/7] Extend `lmul!` for `LowerTriangular` --- src/linalg.jl | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 1d2b0d5..52d56bc 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -45,7 +45,7 @@ end """ eigen_blockwise(B::BlockDiagonal, args...; kwargs...) -> values, vectors -Computes the eigendecomposition for each block separately and keeps the block diagonal +Computes the eigendecomposition for each block separately and keeps the block diagonal structure in the matrix of eigenvectors. Hence any parameters given are applied to each eigendecomposition separately, but there is e.g. no global sorting of eigenvalues. """ @@ -58,7 +58,7 @@ function eigen_blockwise(B::BlockDiagonal, args...; kwargs...) values = promote([e.values for e in eigens]...) vectors = promote([e.vectors for e in eigens]...) vcat(values...), BlockDiagonal([vectors...]) -end +end ## This function never keeps the block diagonal structure function LinearAlgebra.eigen(B::BlockDiagonal, args...; kwargs...) @@ -66,8 +66,8 @@ function LinearAlgebra.eigen(B::BlockDiagonal, args...; kwargs...) vectors = Matrix(vectors) # always convert to avoid type instability (also it speeds up the permutation step) @static if VERSION > v"1.2.0-DEV.275" Eigen(LinearAlgebra.sorteig!(values, vectors, kwargs...)...) - else - Eigen(values, vectors) + else + Eigen(values, vectors) end end @@ -157,6 +157,20 @@ function _mul!(C::BlockDiagonal, A::BlockDiagonal, B::BlockDiagonal, α::Number, return C end +function LinearAlgebra.lmul!(B::LowerTriangular{<:Any,<:BlockDiagonal}, vm::StridedVecOrMat) + row_i = 1 + # BlockDiagonals with non-square blocks + if !all(BlockDiagonals.is_square, blocks(parent(B))) + return lmul!(LowerTriangular(Matrix(B)), vm) # Fallback on the generic LinearAlgebra method + end + for block in blocks(parent(B)) + nrow = size(block, 1) + @views lmul!(LowerTriangular(block), vm[row_i:(row_i + nrow - 1), :]) + row_i += nrow + end + vm +end + function LinearAlgebra.:\(B::BlockDiagonal, vm::AbstractVecOrMat) row_i = 1 # BlockDiagonals with non-square blocks From 07a7e69e203a3ce1cb6fc682a8939d966a3b8219 Mon Sep 17 00:00:00 2001 From: Matthew Priddin Date: Wed, 2 Nov 2022 12:04:26 +0000 Subject: [PATCH 2/7] Add tests for `lmul!` for `LowerTriangular` blockdiagonals --- test/linalg.jl | 23 +++++++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) diff --git a/test/linalg.jl b/test/linalg.jl index d9f8977..2558e39 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -78,7 +78,7 @@ end evals, evecs = eigen(Matrix(B)) @test E isa Eigen - @test Matrix(E) ≈ B + @test Matrix(E) ≈ B # There is no test like @test eigen(B) == eigen(Matrix(B)) # 1. this fails in the complex case. Probably a convergence thing that leads to ever so slight differences @@ -88,7 +88,7 @@ end @static if VERSION < v"1.2" # pre-v1.2 we need to sort the eigenvalues consistently # Since eigenvalues may be complex here, I use this function, which works for this test. - # This test is already somewhat fragile w. r. t. degenerate eigenvalues + # This test is already somewhat fragile w. r. t. degenerate eigenvalues # and this just makes this a little worse. perm_bd = sortperm(real.(evals_bd) + 100*imag.(evals_bd)) evals_bd = evals_bd[perm_bd] @@ -131,7 +131,7 @@ end E = Eigen(block_vals, blocks(vecs)[i]) evals_bd, evecs_bd = E evals, evecs = eigen(block) - + @test block ≈ Matrix(E) @static if VERSION < v"1.2" @@ -144,7 +144,7 @@ end evals = evals[perm] evecs = evecs[:, perm] end - + @test evals_bd ≈ evals @test all(min.(abs.(evecs_bd - evecs), abs.(evecs_bd + evecs)) .< 1e-13) end @@ -245,6 +245,21 @@ end end end end # SVD + @testset "Left multiplication" begin + N1 = 20 + N2 = 8 + N3 = 5 + N4 = N1 + N3 - N2 + A = BlockDiagonal([rand(rng, N1, N1), rand(rng, N2, N2)]) + B = BlockDiagonal([rand(rng, N1, N2), rand(rng, N3, N4)]) + x = rand(rng, N1 + N2) + y = rand(rng, N2 + N4) + + @testset "Lower triangular" begin + @test lmul!(LowerTriangular(A), copy(x)) ≈ lmul!(LowerTriangular(Matrix(A)), x) + @test lmul!(LowerTriangular(B), copy(y)) ≈ lmul!(LowerTriangular(Matrix(B)), y) + end + end @testset "Left division" begin N1 = 20 N2 = 8 From dccbbd22838a2fa0fff58b68bf42d3964a899e04 Mon Sep 17 00:00:00 2001 From: Matthew Priddin Date: Wed, 2 Nov 2022 15:40:06 +0000 Subject: [PATCH 3/7] Add fallback to `+(::BlockDiagonal,::Diagonal)` for nonsquare blocks --- src/base_maths.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/base_maths.jl b/src/base_maths.jl index d74ccf1..97ba761 100644 --- a/src/base_maths.jl +++ b/src/base_maths.jl @@ -42,8 +42,11 @@ function Base.:+(B::BlockDiagonal, M::StridedMatrix) return A end -function Base.:+(B::BlockDiagonal, M::Diagonal)::BlockDiagonal +function Base.:+(B::BlockDiagonal, M::Diagonal) size(B) == size(M) || throw(DimensionMismatch("dimensions must match")) + if !all(is_square, blocks(B)) + return Matrix(B) + M # Fallback on the generic Base method + end A = copy(B) d = diag(M) row = 1 From b62ca2cac4146757e594a41175651a2a3fdc8b31 Mon Sep 17 00:00:00 2001 From: Matthew Priddin Date: Wed, 2 Nov 2022 15:40:27 +0000 Subject: [PATCH 4/7] Add tests --- test/base_maths.jl | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/test/base_maths.jl b/test/base_maths.jl index 4ee295d..675d334 100644 --- a/test/base_maths.jl +++ b/test/base_maths.jl @@ -19,6 +19,8 @@ using Test b64 = BlockDiagonal([rand(rng, 2, 2), rand(rng, 2, 2)]) b32 = BlockDiagonal([rand(rng, Float32, 2, 2), rand(rng, Float32, 2, 2)]) + bns = BlockDiagonal([rand(rng, N1, N2), rand(rng, N2, N3), rand(rng, N3, N1)]) + @testset "Addition" begin @testset "BlockDiagonal + BlockDiagonal" begin @test b1 + b1 isa BlockDiagonal @@ -58,6 +60,15 @@ using Test @test D + b1 isa BlockDiagonal @test D + b1 == D + Matrix(b1) @test_throws DimensionMismatch D′ + b1 + + # Non-square blocks + @test D + bns isa Matrix + @test D + bns == D + Matrix(bns) + @test_throws DimensionMismatch D′ + bns + + @test bns + D isa Matrix + @test bns + D == D + Matrix(bns) + @test_throws DimensionMismatch bns + D′ end @testset "BlockDiagonal + UniformScaling" begin @@ -69,11 +80,11 @@ using Test @test 5I + b1 == 5I + Matrix(b1) end end # Addition - + @testset "Subtraction" begin @test -b1 isa BlockDiagonal @test b1 - b1 isa BlockDiagonal - + @test -b1 == -Matrix(b1) @test b1 - b1 == Matrix(b1) - Matrix(b1) @test Matrix(b1) - b2 == Matrix(b1) - Matrix(b2) From d47ccdbf96aa5e4f2b44bdfaa0bea0bedd79e81a Mon Sep 17 00:00:00 2001 From: Matthew Priddin Date: Wed, 2 Nov 2022 15:45:21 +0000 Subject: [PATCH 5/7] Revert "Add tests for `lmul!` for `LowerTriangular` blockdiagonals" This reverts commit 07a7e69e203a3ce1cb6fc682a8939d966a3b8219. --- test/linalg.jl | 23 ++++------------------- 1 file changed, 4 insertions(+), 19 deletions(-) diff --git a/test/linalg.jl b/test/linalg.jl index 2558e39..d9f8977 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -78,7 +78,7 @@ end evals, evecs = eigen(Matrix(B)) @test E isa Eigen - @test Matrix(E) ≈ B + @test Matrix(E) ≈ B # There is no test like @test eigen(B) == eigen(Matrix(B)) # 1. this fails in the complex case. Probably a convergence thing that leads to ever so slight differences @@ -88,7 +88,7 @@ end @static if VERSION < v"1.2" # pre-v1.2 we need to sort the eigenvalues consistently # Since eigenvalues may be complex here, I use this function, which works for this test. - # This test is already somewhat fragile w. r. t. degenerate eigenvalues + # This test is already somewhat fragile w. r. t. degenerate eigenvalues # and this just makes this a little worse. perm_bd = sortperm(real.(evals_bd) + 100*imag.(evals_bd)) evals_bd = evals_bd[perm_bd] @@ -131,7 +131,7 @@ end E = Eigen(block_vals, blocks(vecs)[i]) evals_bd, evecs_bd = E evals, evecs = eigen(block) - + @test block ≈ Matrix(E) @static if VERSION < v"1.2" @@ -144,7 +144,7 @@ end evals = evals[perm] evecs = evecs[:, perm] end - + @test evals_bd ≈ evals @test all(min.(abs.(evecs_bd - evecs), abs.(evecs_bd + evecs)) .< 1e-13) end @@ -245,21 +245,6 @@ end end end end # SVD - @testset "Left multiplication" begin - N1 = 20 - N2 = 8 - N3 = 5 - N4 = N1 + N3 - N2 - A = BlockDiagonal([rand(rng, N1, N1), rand(rng, N2, N2)]) - B = BlockDiagonal([rand(rng, N1, N2), rand(rng, N3, N4)]) - x = rand(rng, N1 + N2) - y = rand(rng, N2 + N4) - - @testset "Lower triangular" begin - @test lmul!(LowerTriangular(A), copy(x)) ≈ lmul!(LowerTriangular(Matrix(A)), x) - @test lmul!(LowerTriangular(B), copy(y)) ≈ lmul!(LowerTriangular(Matrix(B)), y) - end - end @testset "Left division" begin N1 = 20 N2 = 8 From b62b74820a7f9b9f49cf1496d328496da93fad82 Mon Sep 17 00:00:00 2001 From: Matthew Priddin Date: Wed, 2 Nov 2022 15:46:15 +0000 Subject: [PATCH 6/7] Revert "Extend `lmul!` for `LowerTriangular`" This reverts commit 1f4f04440521df91ba8179575e0a263b2f5c05d7. --- src/linalg.jl | 22 ++++------------------ 1 file changed, 4 insertions(+), 18 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 52d56bc..1d2b0d5 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -45,7 +45,7 @@ end """ eigen_blockwise(B::BlockDiagonal, args...; kwargs...) -> values, vectors -Computes the eigendecomposition for each block separately and keeps the block diagonal +Computes the eigendecomposition for each block separately and keeps the block diagonal structure in the matrix of eigenvectors. Hence any parameters given are applied to each eigendecomposition separately, but there is e.g. no global sorting of eigenvalues. """ @@ -58,7 +58,7 @@ function eigen_blockwise(B::BlockDiagonal, args...; kwargs...) values = promote([e.values for e in eigens]...) vectors = promote([e.vectors for e in eigens]...) vcat(values...), BlockDiagonal([vectors...]) -end +end ## This function never keeps the block diagonal structure function LinearAlgebra.eigen(B::BlockDiagonal, args...; kwargs...) @@ -66,8 +66,8 @@ function LinearAlgebra.eigen(B::BlockDiagonal, args...; kwargs...) vectors = Matrix(vectors) # always convert to avoid type instability (also it speeds up the permutation step) @static if VERSION > v"1.2.0-DEV.275" Eigen(LinearAlgebra.sorteig!(values, vectors, kwargs...)...) - else - Eigen(values, vectors) + else + Eigen(values, vectors) end end @@ -157,20 +157,6 @@ function _mul!(C::BlockDiagonal, A::BlockDiagonal, B::BlockDiagonal, α::Number, return C end -function LinearAlgebra.lmul!(B::LowerTriangular{<:Any,<:BlockDiagonal}, vm::StridedVecOrMat) - row_i = 1 - # BlockDiagonals with non-square blocks - if !all(BlockDiagonals.is_square, blocks(parent(B))) - return lmul!(LowerTriangular(Matrix(B)), vm) # Fallback on the generic LinearAlgebra method - end - for block in blocks(parent(B)) - nrow = size(block, 1) - @views lmul!(LowerTriangular(block), vm[row_i:(row_i + nrow - 1), :]) - row_i += nrow - end - vm -end - function LinearAlgebra.:\(B::BlockDiagonal, vm::AbstractVecOrMat) row_i = 1 # BlockDiagonals with non-square blocks From b6b7048dceea748ad5c456eaf66723c18404d35f Mon Sep 17 00:00:00 2001 From: mjp98 Date: Wed, 2 Nov 2022 17:24:07 +0000 Subject: [PATCH 7/7] Bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index f4b67ad..29d8539 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "BlockDiagonals" uuid = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0" authors = ["Invenia Technical Computing Corporation"] -version = "0.1.37" +version = "0.1.38" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"