Request For Comment: Arithemtic between Operators and LazyOperators #86
Conversation
Extend some of the current operation with LazyOperators to allow for smoother interface with e.g. LazyTensor. As an example I defined +(a::LazyTensor,LazyTensor) so that it returns LazySum. Note that tensor(a::Operator,b::LazyTensor) only behaves well if Operator does not contain a CompositeBasis. This could probably be fixed by specifying the type more explicitly and returning error if it does contain CompositeBasis. Other solutions are welcome.
|
For context, Matias is making this pull request because he is developing an extension to QuantumOptics that needs it on https://github.com/mabuni1998/WaveguideQED.jl/ |
|
Matias, let's include this in the library itself so that we can see its effect on the test runs. Could you rename it to something that fits the naming conventions of the rest of the files (lower snakecase) and include it here (under the other lazy ops) https://github.com/qojulia/QuantumOpticsBase.jl/blob/master/src/QuantumOpticsBase.jl#L68 After the static analysis from JET and Aqua is done complaining we can start looking into also adding tests. |
Codecov Report
@@ Coverage Diff @@
## master #86 +/- ##
==========================================
+ Coverage 96.45% 96.49% +0.03%
==========================================
Files 26 26
Lines 3360 3392 +32
==========================================
+ Hits 3241 3273 +32
Misses 119 119
... and 1 file with indirect coverage changes 📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more |
src/LazyOperatorAritmetic.jl
Outdated
| function QuantumOpticsBase.:⊗(a::Operator,b::LazySum) | ||
| btotal_l = a.basis_l ⊗ b.basis_l | ||
| btotal_r = a.basis_r ⊗ b.basis_r | ||
| ops = Vector{AbstractOperator}(undef,length(b.operators)) |
There was a problem hiding this comment.
How feasible is it to (1) use a tuple instead of a Vector and (2) use a type more specific than AbstractOperator?
There was a problem hiding this comment.
Although not perfect, we now use a Tuple:
There might be better solutions, but in this we at least use a Tuple.
|
The changes are now included in the library under the file: operators_lazyaritmetic.jl |
|
Tests are failing because +(op1::LazyProduct,op2::LazyProduct) no longer throws and error and thus the following line in the test fails: I will work on changing this as well as adding a couple of tests ensuring that the arithmetic is working as it should. |
|
We will probably have a few rounds of test fixing. Looking at the first test that is failing, it seems you are currently a bit permissive with what operators can be summed up. See this now-failing test https://github.com/qojulia/QuantumOpticsBase.jl/blob/master/test/test_operators_lazyproduct.jl#L66 : in it op1 and op2 have different bases so probably this should raise an error. It would probably be easiest to run the tests locally a couple of times, as the CI tests stop at the first failure, while the default local settings are to run all tests and show you all failures. |
Added test for addition and substraction of LazyTensor and LazyProduct. Also tensor of LazyTensor, LazySum, and LazyProduct.
|
I have now added test for the addition of LazyTensor and LazyProduct. The change for LazyProduct thus now reads: QuantumOpticsBase.jl/test/test_operators_lazyproduct.jl Lines 65 to 79 in 690fff2 Where we now allow for addition of LazyProduct and the ArgumentError is therefore not thrown. Also, why is some tests repeated in test_abstractdata.jl ? I have here just outcommented the following line again for the same reason as above: |
src/operators_lazyproduct.jl
Outdated
| @@ -74,6 +94,16 @@ permutesystems(op::LazyProduct, perm::Vector{Int}) = LazyProduct(([permutesystem | |||
| identityoperator(::Type{LazyProduct}, ::Type{S}, b1::Basis, b2::Basis) where S<:Number = LazyProduct(identityoperator(S, b1, b2)) | |||
|
|
|||
|
|
|||
| #Assume same basis | |||
There was a problem hiding this comment.
What does "assume same basis" means here? Is it something that would lead to bugs if not followed? Can we have tests for it, even if the test is basically just checking that an error is thrown instead of silently continuing with the operation?
There was a problem hiding this comment.
What I meant is that BL and BR are the same since I'm not sure the methods defined would work if one considered operators with BL != BR. This is already handled by the dispatch, and the user should just get an error saying that addition for these types of operators is not defined. I will look into this and create specific tests for it.
There was a problem hiding this comment.
BL == BR is no longer assumed and code is updated accordingly
test/test_abstractdata.jl
Outdated
| # Addition | ||
| @test_throws ArgumentError op1 + op2 | ||
| @test_throws ArgumentError op1 - op2 | ||
| # Addition Addition of LazyTensor now returns LazySum |
There was a problem hiding this comment.
Before merge I would suggest deleting the commented-out tests. They would not be informative to future readers, long after this merge. Future readers should use git blame if they want to know why a change was made.
Could you add tests for the failure modes of this contribution? E.g., @test_throws for where your code specifically raises an error to warn about something it can not do. And are there any other implicit assumptions about shapes and bases and other things? We do not want the user to obtain nonsense results if they do not happen to know that some combination of operations is not supported.
There was a problem hiding this comment.
Updated code to no longer assume anything about BL == BR and added test to ensure that errors are thrown if addition between operators is not allowed.
src/QuantumOpticsBase.jl
Outdated
| @@ -80,5 +80,6 @@ include("metrics.jl") | |||
| include("spinors.jl") | |||
| include("phasespace.jl") | |||
| include("printing.jl") | |||
| #include("operators_lazyaritmetic.jl") | |||
There was a problem hiding this comment.
this probably should be deleted before merging
There was a problem hiding this comment.
Oh yeah... Done that now ;)
|
Thanks for this @mabuni1998. I am definitely in favor of making laziness "contagious" in the way you have done here. What do you think @david-pl? |
| function +(a::LazyTensor{B1,B2}, b::LazyTensor{B1,B2}) where {B1,B2} | ||
| if length(a.indices) == 1 && a.indices == b.indices | ||
| op = a.operators[1] * a.factor + b.operators[1] * b.factor | ||
| return LazyTensor(a.basis_l, a.basis_r, a.indices, (op,)) |
There was a problem hiding this comment.
This is the only place I am not sure about defaulting to laziness. It's quite a special case, but I encounter it quite a bit. I suppose the reason to do lazy summing here is mainly to be consistent with the laziness-preserving principle. I have some code that makes use of the existing behavior, but of course I can still do this kind of concrete summing manually if I want to, so I'm not arguing hard to keep it. What are your thoughts?
There was a problem hiding this comment.
My experience was that the previous implementation was very limiting. Especially since the custom operators I have been playing around with were not DataOperators but AbstractOperators, where the operation was defined via a function rather than a matrix. Therefore, these cannot be trivially added (except by using LazySum), and the above implementation fails. Also, length(a.indices) ==1 is required, and I could imagine situations where one would like to be able to add LazyTensors containing more than one operator.
However, one could perhaps keep the original behavior by dispatching on LazyTensors containing only one DataOperator. That is adding a function like this (draft, I'm not entirely sure it works):
const single_dataoperator{B1,B2} = LazyTensor{B1,B2,ComplexF64,Vector{Int64},Tuple{T}} where {B1,B2,T<:DataOperator}
function +(a::T1,b::T2) where {T1 <: single_dataoperator{B1,B2},T2 <: single_dataoperator{B1,B2}}
if length(a.indices) == 1 && a.indices == b.indices
op = a.operators[1] * a.factor + b.operators[1] * b.factor
return LazyTensor(a.basis_l, a.basis_r, a.indices, (op,))
end
throw(ArgumentError("Addition of LazyTensor operators is only defined in case both operators act nontrivially on the same, single tensor factor."))
end
There was a problem hiding this comment.
@mabuni1998 I think it's worth trying to keep the original intact as you suggest. If we can handle it via dispatch, we won't lose anything. Or am I missing some case here?
There was a problem hiding this comment.
No I don't think we will lose anything. I have implemented to above as:
There was a problem hiding this comment.
Thanks for finding a way to keep the original behavior. This is not type-stable, but I can't think of an obvious way to make it otherwise, except by letting LazyTensor indices be type parameters. Maybe it doesn't matter too much, as this will not typically be performance-critical?
There was a problem hiding this comment.
Probably won't be performance-critical no, as you are most likely creating the operators once at the beginning of the simulation and then not changing them as you do multiplications etc.
Abstract type LazyOperator added as supertype for LazySum, LazyTensor, and LazyProduct. +(a::LazyOperator,b::LazyOperator) implemented to make LazyOperators contagious. added identityoperator for composite basis since previously the following would return false: ``` b1a = GenericBasis(2) b1b = GenericBasis(3) b2a = GenericBasis(1) b2b = GenericBasis(4) b3a = GenericBasis(6) b3b = GenericBasis(5) I1 = dense(identityoperator(b1a, b1b)) I2 = dense(identityoperator(b2a, b2b)) I3 = dense(identityoperator(b3a, b3b)) a = tensor(I1,I3) a == identityoperator(a) ```
|
I don't know if this was intended but previously I found that: b1a = GenericBasis(2)
b1b = GenericBasis(3)
b2a = GenericBasis(1)
b2b = GenericBasis(4)
b3a = GenericBasis(6)
b3b = GenericBasis(5)
I1 = dense(identityoperator(b1a, b1b))
I2 = dense(identityoperator(b2a, b2b))
I3 = dense(identityoperator(b3a, b3b))
a = tensor(I1,I3)
a == identityoperator(a)would result in false. To fix it, I implemented identityoperator() when having CompositeBasis. See the following for the implementation: This was important because in LazyTensor we check whether the operator is identityoperator and when the above didn't return true it lead to unwanted errors. |
src/operators_dense.jl
Outdated
| @@ -234,6 +234,9 @@ permutesystems(a::AdjointOperator{B1,B2}, perm) where {B1<:CompositeBasis,B2<:Co | |||
|
|
|||
| identityoperator(::Type{S}, ::Type{T}, b1::Basis, b2::Basis) where {S<:DenseOpType,T<:Number} = | |||
| Operator(b1, b2, Matrix{T}(I, length(b1), length(b2))) | |||
| identityoperator(::Type{S}, ::Type{T}, b1::CompositeBasis, b2::CompositeBasis) where {S<:DenseOpType,T<:Number} = | |||
| Operator(b1, b2, kron([Matrix{T}(I, length(b1.bases[i]), length(b2.bases[i])) for i in reverse(1:length(b1.bases))]...)) | |||
|
|
|||
There was a problem hiding this comment.
This is strange... I am probably missing something obvious, but it looks like these two variants should give exactly equal results for CompositeBasis?
There was a problem hiding this comment.
Yes, I agree they should, but from my testing they don't. Try running the following example:
edit: I should specify that with the above addition the below example is obviously true, while in the current release version of QuantumOpticsBase, the below example won't return true.
I don't know if this was intended but previously I found that:
b1a = GenericBasis(2) b1b = GenericBasis(3) b2a = GenericBasis(1) b2b = GenericBasis(4) b3a = GenericBasis(6) b3b = GenericBasis(5) I1 = dense(identityoperator(b1a, b1b)) I2 = dense(identityoperator(b2a, b2b)) I3 = dense(identityoperator(b3a, b3b)) a = tensor(I1,I3) a == identityoperator(a)would result in false. To fix it, I implemented identityoperator() when having CompositeBasis. See the following for the implementation:
This was important because in LazyTensor we check whether the operator is identityoperator and when the above didn't return true it lead to unwanted errors.
There was a problem hiding this comment.
That's odd. Could it be that the eltype is different? Maybe that identityoperator() method is broken?
There was a problem hiding this comment.
It's the data inside; they do not represent the same matrix. Maybe Julias UnitScaling follows another convention than we expect, or maybe the example I'm making shouldn't return true?
To make it even clearer, consider:
b1a = GenericBasis(2)
b1b = GenericBasis(3)
I1 = identityoperator(b1a, b1b)
I2 = identityoperator(b1b, b1a)
a = tensor(I1,I2)
julia> a.data
6×6 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 4 stored entries:
1.0+0.0im ⋅ ⋅ ⋅ ⋅ ⋅
⋅ 1.0+0.0im ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 1.0+0.0im ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 1.0+0.0im ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅
julia> identityoperator(a).data
6×6 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 6 stored entries:
1.0+0.0im ⋅ ⋅ ⋅ ⋅ ⋅
⋅ 1.0+0.0im ⋅ ⋅ ⋅ ⋅
⋅ ⋅ 1.0+0.0im ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 1.0+0.0im ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 1.0+0.0im ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 1.0+0.0imThere was a problem hiding this comment.
Oh... Yeah, I think in this case they actually should not be equal. It's kind of funny to even define identityoperator() in the rectangular case. Clearly it's not the identity!
There was a problem hiding this comment.
To summarize, the two methods of identityoperator() above are only equivalent in case the left and right composite bases both have the same length and have the same number of factors, and the dimensions of the factors match. Otherwise the first definition (the original one) is correct - it actually produces an identity operator in the square case - where the second is not. The second version also fails in case the left basis has more factors than the right.
Should we consider raising a warning if identityoperator() is used in rectangular cases?
There was a problem hiding this comment.
I see good point. But then I would like some input on how we should treat LazyTensor with non-equivalent left and right basises. See below
There was a problem hiding this comment.
Okay, so that method here should be removed again.
src/operators_lazytensor.jl
Outdated
| function tensor(a::LazyTensor{B1,B2},b::AbstractOperator{B3,B4}) where {B1,B2,B3,B4} | ||
| if isequal(b,identityoperator(b)) | ||
| btotal_l = a.basis_l ⊗ b.basis_l | ||
| btotal_r = a.basis_r ⊗ b.basis_r | ||
| LazyTensor(btotal_l,btotal_r,a.indices,(a.operators...,),a.factor) | ||
| elseif B3 <: CompositeBasis || B4 <: CompositeBasis | ||
| throw(ArgumentError("tensor(a::LazyTensor{B1,B2},b::AbstractOperator{B3,B4}) is not implemented for B3 or B4 being CompositeBasis unless b is identityoperator ")) | ||
| else | ||
| a ⊗ LazyTensor(b.basis_l,b.basis_r,[1],(b,),1) | ||
| end |
There was a problem hiding this comment.
How should we treat LazyTensor with non-equivalent LR and BR. The first if statement is only well defined for BL==BR it seems.
There was a problem hiding this comment.
Yeah, so LazyTensor with non-rectangular factors in the composite basis implicitly defines simple isometries, like you constructed in your new implementation of identityoperator(). Honestly I think it's okay to just not handle the b is an identity case specially at all and just include the explicit identity/isometry operator. It's kind of unfortunate to explicitly check that a matrix is the identity here anyway if I'm expecting a lazy operation.
I have an embed_lazy() that I use to achieve this kind of thing.
Perhaps a nicer solution is to introduce a proper LazyIdentity/LazyIsometry operator that gets materialized only when necessary. Then we could implement tensor(a::LazyTensor, b::LazyIdentity) using something like your first conditional branch here.
There was a problem hiding this comment.
@amilsted , could you put up a gist with embed_lazy or just an incomplete pull request?
This discussion and embed_lazy would probably have a bearing on this issue as well: qojulia/QuantumOptics.jl#322
There was a problem hiding this comment.
I'd also say we can ignore the special isequal(b,identityoperator(b)). It should (hopefully) be the exception to run into this, so the impact shouldn't be too large.
One alternative way of solving this would be doing something like #90, but that's just an idea. For now I'd say don't special handle it, remove the new identityoperator method above and we're good to go.
By the way, thanks for all the work @mabuni1998 and thanks for the reviews @amilsted @Krastanov. I appreciate it!
There was a problem hiding this comment.
You're welcome! I now pushed a new version with the above things fixed.
I would like to say that originally I put the check for identityoperator there because, in my own code, I noticed performance drops when doing stuff like the following (WaveguideBasis is a custom basis with operators defined as functions, not as matrices, and they, therefore, have to be Lazy):
bc = FockBasis(2)
bw = WaveguideBasis(2,times)
adw = tensor(create(bc),identityoperator(bc),identityoperator(bc),destroy(bw))
Because here, identityoperator(bc) is explicitly stored in the resulting LazyTensor (there is a custom method creating LazyTensor here). Probably this is better solved by allowing operators spanning multiple basises to be stored in LazyTensor. Right now, we explicitly check for this case and don't allow it. However, this would require extra work (and I don't know how much) in the mul! method so let's leave it. This problem (I think) would not be solved by implementing LazyIdentity since identityoperator(bc) would not return such abjoct.
Also, I'm aware that in the above, one could just use embed or explicitly define LazyTensor. I simply just wanted to make it easier for the user to have better performance.
Furthermore, my use case is probably niece, since I'm only using Lazy operators because my operators are defined as functions, whereas in most cases, you use it because evaluating the tensor product is expensive. I can just define custom methods that take care of this.
Also removed identityopreator method for CompositeBasis
|
@amilsted @Krastanov @david-pl Are we ready to merge this? |
I think I'd like to do one more patch release first. |
|
@mabuni1998 @amilsted It's good to merge from my side. @amilsted please merge (or let me know) once you're done with the patch release. |
Extend some of the current operation with LazyOperators to allow for smoother interface with e.g. LazyTensor. As an example I defined +(a::LazyTensor,LazyTensor) so that it returns LazySum.
Note that tensor(a::Operator,b::LazyTensor) only behaves well if Operator does not contain a CompositeBasis. This could probably be fixed by specifying the type more explicitly and returning error if it does contain CompositeBasis. Other solutions are welcome.