Subtype AbstractArray? - Githubissues

I was thinking it might make sense to have operators subtype abstract array and implement the corresponding interface. I think that would make sense conceptually, and would make it easier to use other Julia functions with them. For example, one could then use QuantumOptics Operators in LazyArrays.jl functions immediately, and not have to worry about implementing lazy functions in QuantumOptics too.

Does that seem like a good idea?

This definitely sounds interesting. However, I don't think it's as simple as it first sounds and I am not sure what the best strategy here is.

To elaborate: there are quite a few different Operator subtypes, and not all of them should be subtyped to AbstractArray. The most straightforward implementation would be for DenseOperator. For SparseOperator I think we essentially need to define our own version of SparseMatrixCSC (unless I am missing something). This should not be too hard though.

The biggest problem, imo, would be the FFTOperator type. This is an Operator in the sense that it maps a vector on a Hilbert space to another one, but it is not represented by an array and I don't think it should be a subtype of AbstractArray. We could in principle separate this type from all others and define the appropriate methods, but I am not sure if this is the right approach here.

Finally, there is also the Lazy types. As you point out, these could be replaced by using LazyArrays. I was not aware of this package, but it would be great to leave optimizations of lazy types to another package. Especially, if we decide to do the subtyping this would make a ton of sense. The problem I have with this, though, is that this would mean a pretty big regression in speed. Consider the following piece of code:

using LazyArrays, BenchmarkTools, QuantumOptics, LinearAlgebra

N = 50
A = randn(ComplexF64, N,N); B = randn(ComplexF64, N,N)
C = Kron(A,B) # lazy kronecker product
b = randn(ComplexF64, N^2); c = similar(b);

bgen = GenericBasis(N)
A_op = DenseOperator(bgen, A)
B_op = DenseOperator(bgen, B)
C_op = LazyTensor(bgen^2, [1, 2], [B_op, A_op])
b_ket = Ket(bgen^2, b)
c_ket = copy(b_ket)

mul!(c, C, b)
operators.gemv!(complex(1.0), C_op, b_ket, 0.0im, c_ket)
@assert c ≈ c_ket.data

Now, if I am not missing something here, it seems that the currently implemented in-place multiplication in QuantumOptics is much faster:

julia> @benchmark mul!(c, C, b)
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     155.144 ms (0.00% GC)
  median time:      155.447 ms (0.00% GC)
  mean time:        155.668 ms (0.00% GC)
  maximum time:     157.180 ms (0.00% GC)
  --------------
  samples:          33
  evals/sample:     1

julia> @benchmark operators.gemv!(complex(1.0), C_op, b_ket, complex(0.0), c_ket)
BenchmarkTools.Trial:
  memory estimate:  592 bytes
  allocs estimate:  11
  --------------
  minimum time:     65.707 ms (0.00% GC)
  median time:      66.104 ms (0.00% GC)
  mean time:        66.711 ms (0.00% GC)
  maximum time:     79.897 ms (0.00% GC)
  --------------
  samples:          75
  evals/sample:     1

I would really like to discuss this idea more in order to figure out the best way to go about things here. I am just not sure how to best restructure the type system to do this (especially with the FFTOperator). As for the LazyArrays, we could just define our own methods for now, until there is better in-place multiplication.

Thanks for considering it! By the way, I'm pretty new to Julia and even to programming more generally so I'm sure there are more qualified people to discuss with.

For sparse operators, it looks like these's some discussion here about general interfaces for sparse arrays: https://github.com/JuliaLang/julia/issues/25760 which might be relevant (but it doesn't look like there was any resolution there). Also, I don't quite understand what you mean by "I think we essentially need to define our own version of SparseMatrixCSC." Couldn't the data field just be a SparseMatrixCSC, as it is now? Or do you mean that when we define getindex and etc for users to access the data directly from the operator object, we will essentially be copying the methods of SparseMatrixCSC?

For FFTOperator's before: I actually haven't used them so didn't think of that at all. It seems like one wouldn't want to define indexing with them, and the array framework seems inappropriate as you point out. I saw another package, LinearMaps.jl that seems to deal with similar situations: they define linear maps without indexing, that only need to have a defined matrix-vector product. Could that be a useful framework? They seem to use lazy methods for composition, but don't have a tensor product / kron feature.

A related aside: I think one other thing worth thinking about for the design of the operators, is that it seems like usually one would not mutate and operator to change the basis on which it acts, so the basis (and in particular the dimension) would be fixed. So the dimensional information could possibly be able to be known at compile-time. The StaticArrays.jl package makes array types that use this information to dramatically speed things up for low-dimensional objects (I can generate random 4x4 density matrices 8 times faster just by subtyping StaticArray). I think that might not be this packages' focus though (since it's called QuantumOptics afterall). I've been using QuantumOptics just to have a nice setup of operators and bases to do small-dimensional things, like tensor together a few qubits and solve an SDP, or look for a counterexample, etc, so I'm interested in this kind of thing. (Thanks for the great package!). But the fact that the dimension could be available at compile time could still be useful for large dimensions; for example, if you tried to multiply operators which are incompatible you could get a compile time error instead of a run time error.

