Add SVD support for BlockArrays

Author: lkdvos

src/BlockArrays.jl

@@ -66,6 +66,7 @@ include("blocks.jl")
 
 include("blockbroadcast.jl")
 include("blockcholesky.jl")
+include("blocksvd.jl")
 include("blocklinalg.jl")
 include("blockproduct.jl")
 include("show.jl")

src/blocksvd.jl

@@ -0,0 +1,21 @@
+#=
+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`)
+=#
+
+LinearAlgebra.eigencopy_oftype(A::AbstractBlockMatrix, S) = BlockedArray(Array{S}(A), blocksizes(A, 1), blocksizes(A, 2))
+
+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

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,73 @@
+module TestBlockSVD
+
+using BlockArrays, Test, LinearAlgebra, Random
+
+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
+
+end # module