diff --git a/src/operators/localoperator.jl b/src/operators/localoperator.jl
index ffe9e34a..65796b88 100644
--- a/src/operators/localoperator.jl
+++ b/src/operators/localoperator.jl
@@ -1,20 +1,62 @@
-
 # Hamiltonian consisting of local terms
 # -------------------------------------
+"""
+    struct LocalOperator{T<:Tuple,S}
+
+A sum of local operators acting on a lattice. The lattice is stored as a matrix of vector spaces,
+and the terms are stored as a tuple of pairs of indices and operators.
+
+# Fields
+
+- `lattice::Matrix{S}`: The lattice on which the operator acts.
+- `terms::T`: The terms of the operator, stored as a tuple of pairs of indices and operators.
+
+# Constructors
+
+    LocalOperator(lattice::Matrix{S}, terms::Pair...)
+    LocalOperator{T,S}(lattice::Matrix{S}, terms::T) where {T,S} # expert mode
+
+# Examples
+
+```julia
+lattice = fill(ℂ^2, 1, 1) # single-site unitcell
+O1 = LocalOperator(lattice, ((1, 1),) => σx, ((1, 1), (1, 2)) => σx ⊗ σx, ((1, 1), (2, 1)) => σx ⊗ σx)
+```
+"""
 struct LocalOperator{T<:Tuple,S}
     lattice::Matrix{S}
     terms::T
-end
-function LocalOperator(lattice::Matrix{S}, terms::Pair...) where {S}
-    lattice′ = PeriodicArray(lattice)
-    for (inds, operator) in terms
-        @assert operator isa AbstractTensorMap
-        @assert numout(operator) == numin(operator) == length(inds)
-        for i in 1:length(inds)
-            @assert space(operator, i) == lattice′[inds[i]]
+    function LocalOperator{T,S}(lattice::Matrix{S}, terms::T) where {T,S}
+        plattice = PeriodicArray(lattice)
+        # Check if the indices of the operator are valid with themselves and the lattice
+        for (inds, operator) in terms
+            @assert operator isa AbstractTensorMap
+            @assert numout(operator) == numin(operator) == length(inds)
+            @assert spacetype(operator) == S
+
+            for i in 1:length(inds)
+                @assert space(operator, i) == plattice[inds[i]]
+            end
         end
+        return new{T,S}(lattice, terms)
     end
-    return LocalOperator{typeof(terms),S}(lattice, terms)
+end
+function LocalOperator(
+    lattice::Matrix,
+    terms::Pair...;
+    atol=maximum(x -> eps(real(scalartype(x[2])))^(3 / 4), terms),
+)
+    allinds = getindex.(terms, 1)
+    alloperators = getindex.(terms, 2)
+
+    relevant_terms = []
+    for inds in unique(allinds)
+        operator = sum(alloperators[findall(==(inds), allinds)])
+        norm(operator) > atol && push!(relevant_terms, inds => operator)
+    end
+
+    terms_tuple = Tuple(relevant_terms)
+    return LocalOperator{typeof(terms_tuple),eltype(lattice)}(lattice, terms_tuple)
 end
 
 """
@@ -27,6 +69,9 @@ while the second version throws an error if the lattices do not match.
 function checklattice(args...)
     return checklattice(Bool, args...) || throw(ArgumentError("Lattice mismatch."))
 end
+function checklattice(::Type{Bool}, H1::LocalOperator, H2::LocalOperator)
+    return H1.lattice == H2.lattice
+end
 function checklattice(::Type{Bool}, peps::InfinitePEPS, O::LocalOperator)
     return size(peps) == size(O.lattice)
 end
@@ -57,3 +102,22 @@ function Base.repeat(O::LocalOperator, m::Int, n::Int)
     end
     return LocalOperator(lattice, terms...)
 end
+
+# Linear Algebra
+# --------------
+function Base.:*(α::Number, O::LocalOperator)
+    scaled_terms = map(((inds, operator),) -> (inds => α * operator), O.terms)
+    return LocalOperator{typeof(scaled_terms),eltype(O.lattice)}(O.lattice, scaled_terms)
+end
+Base.:*(O::LocalOperator, α::Number) = α * O
+
+Base.:/(O::LocalOperator, α::Number) = O * inv(α)
+Base.:\(α::Number, O::LocalOperator) = inv(α) * O
+
+function Base.:+(O1::LocalOperator, O2::LocalOperator)
+    checklattice(O1, O2)
+    return LocalOperator(O1.lattice, O1.terms..., O2.terms...)
+end
+
+Base.:-(O::LocalOperator) = -1 * O
+Base.:-(O1::LocalOperator, O2::LocalOperator) = O1 + (-O2)