diff --git a/Project.toml b/Project.toml index 056f7f3..2a2d8b2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,13 +1,14 @@ name = "BlockSparseArrays" uuid = "2c9a651f-6452-4ace-a6ac-809f4280fbb4" authors = ["ITensor developers and contributors"] -version = "0.2.6" +version = "0.2.7" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" DerivableInterfaces = "6c5e35bf-e59e-4898-b73c-732dcc4ba65f" +DiagonalArrays = "74fd4be6-21e2-4f6f-823a-4360d37c7a77" Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" GradedUnitRanges = "e2de450a-8a67-46c7-b59c-01d5a3d041c5" @@ -31,7 +32,8 @@ Adapt = "4.1.1" Aqua = "0.8.9" ArrayLayouts = "1.10.4" BlockArrays = "1.2.0" -DerivableInterfaces = "0.3.7" +DerivableInterfaces = "0.3.8" +DiagonalArrays = "0.2.2" Dictionaries = "0.4.3" FillArrays = "1.13.0" GPUArraysCore = "0.1.0, 0.2" diff --git a/docs/make.jl b/docs/make.jl index d24e6c3..18b0932 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -2,7 +2,13 @@ using BlockSparseArrays: BlockSparseArrays using Documenter: Documenter, DocMeta, deploydocs, makedocs DocMeta.setdocmeta!( - BlockSparseArrays, :DocTestSetup, :(using BlockSparseArrays); recursive=true + BlockSparseArrays, + :DocTestSetup, + quote + using BlockSparseArrays + using LinearAlgebra: Diagonal + end; + recursive=true, ) include("make_index.jl") @@ -16,7 +22,7 @@ makedocs(; edit_link="main", assets=String[], ), - pages=["Home" => "index.md"], + pages=["Home" => "index.md", "Library" => "library.md"], ) deploydocs(; diff --git a/docs/src/library.md b/docs/src/library.md new file mode 100644 index 0000000..63af0e8 --- /dev/null +++ b/docs/src/library.md @@ -0,0 +1,5 @@ +# Library + +```@autodocs +Modules = [BlockSparseArrays] +``` diff --git a/src/BlockArraysExtensions/BlockArraysExtensions.jl b/src/BlockArraysExtensions/BlockArraysExtensions.jl index 3208508..e0a3c7a 100644 --- a/src/BlockArraysExtensions/BlockArraysExtensions.jl +++ b/src/BlockArraysExtensions/BlockArraysExtensions.jl @@ -598,3 +598,28 @@ macro view!(expr) @capture(expr, array_[indices__]) return :(view!($(esc(array)), $(esc.(indices)...))) end + +# SVD additions +# ------------- +using LinearAlgebra: Algorithm +using BlockArrays: BlockedMatrix + +# svd first calls `eigencopy_oftype` to create something that can be in-place SVD'd +# Here, we hijack this system to determine if there is any structure we can exploit +# default: SVD is most efficient with BlockedArray +function eigencopy_oftype(A::AbstractBlockArray, S) + return BlockedMatrix{S}(A) +end + +function svd!(A::BlockedMatrix; full::Bool=false, alg::Algorithm=default_svd_alg(A)) + F = svd!(parent(A); full, alg) + + # restore block pattern + m = length(F.S) + bax1, bax2, bax3 = axes(A, 1), blockedrange([m]), axes(A, 2) + + u = BlockedArray(F.U, (bax1, bax2)) + s = BlockedVector(F.S, (bax2,)) + vt = BlockedArray(F.Vt, (bax2, bax3)) + return SVD(u, s, vt) +end diff --git a/src/BlockSparseArrays.jl b/src/BlockSparseArrays.jl index 1425f1d..f125e2b 100644 --- a/src/BlockSparseArrays.jl +++ b/src/BlockSparseArrays.jl @@ -7,7 +7,13 @@ export BlockSparseArray, eachblockstoredindex, eachstoredblock +# factorizations +include("factorizations/svd.jl") + +# possible upstream contributions include("BlockArraysExtensions/BlockArraysExtensions.jl") + +# interface functions that don't have to specialize include("blocksparsearrayinterface/blocksparsearrayinterface.jl") include("blocksparsearrayinterface/linearalgebra.jl") include("blocksparsearrayinterface/getunstoredblock.jl") @@ -16,6 +22,8 @@ include("blocksparsearrayinterface/map.jl") include("blocksparsearrayinterface/arraylayouts.jl") include("blocksparsearrayinterface/views.jl") include("blocksparsearrayinterface/cat.jl") + +# functions defined for any abstractblocksparsearray include("abstractblocksparsearray/abstractblocksparsearray.jl") include("abstractblocksparsearray/wrappedabstractblocksparsearray.jl") include("abstractblocksparsearray/abstractblocksparsematrix.jl") @@ -27,8 +35,12 @@ include("abstractblocksparsearray/broadcast.jl") include("abstractblocksparsearray/map.jl") include("abstractblocksparsearray/linearalgebra.jl") include("abstractblocksparsearray/cat.jl") + +# functions specifically for BlockSparseArray include("blocksparsearray/defaults.jl") include("blocksparsearray/blocksparsearray.jl") +include("blocksparsearray/blockdiagonalarray.jl") + include("BlockArraysSparseArraysBaseExt/BlockArraysSparseArraysBaseExt.jl") end diff --git a/src/abstractblocksparsearray/abstractblocksparsematrix.jl b/src/abstractblocksparsearray/abstractblocksparsematrix.jl index 0c2c578..376f408 100644 --- a/src/abstractblocksparsearray/abstractblocksparsematrix.jl +++ b/src/abstractblocksparsearray/abstractblocksparsematrix.jl @@ -1 +1,84 @@ const AbstractBlockSparseMatrix{T} = AbstractBlockSparseArray{T,2} + +# SVD is implemented by trying to +# 1. Attempt to find a block-diagonal implementation by permuting +# 2. Fallback to AbstractBlockArray implementation via BlockedArray + +function eigencopy_oftype(A::AbstractBlockSparseMatrix, T) + if is_block_permutation_matrix(A) + Acopy = similar(A, T) + for bI in eachblockstoredindex(A) + Acopy[bI] = eigencopy_oftype(@view!(A[bI]), T) + end + return Acopy + else + return BlockedMatrix{T}(A) + end +end + +function is_block_permutation_matrix(a::AbstractBlockSparseMatrix) + return allunique(first ∘ Tuple, eachblockstoredindex(a)) && + allunique(last ∘ Tuple, eachblockstoredindex(a)) +end + +function _allocate_svd_output(A::AbstractBlockSparseMatrix, full::Bool, ::Algorithm) + @assert !full "TODO" + bm, bn = blocksize(A) + bmn = min(bm, bn) + + brows = blocklengths(axes(A, 1)) + bcols = blocklengths(axes(A, 2)) + slengths = Vector{Int}(undef, bmn) + + # fill in values for blocks that are present + bIs = collect(eachblockstoredindex(A)) + browIs = Int.(first.(Tuple.(bIs))) + bcolIs = Int.(last.(Tuple.(bIs))) + for bI in eachblockstoredindex(A) + row, col = Int.(Tuple(bI)) + nrows = brows[row] + ncols = bcols[col] + slengths[col] = min(nrows, ncols) + end + + # fill in values for blocks that aren't present, pairing them in order of occurence + # this is a convention, which at least gives the expected results for blockdiagonal + emptyrows = setdiff(1:bm, browIs) + emptycols = setdiff(1:bn, bcolIs) + for (row, col) in zip(emptyrows, emptycols) + slengths[col] = min(brows[row], bcols[col]) + end + + s_axis = blockedrange(slengths) + U = similar(A, axes(A, 1), s_axis) + S = similar(A, real(eltype(A)), s_axis) + Vt = similar(A, s_axis, axes(A, 2)) + + # also fill in identities for blocks that aren't present + for (row, col) in zip(emptyrows, emptycols) + copyto!(@view!(U[Block(row, col)]), LinearAlgebra.I) + copyto!(@view!(Vt[Block(col, col)]), LinearAlgebra.I) + end + + return U, S, Vt +end + +function svd(A::AbstractBlockSparseMatrix; kwargs...) + return svd!(eigencopy_oftype(A, LinearAlgebra.eigtype(eltype(A))); kwargs...) +end + +function svd!( + A::AbstractBlockSparseMatrix; full::Bool=false, alg::Algorithm=default_svd_alg(A) +) + @assert is_block_permutation_matrix(A) "Cannot keep sparsity: use `svd` to convert to `BlockedMatrix" + U, S, Vt = _allocate_svd_output(A, full, alg) + for bI in eachblockstoredindex(A) + bUSV = svd!(@view!(A[bI]); full, alg) + brow, bcol = Tuple(bI) + U[brow, bcol] = bUSV.U + S[bcol] = bUSV.S + Vt[bcol, bcol] = bUSV.Vt + end + + return SVD(U, S, Vt) +end diff --git a/src/blocksparsearray/blockdiagonalarray.jl b/src/blocksparsearray/blockdiagonalarray.jl new file mode 100644 index 0000000..7d0f281 --- /dev/null +++ b/src/blocksparsearray/blockdiagonalarray.jl @@ -0,0 +1,26 @@ +# type alias for block-diagonal +using LinearAlgebra: Diagonal +using DiagonalArrays: DiagonalArrays, diagonal + +const BlockDiagonal{T,A,Axes,V<:AbstractVector{A}} = BlockSparseMatrix{ + T,A,Diagonal{A,V},Axes +} +const BlockSparseDiagonal{T,A<:AbstractBlockSparseVector{T}} = Diagonal{T,A} + +@interface interface::BlockSparseArrayInterface function blocks(a::BlockSparseDiagonal) + return Diagonal(Diagonal.(blocks(a.diag))) +end + +function BlockDiagonal(blocks::AbstractVector{<:AbstractMatrix}) + return BlockSparseArray( + Diagonal(blocks), (blockedrange(size.(blocks, 1)), blockedrange(size.(blocks, 2))) + ) +end + +function DiagonalArrays.diagonal(S::BlockSparseVector) + D = similar(S, (axes(S, 1), axes(S, 1))) + for bI in eachblockstoredindex(S) + D[bI, bI] = diagonal(@view!(S[bI])) + end + return D +end diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl new file mode 100644 index 0000000..3c82f23 --- /dev/null +++ b/src/factorizations/svd.jl @@ -0,0 +1,206 @@ +using LinearAlgebra: + LinearAlgebra, Factorization, Algorithm, default_svd_alg, Adjoint, Transpose +using BlockArrays: AbstractBlockMatrix, BlockedArray, BlockedMatrix, BlockedVector +using BlockArrays: BlockLayout + +# Singular Value Decomposition: +# need new type to deal with U and V having possible different types +# this is basically a carbon copy of the LinearAlgebra implementation. +# additionally, by default we implement a fallback to the LinearAlgebra implementation +# in hope to support as many foreign types as possible that chose to extend those methods. + +# TODO: add this to MatrixFactorizations +# TODO: decide where this goes +# TODO: decide whether or not to restrict types to be blocked. +""" + SVD <: Factorization + +Matrix factorization type of the singular value decomposition (SVD) of a matrix `A`. +This is the return type of [`svd(_)`](@ref), the corresponding matrix factorization function. + +If `F::SVD` is the factorization object, `U`, `S`, `V` and `Vt` can be obtained +via `F.U`, `F.S`, `F.V` and `F.Vt`, such that `A = U * Diagonal(S) * Vt`. +The singular values in `S` are sorted in descending order. + +Iterating the decomposition produces the components `U`, `S`, and `V`. + +# Examples +```jldoctest +julia> A = [1. 0. 0. 0. 2.; 0. 0. 3. 0. 0.; 0. 0. 0. 0. 0.; 0. 2. 0. 0. 0.] +4×5 Matrix{Float64}: + 1.0 0.0 0.0 0.0 2.0 + 0.0 0.0 3.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 + 0.0 2.0 0.0 0.0 0.0 + +julia> F = BlockSparseArrays.svd(A) +BlockSparseArrays.SVD{Float64, Float64, Matrix{Float64}, Vector{Float64}, Matrix{Float64}} +U factor: +4×4 Matrix{Float64}: + 0.0 1.0 0.0 0.0 + 1.0 0.0 0.0 0.0 + 0.0 0.0 0.0 1.0 + 0.0 0.0 -1.0 0.0 +singular values: +4-element Vector{Float64}: + 3.0 + 2.23606797749979 + 2.0 + 0.0 +Vt factor: +4×5 Matrix{Float64}: + -0.0 0.0 1.0 -0.0 0.0 + 0.447214 0.0 0.0 0.0 0.894427 + 0.0 -1.0 0.0 0.0 0.0 + 0.0 0.0 0.0 1.0 0.0 + +julia> F.U * Diagonal(F.S) * F.Vt +4×5 Matrix{Float64}: + 1.0 0.0 0.0 0.0 2.0 + 0.0 0.0 3.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 + 0.0 2.0 0.0 0.0 0.0 + +julia> u, s, v = F; # destructuring via iteration + +julia> u == F.U && s == F.S && v == F.V +true +``` +""" +struct SVD{T,Tr,M<:AbstractArray{T},C<:AbstractVector{Tr},N<:AbstractArray{T}} <: + Factorization{T} + U::M + S::C + Vt::N + function SVD{T,Tr,M,C,N}( + U, S, Vt + ) where {T,Tr,M<:AbstractArray{T},C<:AbstractVector{Tr},N<:AbstractArray{T}} + Base.require_one_based_indexing(U, S, Vt) + return new{T,Tr,M,C,N}(U, S, Vt) + end +end +function SVD(U::AbstractArray{T}, S::AbstractVector{Tr}, Vt::AbstractArray{T}) where {T,Tr} + return SVD{T,Tr,typeof(U),typeof(S),typeof(Vt)}(U, S, Vt) +end +function SVD{T}(U::AbstractArray, S::AbstractVector{Tr}, Vt::AbstractArray) where {T,Tr} + return SVD( + convert(AbstractArray{T}, U), + convert(AbstractVector{Tr}, S), + convert(AbstractArray{T}, Vt), + ) +end + +function SVD{T}(F::SVD) where {T} + return SVD( + convert(AbstractMatrix{T}, F.U), + convert(AbstractVector{real(T)}, F.S), + convert(AbstractMatrix{T}, F.Vt), + ) +end +LinearAlgebra.Factorization{T}(F::SVD) where {T} = SVD{T}(F) + +# iteration for destructuring into components +Base.iterate(S::SVD) = (S.U, Val(:S)) +Base.iterate(S::SVD, ::Val{:S}) = (S.S, Val(:V)) +Base.iterate(S::SVD, ::Val{:V}) = (S.V, Val(:done)) +Base.iterate(::SVD, ::Val{:done}) = nothing + +function Base.getproperty(F::SVD, d::Symbol) + if d === :V + return getfield(F, :Vt)' + else + return getfield(F, d) + end +end + +function Base.propertynames(F::SVD, private::Bool=false) + return private ? (:V, fieldnames(typeof(F))...) : (:U, :S, :V, :Vt) +end + +Base.size(A::SVD, dim::Integer) = dim == 1 ? size(A.U, dim) : size(A.Vt, dim) +Base.size(A::SVD) = (size(A, 1), size(A, 2)) + +function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, F::SVD) + summary(io, F) + println(io) + println(io, "U factor:") + show(io, mime, F.U) + println(io, "\nsingular values:") + show(io, mime, F.S) + println(io, "\nVt factor:") + return show(io, mime, F.Vt) +end + +Base.adjoint(usv::SVD) = SVD(adjoint(usv.Vt), usv.S, adjoint(usv.U)) +Base.transpose(usv::SVD) = SVD(transpose(usv.Vt), usv.S, transpose(usv.U)) + +# Conversion +Base.AbstractMatrix(F::SVD) = (F.U * Diagonal(F.S)) * F.Vt +Base.AbstractArray(F::SVD) = AbstractMatrix(F) +Base.Matrix(F::SVD) = Array(AbstractArray(F)) +Base.Array(F::SVD) = Matrix(F) +SVD(usv::SVD) = usv +SVD(usv::LinearAlgebra.SVD) = SVD(usv.U, usv.S, usv.Vt) + +# functions default to LinearAlgebra +# ---------------------------------- +""" + svd!(A; full::Bool = false, alg::Algorithm = default_svd_alg(A)) -> SVD + +`svd!` is the same as [`svd`](@ref), but saves space by +overwriting the input `A`, instead of creating a copy. See documentation of [`svd`](@ref) for details. +""" +svd!(A; kwargs...) = SVD(LinearAlgebra.svd!(A; kwargs...)) + +""" + svd(A; full::Bool = false, alg::Algorithm = default_svd_alg(A)) -> SVD + +Compute the singular value decomposition (SVD) of `A` and return an `SVD` object. + +`U`, `S`, `V` and `Vt` can be obtained from the factorization `F` with `F.U`, +`F.S`, `F.V` and `F.Vt`, such that `A = U * Diagonal(S) * Vt`. +The algorithm produces `Vt` and hence `Vt` is more efficient to extract than `V`. +The singular values in `S` are sorted in descending order. + +Iterating the decomposition produces the components `U`, `S`, and `V`. + +If `full = false` (default), a "thin" SVD is returned. For an ``M +\\times N`` matrix `A`, in the full factorization `U` is ``M \\times M`` +and `V` is ``N \\times N``, while in the thin factorization `U` is ``M +\\times K`` and `V` is ``N \\times K``, where ``K = \\min(M,N)`` is the +number of singular values. + +`alg` specifies which algorithm and LAPACK method to use for SVD: +- `alg = DivideAndConquer()` (default): Calls `LAPACK.gesdd!`. +- `alg = QRIteration()`: Calls `LAPACK.gesvd!` (typically slower but more accurate) . + +!!! compat "Julia 1.3" + The `alg` keyword argument requires Julia 1.3 or later. + +# Examples +```jldoctest +julia> A = rand(4,3); + +julia> F = BlockSparseArrays.svd(A); # Store the Factorization Object + +julia> A ≈ F.U * Diagonal(F.S) * F.Vt +true + +julia> U, S, V = F; # destructuring via iteration + +julia> A ≈ U * Diagonal(S) * V' +true + +julia> Uonly, = BlockSparseArrays.svd(A); # Store U only + +julia> Uonly == U +true +``` +""" +svd(A; kwargs...) = + SVD(svd!(eigencopy_oftype(A, LinearAlgebra.eigtype(eltype(A))); kwargs...)) + +LinearAlgebra.svdvals(usv::SVD{<:Any,T}) where {T} = (usv.S)::AbstractVector{T} + +# Added here to avoid type-piracy +eigencopy_oftype(A, S) = LinearAlgebra.eigencopy_oftype(A, S) diff --git a/test/Project.toml b/test/Project.toml index 4268ad3..0cc3624 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -4,6 +4,7 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" BlockSparseArrays = "2c9a651f-6452-4ace-a6ac-809f4280fbb4" +DiagonalArrays = "74fd4be6-21e2-4f6f-823a-4360d37c7a77" GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" GradedUnitRanges = "e2de450a-8a67-46c7-b59c-01d5a3d041c5" JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb" diff --git a/test/basics/test_svd.jl b/test/basics/test_svd.jl new file mode 100644 index 0000000..70dc984 --- /dev/null +++ b/test/basics/test_svd.jl @@ -0,0 +1,78 @@ +using Test +using BlockSparseArrays +using BlockSparseArrays: BlockSparseArray, svd, BlockDiagonal, eachblockstoredindex +using BlockArrays +using Random +using DiagonalArrays: diagonal +using LinearAlgebra: LinearAlgebra + +function test_svd(a, usv) + U, S, V = usv + + @test U * diagonal(S) * V' ≈ a + @test U' * U ≈ LinearAlgebra.I + @test V' * V ≈ LinearAlgebra.I +end + +# regular matrix +# -------------- +sizes = ((3, 3), (4, 3), (3, 4)) +eltypes = (Float32, Float64, ComplexF64) +@testset "($m, $n) Matrix{$T}" for ((m, n), T) in Iterators.product(sizes, eltypes) + a = rand(m, n) + usv = @inferred svd(a) + test_svd(a, usv) +end + +# block matrix +# ------------ +blockszs = (([2, 2], [2, 2]), ([2, 2], [2, 3]), ([2, 2, 1], [2, 3]), ([2, 3], [2])) +@testset "($m, $n) BlockMatrix{$T}" for ((m, n), T) in Iterators.product(blockszs, eltypes) + a = mortar([rand(T, i, j) for i in m, j in n]) + usv = svd(a) + test_svd(a, usv) + @test usv.U isa BlockedMatrix + @test usv.Vt isa BlockedMatrix + @test usv.S isa BlockedVector +end + +# Block-Diagonal matrices +# ----------------------- +@testset "($m, $n) BlockDiagonal{$T}" for ((m, n), T) in + Iterators.product(blockszs, eltypes) + a = BlockDiagonal([rand(T, i, j) for (i, j) in zip(m, n)]) + usv = svd(a) + # TODO: `BlockDiagonal * Adjoint` errors + test_svd(a, usv) +end + +# blocksparse +# ----------- +@testset "($m, $n) BlockSparseMatrix{$T}" for ((m, n), T) in + Iterators.product(blockszs, eltypes) + a = BlockSparseArray{T}(m, n) + + # test empty matrix + usv_empty = svd(a) + test_svd(a, usv_empty) + + # test blockdiagonal + for i in LinearAlgebra.diagind(blocks(a)) + I = CartesianIndices(blocks(a))[i] + a[Block(I.I...)] = rand(T, size(blocks(a)[i])) + end + usv = svd(a) + test_svd(a, usv) + + perm = Random.randperm(length(m)) + b = a[Block.(perm), Block.(1:length(n))] + usv = svd(b) + test_svd(b, usv) + + # test permuted blockdiagonal with missing row/col + I_removed = rand(eachblockstoredindex(b)) + c = copy(b) + delete!(blocks(c).storage, CartesianIndex(Int.(Tuple(I_removed)))) + usv = svd(c) + test_svd(c, usv) +end