Putting these things together, I thought Operators could be parametric types, parametrized by their bases, and subtypes of AbstractArray. Then the basis checks can happen at compile time. Also, the basis could determine the storage type of the operator. For example, maybe Fock space is associated to sparse operators, and Operators on a Fock space all have a sparse data-type. Then functions can dispatch on the type of basis and use the right (optimized) operations depending on that. So instead of defining f(::DenseOperator) and f(::SparseOperator), one could have f(::Operator{B}) where {T <: SparseBasis} or something like that. And I could make a StaticBasis or something which itself is parametrized by the dimension (so there's a different type for every dimension, allowing compile-time use of the dimension information), which presumably would mostly be useful for low dimensional things. I don't actually know if this is a good idea; maybe tying together the storage type and Hilbert space is unnecessary, or a bad idea.

Thanks for all the references, they are very helpful.

Yes, you are probably right in that we can just do it via the data field. Never mind my statement there. Also, in the reference you provided they point to the AbstractSparseArray type, which we should probably consider.

LinearMaps.jl is indeed precisely the same as the way we handle the FFTOperator (combined with LazyProduct). We could also consider using LinearMaps here, but it's just for the one type so I'm not sure. In any case, this is what I meant with keeping this type separate from the rest.

I have also though about using StaticArrays as I have rarely found myself wanting to change the dimension of an operator. Also, it's actually not so simple since the bases fields carry the dimension information. Therefore, it might actually make sense to use this as default for the data fields. Alternatively, we can make a new StaticOperator type which does this. If we don't make it the default, I would leave this as a later addition.

This type parametrerization you mention is subject of the PR I started #234. There, I did exactly what you suggest: namely create one Operator{BL,BR,T} type, where the first two are bases and the last one the data field type (either matrix or sparse array). This one type replaces DenseOperator and SparseOperator. I would definitely put the subtyping together with the changes that I started there, though I am no longer sure if it isn't better to leave sparse and dense operators as separate types (e.g. if we want to subtype SparseOperator to AbstractSparseArray).

I don't think the decision whether operators are sparse or dense should be given by the basis, though. For example, the position operator is diagonal in the position basis (thus should be stored as sparse), whereas the momentum operator has only nonzero elements. Any basis should allow for both operator types.

Also note, that these are all pretty essential changes to the type system. So once we decide on the strategy, it will still take some time to do it properly.

Ah, I saw that PR awhile ago but for some reason didn't look at it more carefully. Yeah, that's exactly what I was thinking. In fact, now that you point to that, I surely read the issue that points out the checking bases at compile time at some point-- that must be why that was in my head. Sorry. Anyway, I see what you mean with AbstractSparseMatrix though and splitting back to DenseOperators and SparseOperators, since if these were to subtype AbstractArray, one would want SparseOperator to subtype AbstractSparseMatrix.

One thing to mention though is the type relation AbstractSparseMatrix <: AbstractMatrix (source). So it could still make sense to have just one parametric type, Operator <: AbstractMatrix, and dispatch on sparse data structures differently, as is done for subtypes of AbstractSparseMatrix. Maybe then one has to repeat everything that's done for AbstractSparseMatrix's again for Operator{BL, BR, T} where {BR,BL, T<:SparseDataType}. Maybe that's a bad idea though.

Re: not tying the storage type to the basis-- I see, that makes sense.

For the StaticArrays considerations-- it seems like there are performance (or at least compile-time) penalties for large matrices. In the readme they mention "For example, the performance crossover point for a matrix multiply microbenchmark seems to be about 11x11 in julia 0.5 with default optimizations;" for my experience, I got impatient just waiting to evaluate one 20x20 random density matrix on my laptop using StaticArrays based code, and ended up cancelling it after a few minutes.

Maybe one could have DenseOperator{BL,BR,T} where T is either a StaticArray type or a usual complex matrix? (or even just a subtype of AbstractArray in general). I imagine most of the code wouldn't dispatch on T (since for the most part you can use them the same), but if something special needed to be defined for StaticArrays, it could be (for example, I would get errors on the current release of StaticArray when trying to take real powers of a matrix, so maybe workarounds would sometimes be needed for things like this). That could also open it down the line for other data storage types to be added without breaking things. (Also, this way StaticArray could be a later addition as you mention).

Okay, so to summarize: I like the idea of subtyping operators that are represented by matrices to AbstractArray and would combine it with #234. Everything else (using LazyArrays, implementing StaticArray operators) can be done later.

The remaining open question to me is this: should SparseOperator and DenseOperator be separate types or not? The main reason to keep them separate (imo) is if SparseOperator should have fields such as .nzval, i.e. really be a sparse matrix type. Other than that, I am not sure if there is a big advantage to it and one can always access such fields via the .data field. Of course, the current implementation has them separate so the changes might be easier if we keep it that way.

I know I've let this issue go quite stale, but I thought about it a lot. Finally, I decided that subtyping operators to AbstractArray would lead to quite some unwanted behaviour (such as multiplication with arbitrary vectors), and it would be awkward to "forbid" these methods. So I took an alternative approach and changed the types such that one can use any type that implements the AbstractArray interface as data in an operator, see #265. You can now already use LazyArrays, and at some point in the future we will probably replace the current implementation using that.

qojulia / QuantumOptics.jl

Subtype AbstractArray? #236