Add support for SVD on BlockArray

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)))

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

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")

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