--- a/LazyTensors/src/lazy_operations.jl	Tue Jun 25 16:37:35 2019 +0200
+++ b/LazyTensors/src/lazy_operations.jl	Tue Jun 25 17:26:39 2019 +0200
@@ -60,6 +60,7 @@
         return new{T,D,Op,T1,T2}(a,b)
     end
 end
+# TODO: Move Op to be the first parameter? Compare to Binary operations
 
 Base.size(v::LazyElementwiseOperation) = size(v.a)
 
@@ -143,6 +144,26 @@
 
 
 
+
+struct LazyTensorMappingBinaryOperation{Op,T,R,D,T1<:TensorMapping{T,R,D},T2<:TensorMapping{T,R,D}} <: TensorMapping{T,D,R}
+    A::T1
+    B::T2
+
+    @inline function LazyTensorMappingBinaryOperation{Op,T,R,D}(A::T1,B::T2) where {Op,T,R,D, T1<:TensorMapping{T,R,D},T2<:TensorMapping{T,R,D}}
+        return new{Op,T,R,D,T1,T2}(A,B)
+    end
+end
+
+apply(mb::LazyTensorMappingBinaryOperation{:+,T,R,D}, v::AbstractArray{T,D}, I::Vararg) where {T,R,D} = apply(mb.A, v, I...) + apply(mb.B,v,I...)
+apply(mb::LazyTensorMappingBinaryOperation{:-,T,R,D}, v::AbstractArray{T,D}, I::Vararg) where {T,R,D} = apply(mb.A, v, I...) - apply(mb.B,v,I...)
+
+range_size(mp::LazyTensorMappingBinaryOperation{Op,T,R,D}, domain_size::NTuple{D,Integer}) where {Op,T,R,D} = range_size(mp.A, domain_size)
+domain_size(mp::LazyTensorMappingBinaryOperation{Op,T,R,D}, range_size::NTuple{R,Integer}) where {Op,T,R,D} = domain_size(mp.A, range_size)
+
+Base.:+(A::TensorMapping{T,R,D}, B::TensorMapping{T,R,D}) where {T,R,D} = LazyTensorMappingBinaryOperation{:+,T,R,D}(A,B)
+Base.:-(A::TensorMapping{T,R,D}, B::TensorMapping{T,R,D}) where {T,R,D} = LazyTensorMappingBinaryOperation{:-,T,R,D}(A,B)
+
+
 # TODO: Write tests and documentation for LazyTensorMappingComposition
 # struct LazyTensorMappingComposition{T,R,K,D} <: TensorMapping{T,R,D}
 #     t1::TensorMapping{T,R,K}

--- a/LazyTensors/test/runtests.jl	Tue Jun 25 16:37:35 2019 +0200
+++ b/LazyTensors/test/runtests.jl	Tue Jun 25 17:26:39 2019 +0200
@@ -57,6 +57,32 @@
     @test_broken BoundsError == (m*m*v)[7]
 end
 
+@testset "TensorMapping binary operations" begin
+    struct ScalarMapping{T,R,D} <: TensorMapping{T,R,D}
+        λ::T
+    end
+
+    LazyTensors.apply(m::ScalarMapping{T,R,D}, v, i) where {T,R,D} = m.λ*v[i]
+    LazyTensors.range_size(m::ScalarMapping, domain_size) = domain_size
+    LazyTensors.domain_size(m::ScalarMapping, range_sizes) = range_sizes
+
+    A = ScalarMapping{Float64,1,1}(2.0)
+    B = ScalarMapping{Float64,1,1}(3.0)
+
+    v = [1.1,1.2,1.3]
+
+    for i ∈ eachindex(v)
+        @test ((A+B)*v)[i] == 2*v[i] + 3*v[i]
+    end
+
+    for i ∈ eachindex(v)
+        @test ((A-B)*v)[i] == 2*v[i] - 3*v[i]
+    end
+
+    @test range_size(A+B, (3,)) == range_size(A, (3,)) == range_size(B,(3,))
+    @test domain_size(A+B, (3,)) == domain_size(A, (3,)) == domain_size(B,(3,))
+end
+
 @testset "LazyArray" begin
     struct DummyArray{T,D, T1<:AbstractArray{T,D}} <: LazyArray{T,D}
         data::T1