From 09d87405528a1dcacd27f63addeb8da2897104c0 Mon Sep 17 00:00:00 2001 From: Theo Galy-Fajou Date: Tue, 23 Nov 2021 21:33:28 +0100 Subject: [PATCH 01/10] Add left division operator for block diag --- src/base_maths.jl | 11 +++++++++++ test/base_maths.jl | 5 +++++ 2 files changed, 16 insertions(+) diff --git a/src/base_maths.jl b/src/base_maths.jl index c4eb24e..945e4d7 100644 --- a/src/base_maths.jl +++ b/src/base_maths.jl @@ -166,3 +166,14 @@ end ## Division Base.:/(B::BlockDiagonal, n::Number) = BlockDiagonal(blocks(B) ./ n) + +function LinearAlgebra.:\(A::BlockDiagonal, B::AbstractVecOrMat) + i = 1 + c = similar(B) + for a in blocks(A) + d = size(a, 1) + c[i:i+d-1, :] = a \ B[i:i+d-1, :] + i += d + end + return c +end diff --git a/test/base_maths.jl b/test/base_maths.jl index c1cbeab..7a33441 100644 --- a/test/base_maths.jl +++ b/test/base_maths.jl @@ -144,5 +144,10 @@ using Test @test sum(size.(b4.blocks, 1)) == size(b4 * b5, 1) @test sum(size.(b5.blocks, 2)) == size(b4 * b5, 2) end + + @testset "Left division" begin + x = rand(rng, N1 + N2 + N3) + @test b1 \ x ≈ inv(b1) * x + end end # Multiplication end From ffa4504f43c8b32b81c70a07dbd060993cff0d80 Mon Sep 17 00:00:00 2001 From: Theo Galy-Fajou Date: Tue, 23 Nov 2021 21:48:03 +0100 Subject: [PATCH 02/10] Moved everything to linalg instead --- src/base_maths.jl | 10 ---------- src/linalg.jl | 11 +++++++++++ test/base_maths.jl | 5 ----- test/linalg.jl | 4 ++++ 4 files changed, 15 insertions(+), 15 deletions(-) diff --git a/src/base_maths.jl b/src/base_maths.jl index 945e4d7..8bc1b20 100644 --- a/src/base_maths.jl +++ b/src/base_maths.jl @@ -167,13 +167,3 @@ end ## Division Base.:/(B::BlockDiagonal, n::Number) = BlockDiagonal(blocks(B) ./ n) -function LinearAlgebra.:\(A::BlockDiagonal, B::AbstractVecOrMat) - i = 1 - c = similar(B) - for a in blocks(A) - d = size(a, 1) - c[i:i+d-1, :] = a \ B[i:i+d-1, :] - i += d - end - return c -end diff --git a/src/linalg.jl b/src/linalg.jl index 2ea9b42..08ec318 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -149,3 +149,14 @@ function _mul!(C::BlockDiagonal, A::BlockDiagonal, B::BlockDiagonal, α::Number, return C end + +function LinearAlgebra.:\(A::BlockDiagonal, B::AbstractVecOrMat) + i = 1 + c = similar(B) + for a in blocks(A) + d = size(a, 1) + c[i:i+d-1, :] = a \ B[i:i+d-1, :] + i += d + end + return c +end diff --git a/test/base_maths.jl b/test/base_maths.jl index 7a33441..c1cbeab 100644 --- a/test/base_maths.jl +++ b/test/base_maths.jl @@ -144,10 +144,5 @@ using Test @test sum(size.(b4.blocks, 1)) == size(b4 * b5, 1) @test sum(size.(b5.blocks, 2)) == size(b4 * b5, 2) end - - @testset "Left division" begin - x = rand(rng, N1 + N2 + N3) - @test b1 \ x ≈ inv(b1) * x - end end # Multiplication end diff --git a/test/linalg.jl b/test/linalg.jl index 86634a1..ff1639c 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -229,4 +229,8 @@ using Test end end end # SVD + @testset "Left division" begin + x = rand(rng, N1 + N2 + N3) + @test b1 \ x ≈ inv(b1) * x + end end From 0a70faf4c285b2cbf9e058d4d0ff3c21f73d7e35 Mon Sep 17 00:00:00 2001 From: Theo Galy-Fajou Date: Tue, 23 Nov 2021 21:50:00 +0100 Subject: [PATCH 03/10] Unncessary endline --- src/base_maths.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/base_maths.jl b/src/base_maths.jl index 8bc1b20..c4eb24e 100644 --- a/src/base_maths.jl +++ b/src/base_maths.jl @@ -166,4 +166,3 @@ end ## Division Base.:/(B::BlockDiagonal, n::Number) = BlockDiagonal(blocks(B) ./ n) - From 5d0574d25dc3423a7984851138cc6772b3074be3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Wed, 24 Nov 2021 12:20:29 +0100 Subject: [PATCH 04/10] Change variable names Co-authored-by: Miha Zgubic --- src/linalg.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 08ec318..574122c 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -150,13 +150,13 @@ function _mul!(C::BlockDiagonal, A::BlockDiagonal, B::BlockDiagonal, α::Number, return C end -function LinearAlgebra.:\(A::BlockDiagonal, B::AbstractVecOrMat) - i = 1 - c = similar(B) - for a in blocks(A) - d = size(a, 1) - c[i:i+d-1, :] = a \ B[i:i+d-1, :] - i += d +function LinearAlgebra.:\(B::BlockDiagonal, vm::AbstractVecOrMat) + row_i = 1 + result = similar(vm) + for block in blocks(B) + nrow = size(block, 1) + result[row_i:row_i+nrow-1, :] = block \ vm[row_i:row_i+nrow-1, :] + row_i += nrow end - return c + return result end From 41e16b29ed7bd9976e192f7665bd68307cd8fc2b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Wed, 24 Nov 2021 12:29:01 +0100 Subject: [PATCH 05/10] Add check for square matrices --- src/linalg.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/linalg.jl b/src/linalg.jl index 574122c..7b24f66 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -150,8 +150,15 @@ function _mul!(C::BlockDiagonal, A::BlockDiagonal, B::BlockDiagonal, α::Number, return C end +# There is no restriction on x to be able to pass Cholesky objexts +function issquare(x) + indsm, indsn = size(x) + n1 == n2 +end + function LinearAlgebra.:\(B::BlockDiagonal, vm::AbstractVecOrMat) row_i = 1 + all(issquare, blocks(B)) || error("left division is not compatible with non-square matrices") result = similar(vm) for block in blocks(B) nrow = size(block, 1) From 139369ef9b35da210d02f9bd19610ce1268e6d7d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Wed, 24 Nov 2021 13:01:40 +0100 Subject: [PATCH 06/10] Add suggestions Co-authored-by: Miha Zgubic --- src/linalg.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 7b24f66..241d66a 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -150,19 +150,19 @@ function _mul!(C::BlockDiagonal, A::BlockDiagonal, B::BlockDiagonal, α::Number, return C end -# There is no restriction on x to be able to pass Cholesky objexts -function issquare(x) - indsm, indsn = size(x) - n1 == n2 -end - function LinearAlgebra.:\(B::BlockDiagonal, vm::AbstractVecOrMat) row_i = 1 - all(issquare, blocks(B)) || error("left division is not compatible with non-square matrices") + # BlockDiagonals with non-square blocks + all(b -> size(b, 1) == size(b, 2), blocks(B)) || throw( + ArgumentError( + "Left division (\\) is not defined for BlockDiagonals with " * + "non-square blocks (they are singular matrices).", + ), + ) result = similar(vm) for block in blocks(B) nrow = size(block, 1) - result[row_i:row_i+nrow-1, :] = block \ vm[row_i:row_i+nrow-1, :] + result[row_i:(row_i + nrow - 1), :] = block \ vm[row_i:(row_i + nrow - 1), :] row_i += nrow end return result From b2b5704ff3a923f1cee9c19beada1506f3311840 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Wed, 24 Nov 2021 13:05:29 +0100 Subject: [PATCH 07/10] Add test for error throws and reuse existing vectors --- test/linalg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/linalg.jl b/test/linalg.jl index ff1639c..508c39b 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -230,7 +230,7 @@ using Test end end # SVD @testset "Left division" begin - x = rand(rng, N1 + N2 + N3) - @test b1 \ x ≈ inv(b1) * x + @test b1 \ a ≈ inv(b1) * a + @test_throws ArgumentError BlockDiagonal([A, B]) \ rand(rng, 2N + N1 + N2) end end From de92bd9f509ec812533b293cc9f17fe797b2b85c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Wed, 24 Nov 2021 13:05:59 +0100 Subject: [PATCH 08/10] Version bump --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 1c5d79b..73c9639 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.23" +version = "0.1.24" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" From 3d0d461e20c074999601f6210912cc6a1240ee1e Mon Sep 17 00:00:00 2001 From: Theo Galy-Fajou Date: Wed, 24 Nov 2021 13:27:00 +0100 Subject: [PATCH 09/10] Rename B variable in Cholesky and SVD tests --- test/linalg.jl | 24 +++++++++++------------- 1 file changed, 11 insertions(+), 13 deletions(-) diff --git a/test/linalg.jl b/test/linalg.jl index 508c39b..aafd705 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -154,7 +154,6 @@ using Test @test logdet(b̂1) ≈ logdet(Matrix(b̂1)) end end # Unary - @testset "Cholesky decomposition" begin X = [ 4 12 -16 12 37 -43 @@ -163,10 +162,10 @@ using Test 0.0 1.0 5.0 0.0 0.0 3.0] - B = BlockDiagonal([X, X]) - C = cholesky(B) + BD = BlockDiagonal([X, X]) + C = cholesky(BD) @test C isa Cholesky{Float64, <:BlockDiagonal{Float64}} - @test C.U ≈ cholesky(Matrix(B)).U + @test C.U ≈ cholesky(Matrix(BD)).U @test C.U ≈ BlockDiagonal([U, U]) @test C.L ≈ BlockDiagonal([U', U']) @test C.UL ≈ C.U @@ -175,26 +174,25 @@ using Test M = BlockDiagonal(map(Matrix, blocks(C.L))) C = Cholesky(M, 'L', 0) - @test C.U ≈ cholesky(Matrix(B)).U + @test C.U ≈ cholesky(Matrix(BD)).U @test C.U ≈ BlockDiagonal([U, U]) @test C.L ≈ BlockDiagonal([U', U']) @test C.UL ≈ C.L @test C.uplo === 'L' @test C.info == 0 end # Cholesky - @testset "Singular Value Decomposition" begin X = [ 4 12 -16 12 37 -43 -16 -43 98] - B = BlockDiagonal([X, X]) + BD = BlockDiagonal([X, X]) @testset "full=$full" for full in (true, false) @testset "svd_blockwise" begin - U, S, Vt = svd_blockwise(B; full=full) + U, S, Vt = svd_blockwise(BD; full=full) F = SVD(U, S, Vt) - @test B ≈ F.U * Diagonal(F.S) * F.Vt + @test BD ≈ F.U * Diagonal(F.S) * F.Vt # Matrices should be BlockDiagonal @test F isa SVD{Float64, Float64, <:BlockDiagonal{Float64}} @@ -203,7 +201,7 @@ using Test @test F.Vt isa BlockDiagonal # Should have same values, but not sorted so as to keep BlockDiagonal structure - F_ = svd(Matrix(B), full=full) + F_ = svd(Matrix(BD), full=full) for fname in fieldnames(SVD) @test sort(vec(getfield(F, fname))) ≈ sort(vec(getfield(F_, fname))) end @@ -213,11 +211,11 @@ using Test end @testset "svd" begin - F = svd(B; full=full) - F_ = svd(Matrix(B), full=full) + F = svd(BD; full=full) + F_ = svd(Matrix(BD), full=full) @test F isa SVD - @test B ≈ F.U * Diagonal(F.S) * F.Vt + @test BD ≈ F.U * Diagonal(F.S) * F.Vt @test F == F_ for fname in fieldnames(SVD) From 7127e3f76ad70403eac7f98e79e2ae07cd147457 Mon Sep 17 00:00:00 2001 From: Theo Galy-Fajou Date: Wed, 24 Nov 2021 15:17:46 +0100 Subject: [PATCH 10/10] Implemented suggestions --- src/linalg.jl | 9 +++------ test/linalg.jl | 11 +++++++++-- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 241d66a..82d56b8 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -153,12 +153,9 @@ end function LinearAlgebra.:\(B::BlockDiagonal, vm::AbstractVecOrMat) row_i = 1 # BlockDiagonals with non-square blocks - all(b -> size(b, 1) == size(b, 2), blocks(B)) || throw( - ArgumentError( - "Left division (\\) is not defined for BlockDiagonals with " * - "non-square blocks (they are singular matrices).", - ), - ) + if !all(is_square, blocks(B)) + return Matrix(B) \ vm # Fallback on the generic LinearAlgebra method + end result = similar(vm) for block in blocks(B) nrow = size(block, 1) diff --git a/test/linalg.jl b/test/linalg.jl index aafd705..457fd0c 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -228,7 +228,14 @@ using Test end end # SVD @testset "Left division" begin - @test b1 \ a ≈ inv(b1) * a - @test_throws ArgumentError BlockDiagonal([A, B]) \ rand(rng, 2N + N1 + N2) + N1 = 20 + N2 = 8 + N3 = 5 + A = BlockDiagonal([rand(rng, N1, N1), rand(rng, N2, N2)]) + B = BlockDiagonal([rand(rng, N1, N2), rand(rng, N3, N1)]) + x = rand(rng, N1 + N2) + y = rand(rng, N1 + N3) + @test A \ x ≈ inv(A) * x + @test B \ y ≈ Matrix(B) \ y end end