Add support for SVD on BlockArray

[lkdvos]

This PR implements support for LinearAlgebra.svd on block arrays, which tries to retain the block structure as much as possible, without making choices for structures that are not canonical.

In particular, this means that for U, S, V = svd(A), the rows of U should have the same structure as the rows of A, which allows U' * A to behave as expected, and similarly the rows of V should have the same structure as the columns of A to support A * V. There is no real block structure that carries over to the rows and columns of S however, resulting in a single block.

As far as I know, there is no real way of implementing this efficiently by making use of the block structure, so the implementation first maps the arrays to a BlockedArray, which is then used to perform the decomposition. Similarly, the resulting U, S and Vt are always BlockedArray, to reflect this.

For BlockDiagonal matrices however, these can be considered as linear maps that conserve the block structure, and in particular the SVD can efficiently be implemented block-by-block. For these cases, a specialized implementation is provided, which does carry over the block structure to S.

Questions and comments

• Because of how LinearAlgebra defines the SVD struct, U and Vt must have the same structure. In particular, this means that the original block matrix should have the same type of blocksizes for its rows and columns. There is currently nothing handles this, or attempts to promote the types when they are not the same.
• There is a type instability related to Type instability of blocksizes(y, i) #425
• For BlockDiagonal matrices, I think it is quite natural that U and V are BlockDiagonal, but I am not so sure how to define S. Currently, it is a BlockVector, which is in line with how LinearAlgebra requires this to be a vector, but it does make it a bit cumbersome to work with: U * Diagonal(S) * Vt does not actually work, as BlockArray(Diagonal) is not the same as Diagonal(BlockVector). It might be reasonable to specifically define Diagonal(::BlockVector) to return a BlockDiagonal of Diagonals
• It might be better to think a bit more carefully about how to define eigencopy_oftype, to avoid unnecessary copies, and to improve compatibility with GPU arrays.
• If this passes comments and checks, I'm happy to implement similar functionality for the other factorizations.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 93.73%. Comparing base (1e5feaa) to head (7ffbf03).

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #426      +/-   ##
==========================================
+ Coverage   93.67%   93.73%   +0.06%     
==========================================
  Files          18       19       +1     
  Lines        1643     1661      +18     
==========================================
+ Hits         1539     1557      +18     
  Misses        104      104              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[mtfishman]

Because of how LinearAlgebra defines the SVD struct, U and Vt must have the same structure. In particular, this means that the original block matrix should have the same type of blocksizes for its rows and columns. There is currently nothing handles this, or attempts to promote the types when they are not the same.

Maybe it would be better to define a BlockSVD struct that is like SVD but allows U and Vt to have different types. (@jishnub @dlfivefifty curious to hear your thoughts on that.)

[dlfivefifty]

Ive come across the issue of LinearAlgebra imposing different components have the same type many times. My solution has been to copy the code with an extra templated variable, see eg LazyBandedMatrices.Tridiagonal and LazyBandedMatrices.SymTridiagonal

If we follow this pattern we could create a BlockArrays.SVD (or possibly MatrixFactorizations.SVD) with an extra templated variable. To do this is actually pretty easy: copy the code and test from LinearAlgebra and just modify any of the templates

[mtfishman]

For BlockDiagonal matrices, I think it is quite natural that U and V are BlockDiagonal, but I am not so sure how to define S. Currently, it is a BlockVector, which is in line with how LinearAlgebra requires this to be a vector, but it does make it a bit cumbersome to work with: U * Diagonal(S) * Vt does not actually work, as BlockArray(Diagonal) is not the same as Diagonal(BlockVector). It might be reasonable to specifically define Diagonal(::BlockVector) to return a BlockDiagonal of Diagonals

Regarding this point, I don't think it is a good idea to make Diagonal(::BlockVector) return a BlockDiagonal since constructors should follow the rule/convention that T(args...) isa T. So probably it is better to have Diagonal(S) return a Diagonal but make sure it acts like a BlockDiagonal where the diagonal blocks are Diagonal. It acts like a blocked array thanks to:

BlockArrays.jl/src/blockarrayinterface.jl

Lines 14 to 16 in 1e5feaa

	function axes(D::Diagonal{<:Any,<:AbstractBlockVector})
	a = axes(parent(D),1)
	(a,a)

but I noticed that there may be some issues with that wrapping: #427. EDIT: I also reproduced that issue multiplying Diagonal wrapping an AbstractBlockVector with a BlockDiagonal, I raised an issue here: #428. Probably best to fix that in conjunction with this PR.