WIP: Add new types QuStateVec and QuDensityMatrix.

[acroy]

This adds two new types QuStateVec (for state vectors/wave functions) and QuDensityMatrix (for density matrices). I am not sure we want/need the former, but since I had implemented it anyways I include it here.

Two new functions, populations and purity, are defined, which return the occupation probability and the purity, respectively. I also added trace and sqrtm for AbstractQuMatrix.

Tests are missing so far.

[jrevels]

I'm not sure what QuStateVec gives us that we don't already have with QuVector. Likewise with this implementation of QuDensityMatrix and QuMatrix. Why not just define the functions you provided here on directly on AbstractQuVector and AbstractQuMatrix?

My assumption was that we were going to have a QuDensityMatrix as a lazy construction, with the ability to convert to a full matrix if desired. In fact, I vote we have a lazy construction for outer products in general:

immutable QuOuter{B,KT<:AbstractQuVector, BR<:DualVector} <: AbstractQuMatrix{B}
   ket::KT
   bra::BR
   # this relies on a generalized CTranspose type
   QuDensity(ket::AbstractQuVector{B}, bra::DualVector{B}) = new(ket, bra)
end

And then define getindex as something like:

getindex(qo::QuOuter, i, j) = qo.ket[i] * qo.bra[j]

And then we define things like size and trace accordingly. This would allow for quick lazy construction of outer products and density matrices for common cases like the partial trace, where most users will want to enter something like

