SVDS of a LinearMap

[dkarrasch]

This is intended as a reminder and progress tracking on the issue of allowing LinearMaps in svds, see 21701 and 27480, completely analogously to eigs. It appears to be feasible, even though the last comment in the second issue indicated that there is need for a *(::LinearMap,::AbstractMatrix)?

[dkarrasch]

Thanks, @andreasnoack . Testing svds for a FunctionMap <: LinearMap seems to go through the eigs call, but then fails upon calling qr on the singular vector "matrices" in lines 331 and 334. There are two options:

1. I define rmul!(A::LinearMap,D::Diagonal) via Matrix(A * D), where A * D is a composed LinearMap and Matrix is like the former full, so retrieves a matrix representation. This is fine for usual SDV applications, where D is usually small, but it seems bad as a generic definition. The purpose of the LinearMap package is to avoid computing matrix representations.
2. You call Matrix on rmul!(...) in the above mentioned lines, within the qr call. That should be a no-op for the usual matrix cases. In any case, this is not part of some iterative process, just the single time post-processing of the singular vectors.

I'd prefer very much option 2. Would that be acceptible for Arpack.jl?

[andreasnoack]

I'm not sure I understand. In

Arpack.jl/src/Arpack.jl

Lines 330 to 331 in 2549905

	V = ex[2]
	U = Array(qr(rmul!(X*V, Diagonal(inv.(svals)))).Q)

V is a Matrix so why wouldn't X::FunctionMap*V also be a Matrix?

[andreasnoack]

Btw, if you don't mind, I would like to hear how https://github.com/andreasnoack/TSVD.jl compares.

[dkarrasch]

Because by default LinearMaps are defined by their action on vectors ("only"). The product of a LinearMap with a Matrix is a CompositeMap, interpreting the matrix as a linear operator, and not as a block of vectors, where you would apply X to each column of V. Imagine X and V to be of large size (generally, not in the usual SVDS case), then that is expected to be expensive. For when you want this to be interpreted as a block of column vectors and want X*V be a Matrix, that's exactly what we have the Matrix constructor for. I'd hope that having the Matrix call here instead of in rmul!(::LinearMap,::Diagonal) would hurt less (if at all). Or are there type-issues when calling Matrix (or Array, which is both defined for LinearMaps).

Sure, I'll take a look at the TSVD package.

[dkarrasch]

Re TSVD.jl: On the LinearMaps side, we would need instances for
A_mul_B!(α::Number, A, x::Vector, β::Number, y::Vector) and Ac_mul_B!(α::Number, A, x::Vector, β::Number, y::Vector) (or Julia 0.7 compatible versions thereof), which we currently don't have. I assume that the mutated vector there is y?

Re Arpack.svds: I realized I cannot define rmul! as a mutating function if rmul!(::LinearMap,::Matrix)::LinearMap because of the additional composition wrapper. So the only solution that comes to my mind is to change the definition to *(::LinearMap,::Matrix)::Matrix, which is not ideal as explained above. Or we could recommend using TSVD or IterativeSolvers.jl.

[andreasnoack]

I'm not sure I agree that it is the right default not to materialize the result of LinearMap*Matrix. When you work with a Matrix, you have implicitly said that it is not too large to handle. Sure (A*B)*x will be more costly than A*(B*x) but if we go down that route then Matrix*Matrix also shouldn't materialize. Anyway, I'm fine changing this to convert(Matrix, X*V) since it is a noop in the normal case. Do you have the convert method defined or only the constructor?

Regarding TSVD then I'd forgotten that it needs to be updated to 0.7. Hopefully, I'll get around to do that soon.

[dkarrasch]

We currently only have an Array and a Matrix constructor. How would we write define convert?

I guess it is a design decision, so there is no right or wrong. A LinearMaps.jl user proposed to distinguish the two cases by having *(::LinearMap,::Matrix) materialize and *(::LinearMap,::LinearMap) compose. I think that should be possible, but the change is "breaking" in a sense, because output types will change. Anyway, there is the additional issue with rmul! that would prevent us from being able to use Arpack.svds for LinearMaps. So, there is no need to change anything right now. Many thanks for your consideration.

[andreasnoack]

In your case, the convert method should just call the constructor, i.e. convert(::Type{Matrix}, A:: LinearMap) = Matrix(A)

Anyway, there is the additional issue with rmul!

What would that be? Once, X*V has been converted to a Matrix I can't see how LinearMaps can influence the computation.

[andreasnoack]

#23 Should hopefully be sufficient once you've defined the convert method.