diff --git a/src/BlockArrays.jl b/src/BlockArrays.jl index 58fa689b..4c8c9aea 100644 --- a/src/BlockArrays.jl +++ b/src/BlockArrays.jl @@ -73,6 +73,7 @@ include("blockreduce.jl") include("blockdeque.jl") include("blockarrayinterface.jl") include("blockbanded.jl") +include("blocksvd.jl") @deprecate getblock(A::AbstractBlockArray{T,N}, I::Vararg{Integer, N}) where {T,N} view(A, Block(I)) @deprecate getblock!(X, A::AbstractBlockArray{T,N}, I::Vararg{Integer, N}) where {T,N} copyto!(X, view(A, Block(I))) diff --git a/src/blocksvd.jl b/src/blocksvd.jl new file mode 100644 index 00000000..4b964431 --- /dev/null +++ b/src/blocksvd.jl @@ -0,0 +1,35 @@ +#= +SVD on blockmatrices: +Interpret the matrix as a linear map acting on vector spaces with a direct sum structure, which is reflected in the structure of U and V. +In the generic case, the SVD does not preserve this structure, and can mix up the blocks, so S becomes a single block. +(Implementation-wise, also most efficiently done by first mapping to a `BlockedArray`) +In the case of `BlockDiagonal` however, the structure is preserved and carried over to the structure of `S`. +=# + +LinearAlgebra.eigencopy_oftype(A::AbstractBlockMatrix, S) = BlockedArray{S}(A) + +function LinearAlgebra.eigencopy_oftype(A::BlockDiagonal, S) + diag = map(Base.Fix2(LinearAlgebra.eigencopy_oftype, S), A.blocks.diag) + return BlockDiagonal(diag) +end + +function LinearAlgebra.svd!(A::BlockedMatrix; full::Bool=false, alg::LinearAlgebra.Algorithm=default_svd_alg(A)) + USV = svd!(parent(A); full, alg) + + # restore block pattern + m = length(USV.S) + bsz1, bsz2, bsz3 = blocksizes(A, 1), [m], blocksizes(A, 2) + + u = BlockedArray(USV.U, bsz1, bsz2) + s = BlockedVector(USV.S, bsz2) + vt = BlockedArray(USV.Vt, bsz2, bsz3) + return SVD(u, s, vt) +end + +function LinearAlgebra.svd!(A::BlockDiagonal; full::Bool=false, alg::LinearAlgebra.Algorithm=default_svd_alg(A)) + USVs = map(a -> svd!(a; full, alg), A.blocks.diag) + Us = map(Base.Fix2(getproperty, :U), USVs) + Ss = map(Base.Fix2(getproperty, :S), USVs) + Vts = map(Base.Fix2(getproperty, :Vt), USVs) + return SVD(BlockDiagonal(Us), mortar(Ss), BlockDiagonal(Vts)) +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index f6ce9823..cfee2ded 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -22,6 +22,7 @@ include("test_blockrange.jl") include("test_blockarrayinterface.jl") include("test_blockbroadcast.jl") include("test_blocklinalg.jl") +include("test_blocksvd.jl") include("test_blockproduct.jl") include("test_blockreduce.jl") include("test_blockdeque.jl") diff --git a/test/test_blocksvd.jl b/test/test_blocksvd.jl new file mode 100644 index 00000000..cc86b7a5 --- /dev/null +++ b/test/test_blocksvd.jl @@ -0,0 +1,107 @@ +module TestBlockSVD + +using BlockArrays, Test, LinearAlgebra, Random +using BlockArrays: BlockDiagonal + +Random.seed!(0) + +eltypes = (Float32, Float64, ComplexF32, ComplexF64, Int) + +@testset "Block SVD ($T)" for T in eltypes + x = rand(T, 4, 4) + USV = svd(x) + U, S, Vt = USV.U, USV.S, USV.Vt + + y = BlockArray(x, [2, 2], [2, 2]) + # https://github.com/JuliaArrays/BlockArrays.jl/issues/425 + # USV_blocked = @inferred svd(y) + USV_block = svd(y) + U_block, S_block, Vt_block = USV_block.U, USV_block.S, USV_block.Vt + + # test types + @test U_block isa BlockedMatrix + @test eltype(U_block) == float(T) + @test S_block isa BlockedVector + @test eltype(S_block) == real(float(T)) + @test Vt_block isa BlockedMatrix + @test eltype(Vt_block) == float(T) + + # test structure + @test blocksizes(U_block, 1) == blocksizes(y, 1) + @test length(blocksizes(U_block, 2)) == 1 + @test blocksizes(Vt_block, 2) == blocksizes(y, 2) + @test length(blocksizes(Vt_block, 1)) == 1 + + # test correctness + @test U_block ≈ U + @test S_block ≈ S + @test Vt_block ≈ Vt + @test U_block * Diagonal(S_block) * Vt_block ≈ y +end + +@testset "Blocked SVD ($T)" for T in eltypes + x = rand(T, 4, 4) + USV = svd(x) + U, S, Vt = USV.U, USV.S, USV.Vt + + y = BlockedArray(x, [2, 2], [2, 2]) + # https://github.com/JuliaArrays/BlockArrays.jl/issues/425 + # USV_blocked = @inferred svd(y) + USV_blocked = svd(y) + U_blocked, S_blocked, Vt_blocked = USV_blocked.U, USV_blocked.S, USV_blocked.Vt + + # test types + @test U_blocked isa BlockedMatrix + @test eltype(U_blocked) == float(T) + @test S_blocked isa BlockedVector + @test eltype(S_blocked) == real(float(T)) + @test Vt_blocked isa BlockedMatrix + @test eltype(Vt_blocked) == float(T) + + # test structure + @test blocksizes(U_blocked, 1) == blocksizes(y, 1) + @test length(blocksizes(U_blocked, 2)) == 1 + @test blocksizes(Vt_blocked, 2) == blocksizes(y, 2) + @test length(blocksizes(Vt_blocked, 1)) == 1 + + # test correctness + @test U_blocked ≈ U + @test S_blocked ≈ S + @test Vt_blocked ≈ Vt + @test U_blocked * Diagonal(S_blocked) * Vt_blocked ≈ y +end + +@testset "BlockDiagonal SVD ($T)" for T in eltypes + blocksz = (2, 3, 1) + y = BlockDiagonal([rand(T, d, d) for d in blocksz]) + x = Array(y) + + USV = svd(x) + U, S, Vt = USV.U, USV.S, USV.Vt + + # https://github.com/JuliaArrays/BlockArrays.jl/issues/425 + # USV_blocked = @inferred svd(y) + USV_block = svd(y) + U_block, S_block, Vt_block = USV_block.U, USV_block.S, USV_block.Vt + + # test types + @test U_block isa BlockDiagonal + @test eltype(U_block) == float(T) + @test S_block isa BlockVector + @test eltype(S_block) == real(float(T)) + @test Vt_block isa BlockDiagonal + @test eltype(Vt_block) == float(T) + + # test structure + @test blocksizes(U_block, 1) == blocksizes(y, 1) + @test length(blocksizes(U_block, 2)) == length(blocksz) + @test blocksizes(Vt_block, 2) == blocksizes(y, 2) + @test length(blocksizes(Vt_block, 1)) == length(blocksz) + + # test correctness: SVD is not unique, so cannot compare to dense + @test U_block * BlockDiagonal(Diagonal.(S_block.blocks)) * Vt_block ≈ y + @test U_block' * U_block ≈ LinearAlgebra.I + @test Vt_block * Vt_block' ≈ LinearAlgebra.I +end + +end # module