Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add left division operator #82

Merged
merged 10 commits into from
Nov 24, 2021
Merged

Add left division operator #82

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
09d8740
Add left division operator for block diag
theogf Nov 23, 2021
ffa4504
Moved everything to linalg instead
theogf Nov 23, 2021
0a70faf
Unncessary endline
theogf Nov 23, 2021
5d0574d
Change variable names
theogf Nov 24, 2021
41e16b2
Add check for square matrices
theogf Nov 24, 2021
139369e
Add suggestions
theogf Nov 24, 2021
b2b5704
Add test for error throws and reuse existing vectors
theogf Nov 24, 2021
de92bd9
Version bump
theogf Nov 24, 2021
3d0d461
Rename B variable in Cholesky and SVD tests
theogf Nov 24, 2021
7127e3f
Implemented suggestions
theogf Nov 24, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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"
Expand Down
15 changes: 15 additions & 0 deletions src/linalg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -149,3 +149,18 @@ function _mul!(C::BlockDiagonal, A::BlockDiagonal, B::BlockDiagonal, α::Number,

return C
end

function LinearAlgebra.:\(B::BlockDiagonal, vm::AbstractVecOrMat)
row_i = 1
# BlockDiagonals with non-square blocks
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)
result[row_i:(row_i + nrow - 1), :] = block \ vm[row_i:(row_i + nrow - 1), :]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i've not thought too much, but perhaps you know: is there any rewriting of this that would allow us to access vm col-wise rather than row-wise?

also does using views of vm help performance here at all?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think accessing vm colwise is possible for this problem.
views could help I suppose. I am going to run some benchmark and see

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Views does not seem to help that much in this simple example

julia> using BenchmarkTools
julia> N = 400
julia> A = BlockDiagonal([rand(N, N) for _ in 1:5])
julia> x = rand(5N, 100)
julia> A \ x
julia> @btime $A \ $x # No views
  7.085 ms (37 allocations: 10.70 MiB)
julia> @btime $A \ $x # With views
  7.403 ms (27 allocations: 9.17 MiB)

Copy link
Contributor

@nickrobinson251 nickrobinson251 Nov 24, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i think the reduced allocations suggests let's use views:

e.g. on a matrix of the same size, but more blocks:

julia> M = 40
julia> B = BlockDiagonal([rand(M, M) for _ in 1:5])
julia> x = rand(5N, 100)
julia> @btime $A \ $x # No views
  2.155 ms (302 allocations: 5.22 MiB)
julia> @btime $A \ $x # With views
  2.264 ms (202 allocations: 3.69 MiB)

(i know the times are ~equal but i usually trust allocations as at least as good a guide in practice)

row_i += nrow
end
return result
end
35 changes: 22 additions & 13 deletions test/linalg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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}}
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -229,4 +227,15 @@ using Test
end
end
end # SVD
@testset "Left division" begin
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