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.

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

RFC: Allows blocks to be <:AbstractVector or <:Tuple #56

Closed
wants to merge 5 commits into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
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
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.14"
version = "0.2.0"

[deps]
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
Expand Down
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@

[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://invenia.github.io/BlockDiagonals.jl/stable)
[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://invenia.github.io/BlockDiagonals.jl/dev)
[![CI](https://github.com/invenia/BlockDiagonals.jl/workflows/CI/badge.svg)](https://github.com/Invenia/BlockDiagonals.jl/actions?query=workflow%3ACI)
[![CI](https://github.com/invenia/BlockDiagonals.jl/workflows/CI/badge.svg)](https://github.com/Invenia/BlockDiagonals.jl/actions?query=workflow:CI)
[![Codecov](https://codecov.io/gh/invenia/BlockDiagonals.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/invenia/BlockDiagonals.jl)
[![code style blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle)

Expand Down
63 changes: 47 additions & 16 deletions docs/Manifest.toml
Original file line number Diff line number Diff line change
Expand Up @@ -4,48 +4,71 @@
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[BlockDiagonals]]
deps = ["FillArrays", "LinearAlgebra"]
deps = ["ChainRulesCore", "FillArrays", "LinearAlgebra"]
path = ".."
uuid = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0"
version = "0.1.0"
version = "0.2.0"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "de4f08843c332d355852721adb1592bce7924da3"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "0.9.29"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "919c7f3151e79ff196add81d7f4e45d91bbf420b"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.25.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[DocStringExtensions]]
deps = ["LibGit2", "Markdown", "Pkg", "Test"]
git-tree-sha1 = "0513f1a8991e9d83255e0140aace0d0fc4486600"
git-tree-sha1 = "50ddf44c53698f5e784bbebb3f4b21c5807401b1"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.0"
version = "0.8.3"

[[Documenter]]
deps = ["Base64", "DocStringExtensions", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
git-tree-sha1 = "4a84478277020abfff208cde31ba1aa68a5bc572"
deps = ["Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
git-tree-sha1 = "21fb992ef1b28ff8f315354d3808ebf4a8fa6e45"
uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
version = "0.23.0"
version = "0.26.2"

[[FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Test"]
git-tree-sha1 = "9ab8f76758cbabba8d7f103c51dce7f73fcf8e92"
deps = ["LinearAlgebra", "Random", "SparseArrays"]
git-tree-sha1 = "ff537e5a3cba92fb48f30fec46723510450f2c0e"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.6.3"
version = "0.10.2"

[[IOCapture]]
deps = ["Logging"]
git-tree-sha1 = "377252859f740c217b936cebcd918a44f9b53b59"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.1.1"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "b34d7cef7b337321e97d22242c3c2b91f476748e"
git-tree-sha1 = "81690084b6198a2e1da36fcfda16eeca9f9f24e4"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.0"
version = "0.21.1"

[[LibGit2]]
deps = ["Printf"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[Libdl]]
Expand All @@ -66,13 +89,13 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Parsers]]
deps = ["Dates", "Test"]
git-tree-sha1 = "db2b35dedab3c0e46dc15996d170af07a5ab91c9"
deps = ["Dates"]
git-tree-sha1 = "50c9a9ed8c714945e01cd53a21007ed3865ed714"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "0.3.6"
version = "1.0.15"

[[Pkg]]
deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[Printf]]
Expand All @@ -93,13 +116,21 @@ uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[Test]]
deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
Expand Down
2 changes: 1 addition & 1 deletion docs/Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -3,4 +3,4 @@ BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0"
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"

[compat]
Documenter = "~0.23"
Documenter = "0.26"
10 changes: 7 additions & 3 deletions docs/make.jl
Original file line number Diff line number Diff line change
@@ -1,16 +1,20 @@
using Documenter, BlockDiagonals
using BlockDiagonals
using Documenter

makedocs(;
modules=[BlockDiagonals],
format=Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
format=Documenter.HTML(prettyurls=false),
pages=[
"Home" => "index.md",
],
repo="https://github.com/invenia/BlockDiagonals.jl/blob/{commit}{path}#L{line}",
repo="https://github.com/invenia/BlockDiagonals.jl/blob/{commit}{path}#{line}",
sitename="BlockDiagonals.jl",
authors="Invenia Technical Computing",
strict=true,
checkdocs=:exports,
)

deploydocs(;
repo="github.com/invenia/BlockDiagonals.jl",
push_preview=true,
)
5 changes: 2 additions & 3 deletions docs/src/index.md
Original file line number Diff line number Diff line change
@@ -1,11 +1,10 @@
# BlockDiagonals.jl

[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://invenia.github.io/BlockDiagonals.jl/stable)
[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://invenia.github.io/BlockDiagonals.jl/dev)
![CI](https://github.com/invenia/BlockDiagonals.jl/workflows/CI/badge.svg)
[![CI](https://github.com/invenia/BlockDiagonals.jl/workflows/CI/badge.svg)](https://github.com/Invenia/BlockDiagonals.jl/actions?query=workflow:CI)
[![Codecov](https://codecov.io/gh/invenia/BlockDiagonals.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/invenia/BlockDiagonals.jl)

Functionality for working efficiently with [block diagonal matrices](https://en.wikipedia.org/wiki/Block_matrix#Block_diagonal_matrices).
[BlockDiagonals.jl](https://github.com/invenia/BlockDiagonals.jl) provides functionality for working efficiently with [block diagonal matrices](https://en.wikipedia.org/wiki/Block_matrix#Block_diagonal_matrices).

```@autodocs
Modules = [BlockDiagonals]
Expand Down
27 changes: 16 additions & 11 deletions src/blockdiagonal.jl
Original file line number Diff line number Diff line change
@@ -1,19 +1,22 @@
# Core functionality for the `BlockDiagonal` type

"""
BlockDiagonal{T, V<:AbstractMatrix{T}} <: AbstractMatrix{T}
BlockDiagonal{T, V} <: AbstractMatrix{T}
BlockDiagonal(blocks::V) -> BlockDiagonal{T,V}

A matrix with matrices on the diagonal, and zeros off the diagonal.
"""
struct BlockDiagonal{T, V<:AbstractMatrix{T}} <: AbstractMatrix{T}
blocks::Vector{V}

function BlockDiagonal{T, V}(blocks::Vector{V}) where {T, V<:AbstractMatrix{T}}
return new{T, V}(blocks)
end
!!! info "`V` type"
`blocks::V` should be a `Tuple` or `AbstractVector` where each component (each block) is
`<:AbstractMatrix{T}` for some common element type `T`.
Copy link
Collaborator

Choose a reason for hiding this comment

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

do we want to constrain T to at least a Number, or are we happy with anything?

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 think if T isn't Number-like, some methods will fail. But 🤷 Would there be any gain from constraining things?

"""
struct BlockDiagonal{T, V} <: AbstractMatrix{T}
blocks::V
end

function BlockDiagonal(blocks::Vector{V}) where {T, V<:AbstractMatrix{T}}
function BlockDiagonal(blocks::V) where {
T, V<:Union{Tuple{Vararg{<:AbstractMatrix{T}}}, AbstractVector{<:AbstractMatrix{T}}}
Copy link
Contributor

Choose a reason for hiding this comment

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

Do you really want all block matrices to have the same eltype? Or could you live with them having different eltypes, and T being promote_type(map(eltype, blocks)...)? See LinearMaps.jl. If you need them to be homogeneous, then you may wish to include some promotion mechanism.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

could you live with them having different eltypes, and T being promote_type(map(eltype, blocks)...)?

that's a great suggestion

to be honest, i hadn't considered it. This currently just keeps T with the same meaning as on master.
Are there any downsides to T being promote_type(map(eltype, blocks)...)?

Copy link
Contributor Author

@nickrobinson251 nickrobinson251 Feb 16, 2021

Choose a reason for hiding this comment

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

trying this locally, and running tests, it does cause some inference issues for rrule

BlockDiagonal * Vector: Error During Test at /Users/npr/repos/BlockDiagonals.jl/test/chainrules.jl:12
  Got exception outside of a @test
  return type Array{Float64,1} does not match inferred return type Any
  Stacktrace:
   [1] error(::String) at ./error.jl:33
   [2] _test_inferred(::Function, ::InplaceableThunk{Thunk{BlockDiagonals.var"#13#19"{Array{Float64,1},BlockDiagonal{Float64,Array{Array{Float64,2},1}}}},BlockDiagonals.var"#14#20"{Array{Float64,1},BlockDiagonal{Float64,Array{Array{Float64,2},1}}}}; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /Users/npr/.julia/packages/ChainRulesTestUtils/2QRer/src/testers.jl:228
   [3] _test_inferred(::Function, ::InplaceableThunk{Thunk{BlockDiagonals.var"#13#19"{Array{Float64,1},BlockDiagonal{Float64,Array{Array{Float64,2},1}}}},BlockDiagonals.var"#14#20"{Array{Float64,1},BlockDiagonal{Float64,Array{Array{Float64,2},1}}}}) at /Users/npr/.julia/packages/ChainRulesTestUtils/2QRer/src/testers.jl:227
   [4] test_rrule(::typeof(*), ::BlockDiagonal{Float64,Array{Array{Float64,2},1}}, ::Vararg{Any,N} where N; output_tangent::ChainRulesTestUtils.Auto, fdm::FiniteDifferences.AdaptedFiniteDifferenceMethod{5,1,FiniteDifferences.UnadaptedFiniteDifferenceMethod{7,5}}, check_inferred::Bool, fkwargs::NamedTuple{(),Tuple{}}, rtol::Float64, atol::Float64, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /Users/npr/.julia/packages/ChainRulesTestUtils/2QRer/src/testers.jl:187
   [5] test_rrule(::Function, ::BlockDiagonal{Float64,Array{Array{Float64,2},1}}, ::Vararg{Any,N} where N) at /Users/npr/.julia/packages/ChainRulesTestUtils/2QRer/src/testers.jl:153
   [6] top-level scope at /Users/npr/repos/BlockDiagonals.jl/test/chainrules.jl:15

perhaps someone more familiar with ChainRules has some insight on why this might be (@willtebbutt , @mzgubic)?

Copy link
Collaborator

Choose a reason for hiding this comment

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

I don't know why this is (I assume the answer contains the word compiler). But I'd rather live with the same eltype than with this I think

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Okay, this change would be non-breaking anyway, so let's leave it as a potentially follow-up #62

Copy link
Contributor

Choose a reason for hiding this comment

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

I think that's perfectly fine. But then you may wish to consider to include a constructor that by itself does not yet require equal eltypes, but first determines T via promotion and then map(B -> convert(AbstractMatrix{T}, B), blocks). This should be a no-op if all matrices have the same eltype, but otherwise relieve you from annoying method errors. Just a suggestion.

}
return BlockDiagonal{T, V}(blocks)
end

Expand All @@ -22,15 +25,15 @@ BlockDiagonal(B::BlockDiagonal) = B
is_square(A::AbstractMatrix) = size(A, 1) == size(A, 2)

"""
blocks(B::BlockDiagonal{T, V}) -> Vector{V}
blocks(B::BlockDiagonal{T, V}) -> V

Return the on-diagonal blocks of B.
"""
blocks(B::BlockDiagonal) = B.blocks

# BlockArrays-like functions
"""
blocksizes(B::BlockDiagonal) -> Vector{Tuple}
blocksizes(B::BlockDiagonal{T, V}) -> V
Copy link
Collaborator

Choose a reason for hiding this comment

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

Do we want this to return a Vector or a Tuple depending on the block type? Or should we force it to return a vector/tuple always?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

my thinking is that we'd use blocks::Tuple for performance reasons when you've a small number of heterogenous blocks, in which case getting small Tuples from methods that operate on blocks will also be slightly better for performance, so no reason to force a Vector to be returned. It doesn't make sense to force a Tuple return in general, since there could be many elements.


Return the size of each on-diagonal block in order.

Expand Down Expand Up @@ -151,7 +154,9 @@ function _block_indices(B::BlockDiagonal, i::Integer, j::Integer)
p += 1
j -= ncols[p]
end
i -= sum(nrows[1:(p-1)])
if !isempty(nrows[1:(p-1)])
i -= sum(nrows[1:(p-1)])
end
# if row `i` outside of block `p`, set `p` to place-holder value `-1`
if i <= 0 || i > nrows[p]
p = -1
Expand Down
4 changes: 2 additions & 2 deletions src/chainrules.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# constructor
function ChainRulesCore.rrule(::Type{<:BlockDiagonal}, blocks::Vector{V}) where {V}
function ChainRulesCore.rrule(::Type{<:BlockDiagonal}, blocks::V) where {V<:AbstractVector}
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 think AbstractVector is fine here, but i'll rely on @willtebbutt to tell me if not

Copy link
Collaborator

Choose a reason for hiding this comment

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

If we exclude the tuple we won't have an rrule for constructors of BlockDiagonals with tuple blocks. Are you asking if it's ok not to have rrule for tuple blocks, or am I missing the point of the question?

Copy link
Contributor Author

@nickrobinson251 nickrobinson251 Feb 16, 2021

Choose a reason for hiding this comment

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

we won't have an rrule for constructors of BlockDiagonals with tuple blocks

correct, that'd be a nice feature to add as a follow-up, but out-of-scope for this PR

But my question is a different one. Before this PR, blocks were always a Vector, so we had rrule(::Type{<:BlockDiagonal}, blocks::Vector) whereas this PR allows them to be any AbstractVector (or Tuple but ignore that, since we're not adding that rrule here). So my question is, now that blocks can be any AbstractVector, are we okay with broadening the rrule from rrule(::Type{<:BlockDiagonal}, blocks::Vector) to rrule(::Type{<:BlockDiagonal}, blocks::AbstractVector) ?

(it's a question about how tightly to constrain rrules #56 (comment))

BlockDiagonal_pullback(Δ::Composite) = (NO_FIELDS, Δ.blocks)
return BlockDiagonal(blocks), BlockDiagonal_pullback
end
Expand Down Expand Up @@ -27,7 +27,7 @@ function ChainRulesCore.rrule(
::typeof(*),
bm::BlockDiagonal{T, V},
v::StridedVector{T}
) where {T<:Union{Real, Complex}, V<:Matrix{T}}
) where {T<:Union{Real, Complex}, V<:Vector{Matrix{T}}}
Copy link
Contributor Author

Choose a reason for hiding this comment

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

this okay, @willtebbutt ?


y = bm * v

Expand Down
13 changes: 12 additions & 1 deletion src/linalg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ svdvals_blockwise(B::BlockDiagonal) = mapreduce(svdvals, vcat, blocks(B))
LinearAlgebra.svdvals(B::BlockDiagonal) = sort!(svdvals_blockwise(B); rev=true)

# `B = U * Diagonal(S) * Vt` with `U` and `Vt` `BlockDiagonal` (`S` only sorted block-wise).
function svd_blockwise(B::BlockDiagonal{T}; full::Bool=false) where T
function svd_blockwise(B::BlockDiagonal{T, <:AbstractVector}; full::Bool=false) where T
U = Matrix{float(T)}[]
S = Vector{float(T)}()
Vt = Matrix{float(T)}[]
Expand All @@ -88,6 +88,17 @@ function svd_blockwise(B::BlockDiagonal{T}; full::Bool=false) where T
end
return BlockDiagonal(U), S, BlockDiagonal(Vt)
end
function svd_blockwise(B::BlockDiagonal{T, <:Tuple}; full::Bool=false) where T
S = Vector{float(T)}()
U_Vt = ntuple(length(blocks(B))) do i
F = svd(getblock(B, i), full=full)
append!(S, F.S)
(F.U, F.Vt)
Copy link
Collaborator

Choose a reason for hiding this comment

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

How come we aren't push!ing like above? Or vice versa if first.()/last.() is better?

Copy link
Contributor Author

@nickrobinson251 nickrobinson251 Feb 16, 2021

Choose a reason for hiding this comment

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

because above we want U (likewise Vt) to be a Vector, so create a Vector and add the elements. Whereas here we want it to be a Tuple hence ntuple, and first.(ntup) will still give us a tuple.

rather than push! we should probably be preallocating a Vector of the right length and overwriting the values, but the two functions would still be different (since Tuples are immutable)

end
U = first.(U_Vt)
Vt = last.(U_Vt)
return BlockDiagonal(U), S, BlockDiagonal(Vt)
end

function LinearAlgebra.svd(B::BlockDiagonal; full::Bool=false)::SVD
U, S, Vt = svd_blockwise(B, full=full)
Expand Down
Loading