diff --git a/examples/bivariate.jl b/examples/bivariate.jl
new file mode 100644
index 0000000..4321d31
--- /dev/null
+++ b/examples/bivariate.jl
@@ -0,0 +1,48 @@
+# This file is a part of StarAlgebras.jl. License is MIT: https://github.com/JuliaAlgebra/StarAlgebras.jl/blob/main/LICENSE
+# Copyright (c) 2021-2025: Marek Kaluba, Benoît Legat
+
+# Example implementation of Bivariate polynomials
+# See MultivariatePolynomials.jl for general implementation
+struct ExponentsIterator end
+Base.eltype(::Type{ExponentsIterator}) = NTuple{2,Int}
+Base.IteratorSize(::Type{ExponentsIterator}) = Base.IsInfinite()
+
+function Base.iterate(::ExponentsIterator)
+    z = (0, 0)
+    return z, (z, 0)
+end
+
+function Base.iterate(::ExponentsIterator, state)
+    z, deg = state
+    if iszero(z[2])
+        deg += 1
+        z = (0, deg)
+    else
+        z = (z[1] + 1, z[2] - 1)
+    end
+    return z, (z, deg)
+end
+
+function grlex(a::NTuple{2,Int}, b::NTuple{2,Int})
+    return isless((sum(a), a), (sum(b), b))
+end
+
+struct Monomial
+    exponents::NTuple{2,Int}
+end
+
+Base.one(::Monomial) = Monomial((0, 0))
+Base.:*(a::Monomial, b::Monomial) = Monomial(a.exponents .+ b.exponents)
+
+monomial(exp) = Monomial(exp)
+exponents(mono::Monomial) = mono.exponents
+
+SA.comparable(::ExponentsIterator) = grlex
+
+function bivariate_algebra()
+    exps = ExponentsIterator()
+    basis = SA.MappedBasis(exps, monomial, exponents)
+    mstr = SA.DiracMStructure(basis, *)
+    object = Monomial((0, 0))
+    return SA.StarAlgebra(object, mstr)
+end
diff --git a/src/bases.jl b/src/bases.jl
index aa6e6b9..4e1e544 100644
--- a/src/bases.jl
+++ b/src/bases.jl
@@ -38,11 +38,15 @@ key_type(b::AbstractBasis) = key_type(typeof(b))
 Implicit bases are bases that contains the product of all its elements.
 This makes these bases particularly useful to work with [`AlgebraElement`](@ref)s with supports that can not be reasonably bounded.
 Note that these bases may not explictly store its elements in memory as they may be potentially infinite.
+The elements of the basis are iterated in increasing order according to `comparable(object(b))`.
 """
 abstract type ImplicitBasis{T,I} <: AbstractBasis{T,I} end
 
-function zero_coeffs(::Type{S}, ::ImplicitBasis{T,I}) where {S,T,I}
-    return SparseCoefficients(I[], S[])
+comparable(::Type) = isless
+comparable(object) = comparable(eltype(object))
+
+function zero_coeffs(::Type{S}, basis::ImplicitBasis{T,I}) where {S,T,I}
+    return SparseCoefficients(I[], S[], comparable(object(basis)))
 end
 
 """
diff --git a/src/sparse_coeffs.jl b/src/sparse_coeffs.jl
index e7df360..9750424 100644
--- a/src/sparse_coeffs.jl
+++ b/src/sparse_coeffs.jl
@@ -1,32 +1,33 @@
 # This file is a part of StarAlgebras.jl. License is MIT: https://github.com/JuliaAlgebra/StarAlgebras.jl/blob/main/LICENSE
 # Copyright (c) 2021-2025: Marek Kaluba, Benoît Legat
 
-struct SparseCoefficients{K,V,Vk,Vv} <: AbstractCoefficients{K,V}
+struct SparseCoefficients{K,V,Vk,Vv,L} <: AbstractCoefficients{K,V}
     basis_elements::Vk
     values::Vv
