Indicate QuArray as transposed/dual

[acroy]

We should have some way to indicate that a QuVector or a QuMatrix is transposed/dual. I am thinking of something like JuliaLang/julia#6837, i.e. introducing a dual/transposed type. This would make it easier to define inner products and such.

references: JuliaLang/julia#6837, JuliaLang/julia#4774

[jiahao]

+1

It is safe to assume that QuArrays represent objects in Hilbert space(s), so the dual/transposed type makes a lot of sense. The question in base Julia is whether we want to assume so much structure for the vector space defined by the base Array type.

[i2000s]

@jiahao I would promote the idea of implementing dual vector space in the base Julia. To me, Julia makes vectors as pure one-dimensional objects that is really a good thing. Yet, I feel the vectors in Julia are lack of some native mathematical structures, which might have lead to the problems of vector transposes and so forth. There have been various solutions proposed by @StefanKarpinski, @andreasnoack, @Jutho and others. I think maybe we can use @andreasnoack's idea to implement something for QuArray as a test bed for the vector and matrix algebra in Julia. As @Jutho's work moving on, we shall also see how the tensor algebra integrate with quantum packages. I think those features are useful and crucial for the bra-ket linear algebra in quantum mechanics. We can test those features here immediately while the QuBase.jl is being filled up. Then we would have more confidence to decide if we should give the base Julia a shot.

[Jutho]

I don't think the basic Vector's in Julia need a mathematical structure, they are programming structures, i.e. lists, the computer theory type of vector. Sometimes they represent mathematical vectors, sometimes something different, and sometimes it is more natural to represent mathematical vectors in a different way that doesn't fit in the AbstractVector type hierarchy. If one tries to force a specific mathematical structure in Julia base, it will never be possible to satisfy everybody's need. That being said, there does need to be a simple wrapper type for representing the transpose of vector or matrix, in order to be able to support writing compact matrix algebra expressions involving products of vectors and matrices and their transposes. I think that's also the conclusion from #6837.

[jrevels]

If we use a wrapper type for transposition, how do we handle other wrapper types with QuArrays as fields? This is relevant because DiracArrays in QuDirac are wrappers around QuArrays (I want to avoid having to duplicate functionality that already exists in QuBase).

For example, we can define the following right now:

type QuWrapper{B,T,N,A} <: AbstractQuArray{B,T,N}
    qa::QuArray{B,T,N,A}
end

...and delegate the transpose to qa (e.g. transpose(qw::QuWrapper) = QuWrapper(transpose(qa))). This would not be so simple if transposition wrapped the QuArray in another type, since qa's type annotation in the above definition does not encompass the transpose wrapper type.

This could be a non-issue if we defined things like so:

type QuTranspose{B,T,N,Conj,Q<:AbstractQuArray{B,T,N}} <: AbstractQuArray{B,T,N}
    data::Q
end

type QuWrapper{B,T,N,Q<:AbstractQuArray{B,T,N}} <: AbstractQuArray{B,T,N}
    qa::Q
end

This allows qa to be an instance of QuTranspose, but it's messier than I'd like, mostly since the constraints on the behavior of QuArrays are more well-defined than those on AbstractQuArrays.

[jrevels]

The other potential way forward to solve the above issue is to try to write QuArray methods more generically than they have been written so far. We could define basic accessor methods on QuArray, and then make all of our other methods act on AbstractQuArray. Or even further, we could go the duck typing route for such methods, only using type annotations to avoid ambiguity (which, I believe, is the "Julian" way anyway). I think lessening the type restrictions on the functions we have could put us in a better position to tackle future work.

I'll try to get a PR in this weekend so that we have a more practical basis for discussing this.

[jrevels]

After spending a good amount of time playing around with QuCoeffs and reviewing the code from JuliaLang/julia#6837, I've gone back and implemented lazy Transpose/CTranspose/Conjugate wrapper types for QuArrays and tossed away the QuCoeffs work. In an attempt to kick the tires a bit, the new stuff is currently on master (see 87bb3a652fb5f13fd47ad24ee80bbf11907b1032). It's mostly inspired by JuliaLang/julia@1919d591c892fddae08e2ed9467986954a9f95fe. However, without the limitation on having these types be subtyped to AbstractArray, our implementation becomes a little simpler. I've also included a Conjugate type for the same reasons of view vs. copy consistency that have been previously discussed.