julia> ptrace(ψ*ψ', 1) # trace out the first subsystem

e.g., a calculation that only takes the outer product once and does not require repeated access to the resulting elements.

We can overload full to take the eager outer product, and then define things like

setindex!(::QuOuter, i...) = error("QuOuter is an immutable representation of an outer product of states. Mutating functions can only be called on mutable representations, which can be obtained using full(::QuOuter)."

Anyway, I was gonna tackle #31 today, then move on to making a PR that defines the above stuff, but it seems like that PR would be competing with this one. Thoughts on the above?

[acroy]

Yeah, I am not sure about QuStateVec, but I think some sort of density matrix type is useful/necessary. Firstly, note that your QuOuter only covers pure states. A general (mixed) state is a sum of outer products (which I also would like to support in some way; this connects to the density matrix decompositions I mentioned in #7). Secondly, the norm of a density matrix is (unfortunately) defined by the trace not the matrix dot product (although that is probably only a minor point). The question is if we want to have a distinction between states and operators? Purity, populations, fidelity and so on only make sense for states and typically one treats states and operators differently (for example, the time evolution is different), so it might make sense to be able to dispatch on states?

Independent of that, I think your suggestions regarding a lazy outer product sound like a great idea and there is probably no conflict with this PR.

[amitjamadagni]

I was thinking whether we could have something like this

type QuDensityMatrix
            prob::Array{Float64,1}
            state::QuStateVec
            QuDensityMatrix(prob::Array{Float64,1}, state::QuStateVec) = new(prob, state)
end

QuDensityMatrix(qa::QuDensityMatrix) = QuMatrix(#summation_of_states_outerproduct given by coeffs(qa.prob * (qa.state*(qa.StateVec)'))

I guess using this we could even fix conversion of state to density matrix using this. Please do let me know if I am missing something on this.

[jrevels]

Firstly, note that your QuOuter only covers pure states. A general (mixed) state is a sum of outer products (which I also would like to support in some way; this connects to the density matrix decompositions I mentioned in #7).

What I'm saying is that we already have a type which explicitly covers this behavior: AbstractQuMatrix. So, following the example I gave above, ptrace(ψ*ψ', 1) would return a QuMatrix. Any decompositions you might want to define on QuDensityMatrix are also definable on AbstractQuMatrix. I also think, API-wise, it would be strange not to define them this way - everything you'd want to define on QuDensityMatrix should be defined on AbstractQuMatrix as well, since not doing so would likely be confusing to the end-user (e.g. "I constructed this operator that I know is a density operator, and acts exactly like a density operator, but can't be treated like a density operator?"). At that point, the only effect of having a QuDensityMatrix type at all is to have the obligation to maintain more conversion rules.

Secondly, the norm of a density matrix is (unfortunately) defined by the trace not the matrix dot product (although that is probably only a minor point).

...and the trace over a QuOuter can be defined as:

function trace(qo::QuOuter)
   over =  min(size(qo, 1), size(qo, 2))
   return sum(i -> getindex(qo, i, i), 1:over)
end

Purity, populations, fidelity and so on only make sense for states and typically one treats states and operators differently (for example, the time evolution is different), so it might make sense to be able to dispatch on states?

The point, I think, is that any AbstractQuVector can act like and be considered a state (with CTranspose indicating it's dual). Any AbstractQuMatrix can act like and can be considered an operator (with CTranspose indicating it's adjoint). I can't think of any reason to implement types that, mechanically, are the exact same as the type we already have, which covers all of these cases and can be dispatched on just as easily (using AbstractQuVector and AbstractQuMatrix).

[jrevels]

@amitjamadagni If I understand the type you gave (which I might not, correct me if I'm wrong), it is equivalent to an outer product of pure states - states that can be straightforwardly represented as a superposition of other states (e.g., a QuVector), are always pure. Thus, it couldn't be used to represent a mixed state (e.g. a reduced density matrix), which is more general and lacks the cartesian product structure you get with a simple outer product of states.

[jrevels]

@acroy I may have misinterpreted the intent of your statement about the norm being defined by the trace - so you're saying we'd have conflicting definitions for norm for a QuMatrix if we use QuMatrix for density operators?

[jrevels]

@acroy Also, an addendum to this statement:

Any AbstractQuMatrix can act like and can be considered an operator (with CTranspose indicating it's adjoint).

The way I've been thinking of it is that a density operator does represent a state, but what it is is an operator.

So, re-reading your comment, your claim is that there will be functions for which we want a different definition for QuMatrix than for QuDensityMatrix (e.g. norm)? If that is the case, I think a list of such functions would be really beneficial (even if they don't necessarily get implemented immediately).

[acroy]

@jrevels : I basically agree with what you say and, in particular, I would also like to avoid redundant types. Unfortunately, states (= density matrices) and operators do not always behave in the same way. The first example is norm, like you noticed above for operators (=QuMatrix) one uses trace(qa'*qa) (which is consistent with the inner product for operators), but for density matrices one uses trace(dm). The second example is (unitary) time evolution. If U=expm(-im*H*t) then dm(t)=U*dm(0)*U' for density matrices and o(t) = U'*o(0)*U for operators.

But we could also just use the trace norm for operators as well and accept the inconsistency with the inner product. However, I am still wondering if it would be convenient in some cases to say that an argument has to be a density matrix/state (like for propagation)?

[amitjamadagni]

@jrevels I have been thinking on these lines. Like given the probability and states we have the summation over the product of respective probabilities of the states and the outer products of the states.

[jrevels]

Hmm. I've normally thought of the density matrix as a representation of an operator describing an ensemble of states, rather than a representation of a state in and of itself, but I guess that both concepts are equally valid when referring to mixed states. The separate time evolution behavior makes sense to me after looking at a derivation from the Von Neumann equation. I couldn't find anything on defining the norm to be the square root of the trace, but it makes sense to me if the diagonal elements of the matrix correspond to the probability amplitudes for the basis states.

That brings up the point: Is it feasible that one would actually want to treat a density matrix as an operator or a state, given the context of different calculations? For example, would a user ever want to take the traditional Frobenius norm for the density operator, rather than the norm of the state that operator describes?

If so, it might make more sense to distinguish the behavior on a functional level (e.g. statenorm vs. norm) rather than on the type level. Or, having types that capture each behavior (QuDensityState, QuDensityOp).

Anyway, this PR makes a bit more sense to me now, thanks for your patience.

@amitjamadagni So the type you listed represents the sum (mixing Julia notation and math) ∑ prob[i] * | i >< i |? This could only work in the case where the off-diagonal elements of the matrix are zero, i.e., there are no coherences in the mixed state.

[acroy]

I've normally thought of the density matrix as a representation of an operator describing an ensemble of states, rather than a representation of a state in and of itself

I think this is actually a pretty subtle point, which touches the ontology of density matrices and state vectors. Some consider the latter as fundamental, while others argue for the opposite. Fortunately we don't have to decide this here :-)

Is it feasible that one would actually want to treat a density matrix as an operator or a state, given the context of different calculations?

Good question. One certainly wants the DM to be an operator in most situations, just for some (interpretation related) questions one needs special behavior. The issue with norm is only minor I would say, after all there is vecnorm and norm(A,p) in Base, so statenorm or such is probably reasonable. I think the question is, whether we want to distinguish between the concepts "state" and "operator" and reflect this by different subtypes. I always thought this is useful, but I can see now that this might cause confusion.

Regarding Amit's suggestion: the sum is like you say, but |i> is actually an eigenstate of some density matrix. His proposal corresponds basically to the Eigen type returned by eigfact. This kind of decomposition can be really useful for Monte-Carlo like calculations.

[jrevels]

After being absent from QuBase.jl for a while, I'm fairly convinced that traits are the correct way to do this. As we've already discussed, density matrices can be treated like states or operators depending on context. Thus, a density matrix could be defined as a subtype of an AbstractQuMatrix that supports both a QuState interface and a QuOperator interface.