+    isless::L
 end
 
-function SparseCoefficients(elts::Ks, vals::Vs) where {Ks,Vs}
-    return SparseCoefficients{eltype(elts),eltype(vals),Ks,Vs}(elts, vals)
+function SparseCoefficients(elts::Ks, vals::Vs, isless = isless) where {Ks,Vs}
+    return SparseCoefficients{eltype(elts),eltype(vals),Ks,Vs,typeof(isless)}(elts, vals, isless)
 end
 
 Base.keys(sc::SparseCoefficients) = sc.basis_elements
 Base.values(sc::SparseCoefficients) = sc.values
 function Base.copy(sc::SparseCoefficients)
-    return SparseCoefficients(copy(keys(sc)), copy(values(sc)))
+    return SparseCoefficients(copy(keys(sc)), copy(values(sc)), sc.isless)
 end
 
-function _search(keys::Tuple, key)
+function _search(keys::Tuple, key; lt)
     # `searchsortedfirst` is not defined for `Tuple`
     return findfirst(isequal(key), keys)
 end
 
-function _search(keys, key::K) where {K}
-    return searchsortedfirst(keys, key; lt = comparable(K))
+function _search(keys, key::K; lt) where {K}
+    return searchsortedfirst(keys, key; lt)
 end
 
 function Base.getindex(sc::SparseCoefficients{K}, key::K) where {K}
-    k = _search(sc.basis_elements, key)
+    k = _search(sc.basis_elements, key; lt = sc.isless)
     if k in eachindex(sc.basis_elements)
         v = sc.values[k]
         if sc.basis_elements[k] == key
@@ -40,7 +41,7 @@ function Base.getindex(sc::SparseCoefficients{K}, key::K) where {K}
 end
 
 function Base.setindex!(sc::SparseCoefficients{K}, val, key::K) where {K}
-    k = searchsortedfirst(sc.basis_elements, key; lt = comparable(K))
+    k = searchsortedfirst(sc.basis_elements, key; lt = sc.isless)
     if k in eachindex(sc.basis_elements) && sc.basis_elements[k] == key
         sc.values[k] = val
     else
@@ -100,7 +101,7 @@ _first_sparse_coeffs(c::SparseCoefficients, args...) = c
 _first_sparse_coeffs(_, args...) = _first_sparse_coeffs(args...)
 
 function Base.zero(sc::SparseCoefficients)
-    return SparseCoefficients(empty(keys(sc)), empty(values(sc)))
+    return SparseCoefficients(empty(keys(sc)), empty(values(sc)), sc.isless)
 end
 
 _similar(x::Tuple) = _similar(x, typeof(x[1]))
@@ -111,16 +112,16 @@ _similar(x, ::Type{T}) where {T} = similar(x, T)
 _similar_type(::Type{<:Tuple}, ::Type{T}) where {T} = Vector{T}
 _similar_type(::Type{V}, ::Type{T}) where {V,T} = similar_type(V, T)
 
-function similar_type(::Type{SparseCoefficients{K,V,Vk,Vv}}, ::Type{T}) where {K,V,Vk,Vv,T}
-    return SparseCoefficients{K,T,_similar_type(Vk, K),_similar_type(Vv, T)}
+function similar_type(::Type{SparseCoefficients{K,V,Vk,Vv,L}}, ::Type{T}) where {K,V,Vk,Vv,T,L}
+    return SparseCoefficients{K,T,_similar_type(Vk, K),_similar_type(Vv, T),L}
 end
 
 function Base.similar(s::SparseCoefficients, ::Type{T} = value_type(s)) where {T}
-    return SparseCoefficients(collect(s.basis_elements), _similar(s.values, T))
+    return SparseCoefficients(collect(s.basis_elements), _similar(s.values, T), s.isless)
 end
 
 function map_keys(f::Function, s::SparseCoefficients)