At this point, I've tried a bunch of different ways to wrap eagerly conjugated/transposed coefficients in a way that solves the relevant ambiguity issues, but every method I thought of ended up being too inconsistent to be sustainable. My goal was to solve the problem in such a way that we could fall back to the "multiplication zoo" (A_mul_b etc.) on the wrapped coefficients for most operations, the benefit being that a lot of coefficient container types are going to have those methods already overloaded to act efficiently with that type. When I went to define some of these methods for general QuArray operations, though, enforcing type consistency became nightmarish - I couldn't trust operations on eagerly transposed arrays to return arrays of the expected dimensionality...which is, I suppose, the crux of the issue in the first place.

Thus, I've implemented the aforementioned lazy types and have landed on the side of us doing the work to define efficient multiplication methods for special array representations. Though it'll be more work, it will at the very least be consistent and easier to reason about. Thoughts?

[acroy]

Maybe I am making this too simple, but given the implementation in master, couldn't we just have

*{B,T1,T2,N,A1,A2}(qa::QuArray{B,T1,N,A1}, qb::QuArray{B,T2,N,A2}) = QuArray(coeffs(qa)*coeffs(qb), bases(qa))

*{B,T1,T2,N,A1,A2}(qa::QuArray{B,T1,N,A1}, qb::CTranspose{B,T2,N,A2}) = QuArray(A_mul_Bc(coeffs(qa),coeffs(qb)), bases(qa))

and so on (probably plus some checks and taking care of products of basis functions)? Or maybe this is what you meant? I think this would be fine. Of course we want to have mutating versions as well, which might be trickier ...

I noticed that coeffs(dw::DualWrapper) returns the coefficients of the wrapped QuArray. This might be confusing since getindex etc swap indices. Maybe we could have rawcoeffs or something or just use coeffs(data(dw)) to get the original coefficients and have coeffs return the transposed/conjugated data?

[jrevels]

What I was having trouble with was ensuring the proper return type for an operation by falling back to the A_mul_B family. For example, what should the return type of *{B}(tvec::TranVector{B}, qmat::QuMatrix{B}) be?

Mathematically, it's a row vector - currently, a QuArray row vector is explicitly a TranVector/CTranVector. At_mul_B will return an M<:AbstractMatrix, which we'd have to turn into a V<:AbstractVector, then wrap that in a QuArray and take the transpose.

We can't wrap it in a QuArray without vectorizing the result first, since coefficient dimension is enforced to equal the number of bases present, such that a QuVector really stores a V<:AbstractVector and not a M<:AbstractMatrix (which could actually be a Julian row vector). Previously, this was not the case - the dimension of a QuArray was solely defined by the number of bases present, so that one could have a QuArray row vector of type QuArray{B,T,1,M<:AbstractMatrix}. The change to tying coefficient dimension to number of bases was made in order to solidify the definition of a TranVector as a row vector; otherwise, it would not be distinguishable via type whether a TranVector is a row or column vector.

Now that I've actually written that out, perhaps it wouldn't be a bad thing to go back to conflating row vectors and column vectors. It's still better than conflating row vectors and matrices, and might allow us to fall back to A_mul_B functions without a hassle. I'll try it out and submit a PR if it works out.

Also, the switch from coeffs to rawcoeffs makes a lot of sense, I committed that change here: f5f09080880fd4838d4ab095791388db5e981df5

[acroy]

Since we don't store the actual coefficients in CTranspose we will always have to undo the (c)transpose if we want to wrap anything as a CTranspose, right? For the row vector*matrix=r*M we could of course also calculate (M'*r), wrap the result in a QuVector and wrap that as CTranspose ...

[acroy]

So, basically this should work and is type-safe

*{B,T1,T2,A1,A2}(qa::CTranspose{B,T1,1,A1}, qb::QuArray{B,T2,2,A2}) = CTranspose(QuArray(Ac_mul_B(coeffs(qb),coeffs(qa)), bases(data(qa))))

(I am still using coeffs instead of rawcoeffs). This can be generalized for larger N as well, but assumes that Ac_mul_B is efficient.

[jrevels]

^Perfect, made a PR using this method. Hopefully if there aren't any problems with it we can finally close this issue.

[jrevels]

Closed by #14