-    return SparseCoefficients(map(f, s.basis_elements), s.values)
+    return SparseCoefficients(map(f, s.basis_elements), s.values, s.isless)
 end
 
 function MA.mutability(
@@ -144,29 +145,24 @@ function __prealloc(X::SparseCoefficients, Y::SparseCoefficients, op)
     return similar(X, T)
 end
 
-comparable(::Type) = isless
-function MA.operate!(::typeof(canonical), res::SparseCoefficients)
-    return MA.operate!(canonical, res, comparable(key_type(res)))
-end
-
 function unsafe_push!(res::SparseCoefficients, key, value)
     push!(res.basis_elements, key)
     push!(res.values, value)
     return res
 end
 
-# `::C` is needed to force Julia specialize on the function type
+# `{...,L}` is needed to force Julia specialize on the function type
 # Otherwise, we get one allocation when we call `issorted`
 # See https://docs.julialang.org/en/v1/manual/performance-tips/#Be-aware-of-when-Julia-avoids-specializing
-function MA.operate!(::typeof(canonical), res::SparseCoefficients, cmp::C) where {C}
-    sorted = issorted(res.basis_elements; lt = cmp)
+function MA.operate!(::typeof(canonical), res::SparseCoefficients{K,V,Vk,Vv,L}) where {K,V,Vk,Vv,L}
+    sorted = issorted(res.basis_elements; lt = res.isless)
     distinct = allunique(res.basis_elements)
     if sorted && distinct && !any(iszero, res.values)
         return res
     end
 
     if !sorted
-        p = sortperm(res.basis_elements; lt = cmp)
+        p = sortperm(res.basis_elements; lt = res.isless)
         permute!(res.basis_elements, p)
         permute!(res.values, p)
     end
diff --git a/src/star.jl b/src/star.jl
index 9b81d68..212f56b 100644
--- a/src/star.jl
+++ b/src/star.jl
@@ -16,7 +16,7 @@ star(::AbstractBasis, x) = star(x)
 function star(b::AbstractBasis, d::SparseCoefficients)
     k = star.(Ref(b), keys(d))
     v = star.(values(d))
-    return SparseCoefficients(k, v)
+    return SparseCoefficients(k, v, d.isless)
 end
 
 function star(b::FixedBasis, i::Integer)
diff --git a/test/graded_lex.jl b/test/graded_lex.jl
new file mode 100644
index 0000000..4be60be
--- /dev/null
+++ b/test/graded_lex.jl
@@ -0,0 +1,24 @@
+# This file is a part of StarAlgebras.jl. License is MIT: https://github.com/JuliaAlgebra/StarAlgebras.jl/blob/main/LICENSE
+# Copyright (c) 2021-2025: Marek Kaluba, Benoît Legat
+
+# Example with Graded Lex Ordering
+using Test
+import StarAlgebras as SA
+
+@testset "Graded Lex" begin
+    alg = bivariate_algebra()
+    o = one(alg)
+    @test isone(o)
+    @test SA.coeffs(o).isless == grlex
+    a = SA.AlgebraElement(
+        SA.SparseCoefficients(
+            collect(Iterators.take(SA.object(SA.basis(alg)), 3)),
+            [2, -1, 3],
+            grlex,
+        ),
+        alg,
+    )
+    c = a * a
+    @test c.coeffs.values == [4, -4, 12, 1, -6, 9]
+    @test c.coeffs.basis_elements == [(0, 0), (0, 1), (1, 0), (0, 2), (1, 1), (2, 0)]
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 146f64e..07f31bd 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -46,8 +46,7 @@ end
 
     # some applications:
     using Groups
-    function Base.isless(g::Groups.FPGroupElement, h::Groups.FPGroupElement)
-        return isless(Groups.word(g), Groups.word(h))
-    end
+    lexord(a, b) = isless(Groups.word(a), Groups.word(b))
+    SA.comparable(::Type{<:Groups.FPGroupElement}) = lexord
     include("sum_of_squares.jl")
 end