diff --git a/CITATION.cff b/CITATION.cff
index 18ff414f..c2c0a8fb 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -28,3 +28,23 @@ keywords:
   - automatic differentiation
   - julia programming language
 license: MIT
+preferred-citation:
+  authors:
+    - given-names: Alexis
+      family-names: Montoison
+      orcid: 'https://orcid.org/0000-0002-3403-5450'
+    - given-names: Guillaume
+      family-names: Dalle
+      orcid: 'https://orcid.org/0000-0003-4866-1687'
+    - given-names: Assefaw
+      family-names: Gebremedhin
+      orcid: 'https://orcid.org/0000-0001-5383-8032'
+  title: "Revisiting Sparse Matrix Coloring and Bicoloring"
+  year: 2025
+  type: article
+  url: 'https://arxiv.org/abs/2505.07308'
+  identifiers:
+    - type: doi
+      value: 10.48550/arXiv.2505.07308
+      description: Arxiv
+  abstract: "Sparse matrix coloring and bicoloring are fundamental building blocks of sparse automatic differentiation. Bicoloring is particularly advantageous for rectangular Jacobian matrices with at least one dense row and column. Indeed, in such cases, unidirectional row or column coloring demands a number of colors equal to the number of rows or columns. We introduce a new strategy for bicoloring that encompasses both direct and substitution-based decompression approaches. Our method reformulates the two variants of bicoloring as star and acyclic colorings of an augmented symmetric matrix. We extend the concept of neutral colors, previously exclusive to bicoloring, to symmetric colorings, and we propose a post-processing routine that neutralizes colors to further reduce the overall color count. We also present the Julia package SparseMatrixColorings, which includes these new bicoloring algorithms alongside all standard coloring methods for sparse derivative matrix computation. Compared to ColPack, the Julia package also offers enhanced implementations for star and acyclic coloring, vertex ordering, as well as decompression."
diff --git a/README.md b/README.md
index 40b3a2d6..66b10dd3 100644
--- a/README.md
+++ b/README.md
@@ -5,6 +5,7 @@
 [![Dev Documentation](https://img.shields.io/badge/docs-dev-blue.svg)](https://gdalle.github.io/SparseMatrixColorings.jl/dev/)
 [![Coverage](https://codecov.io/gh/gdalle/SparseMatrixColorings.jl/branch/main/graph/badge.svg)](https://app.codecov.io/gh/gdalle/SparseMatrixColorings.jl)
 [![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/JuliaDiff/BlueStyle)
+[![arXiv](https://img.shields.io/badge/arXiv-2505.07308-b31b1b.svg)](https://arxiv.org/abs/2505.07308)
 [![DOI](https://zenodo.org/badge/801999408.svg)](https://zenodo.org/doi/10.5281/zenodo.11314275)
 
 Coloring algorithms for sparse Jacobian and Hessian matrices.
@@ -19,7 +20,11 @@ pkg> add SparseMatrixColorings
 
 ## Background
 
-The algorithms implemented in this package are taken from the following articles:
+The algorithms implemented in this package are described in the following preprint:
+
+- [_Revisiting Sparse Matrix Coloring and Bicoloring_](https://arxiv.org/abs/2505.07308), Montoison et al. (2025)
+
+and inspired by previous works:
 
 - [_What Color Is Your Jacobian? Graph Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005)
 - [_New Acyclic and Star Coloring Algorithms with Application to Computing Hessians_](https://epubs.siam.org/doi/abs/10.1137/050639879), Gebremedhin et al. (2007)
@@ -35,5 +40,18 @@ Some parts of the articles (like definitions) are thus copied verbatim in the do
 
 ## Citing
 
-Please cite this software using the provided `CITATION.cff` file.
+Please cite this software using the provided `CITATION.cff` file or the `.bib` entry below:
+
+```bibtex
+@unpublished{montoison2025revisitingsparsematrixcoloring,
+      title={Revisiting Sparse Matrix Coloring and Bicoloring}, 
+      author={Alexis Montoison and Guillaume Dalle and Assefaw Gebremedhin},
+      year={2025},
+      eprint={2505.07308},
+      archivePrefix={arXiv},
+      primaryClass={math.NA},
+      url={https://arxiv.org/abs/2505.07308}, 
+}
+```
+
 The link <https://zenodo.org/doi/10.5281/zenodo.11314275> resolves to the latest version on Zenodo.
diff --git a/docs/make.jl b/docs/make.jl
index 405e8082..78fba268 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -12,9 +12,14 @@ makedocs(;
     sitename="SparseMatrixColorings.jl",
     format=Documenter.HTML(),
     pages=[
-        "Home" => "index.md", "api.md", "Developer Documentation" => ["dev.md", "vis.md"]
+        "Home" => "index.md",
+        "tutorial.md",
+        "api.md",
+        "Developer Documentation" => ["dev.md", "vis.md"],
     ],
     plugins=[links],
 )
 
-deploydocs(; repo="github.com/gdalle/SparseMatrixColorings.jl", devbranch="main")
+deploydocs(;
+    repo="github.com/gdalle/SparseMatrixColorings.jl", push_preview=true, devbranch="main"
+)
diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md
new file mode 100644
index 00000000..fec54f13
--- /dev/null
+++ b/docs/src/tutorial.md
@@ -0,0 +1,323 @@
+# Tutorial
+
+Here we give a brief introduction to the contents of the package, see the [API reference](@ref) for more details.
+
+```@example tutorial
+using SparseMatrixColorings
+using LinearAlgebra
+using SparseArrays
+using StableRNGs
+using SparseMatrixColorings: show_colors # hide
+using Images # hide
+
+scale=15 # hide
+pad=3 # hide
+border=2 # hide
+nothing # hide
+```
+
+## Coloring problems and algorithms
+
+SparseMatrixColorings.jl is based on the combination of a coloring problem and a coloring algorithm, which will be passed to the [`coloring`](@ref) function.
+
+The problem defines what you want to solve. It is always a [`ColoringProblem`](@ref), and you can select options such as
+
+- the structure of the matrix (`:nonsymmetric` or `:symmetric`)
+- the type of partition you want (`:column`, `:row` or `:bidirectional`).
+
+```@example tutorial
+problem = ColoringProblem()
+```
+
+The algorithm defines how you want to solve it. It can be either a [`GreedyColoringAlgorithm`](@ref) or a [`ConstantColoringAlgorithm`](@ref). For `GreedyColoringAlgorithm`, you can select options such as
+
+- the order in which vertices are processed (a subtype of [`AbstractOrder`](@ref SparseMatrixColorings.AbstractOrder))
+- the type of decompression you want (`:direct` or `:substitution`)
+
+```@example tutorial
+algo = GreedyColoringAlgorithm()
+```
+
+## Coloring results
+
+The [`coloring`](@ref) function takes a matrix, a problem and an algorithm, to return a result subtyping [`AbstractColoringResult`](@ref).
+
+```@example tutorial
+S = sparse([
+    0 0 1 1 0 1
+    1 0 0 0 1 0
+    0 1 0 0 1 0
+    0 1 1 0 0 0
+])
+
+result = coloring(S, problem, algo)
+```
+
+The detailed type and fields of that result are _not part of the public API_.
+To access its contents, you can use the following getters:
+
+- [`sparsity_pattern`](@ref) for the matrix initially provided to `coloring`
+- [`ncolors`](@ref) for the total number of distinct colors
+- [`row_colors`](@ref), [`column_colors`](@ref) for vectors of integer colors (depending on the partition)
+- [`row_groups`](@ref), [`column_groups`](@ref) for vector of row or column indices grouped by color (depending on the partition)
+
+Here, we have a column coloring, so we can try the following:
+
+```@example tutorial
+column_colors(result)
+```
+
+```@example tutorial
+column_groups(result)
+```
+
+```@example tutorial
+ncolors(result)
+```
+
+## Compression and decompression
+
+The functions [`compress`](@ref) and [`decompress`](@ref) efficiently store and retrieve compressed representations of sparse matrices, using the coloring result as a starting point.
+
+Compression sums all columns or rows with the same color:
+
+```@example tutorial
+M = sparse([
+    0 0 4 6 0 9
+    1 0 0 0 7 0
+    0 2 0 0 8 0
+    0 3 5 0 0 0
+])
+
+B = compress(M, result)
+```
+
+Decompression recovers the original matrix from its compressed version:
+
+```@example tutorial
+C = decompress(B, result)
+```
+
+The functions [`decompress!`](@ref) and [`decompress_single_color!`](@ref) are in-place variants of [`decompress`](@ref).
+
+```@example tutorial
+D = [
+    10 14 18
+    11 15 0
+    12 16 0
+    13 17 0
+]
+
+decompress!(C, D, result)
+```
+
+```@example tutorial
+nonzeros(C) .= -1
+decompress_single_color!(C, D[:, 2], 2, result)
+```
+
+## Unidirectional variants
+
+We now illustrate the variants of colorings available, on the following matrix:
+
+```@example tutorial
+S = sparse(Symmetric(sprand(StableRNG(0), Bool, 10, 10, 0.4)))
+```
+
+We start with unidirectional colorings, where only rows or columns are colored and the matrix is not assumed to be symmetric.
+
+### Column coloring
+
+```@example tutorial
+problem = ColoringProblem(; structure=:nonsymmetric, partition=:column)
+algo = GreedyColoringAlgorithm(; decompression=:direct)
+result = coloring(S, problem, algo)
+nothing # hide
+```
+
+```@example tutorial
+ncolors(result), column_colors(result)
+```
+
+Here is the colored matrix
+
+```@example tutorial
+A_img, B_img = show_colors(result; scale, pad, border) # hide
+A_img # hide
+```
+
+and its columnwise compression
+
+```@example tutorial
+B_img # hide
+```
+
+### Row coloring
+
+```@example tutorial
+problem = ColoringProblem(; structure=:nonsymmetric, partition=:row)
+algo = GreedyColoringAlgorithm(; decompression=:direct)
+result = coloring(S, problem, algo)
+nothing # hide
+```
+
+```@example tutorial
+ncolors(result), row_colors(result)
+```
+
+Here is the colored matrix
+
+```@example tutorial
+A_img, B_img = show_colors(result; scale, pad, border) # hide
+A_img # hide
+```
+
+and its rowwise compression
+
+```@example tutorial
+B_img # hide
+```
+
+## Symmetric variants
+
+We continue with unidirectional symmetric colorings, where coloring rows is equivalent to coloring columns.
+Symmetry is leveraged to possibly reduce the number of necessary colors.
+
+### Star coloring
+
+Star coloring is the algorithm used for symmetric matrices with direct decompression.
+
+```@example tutorial
+problem = ColoringProblem(; structure=:symmetric, partition=:column)
+algo = GreedyColoringAlgorithm(; decompression=:direct)
+result = coloring(S, problem, algo)
+nothing # hide
+```
+
+```@example tutorial
+ncolors(result), column_colors(result)
+```
+
+Here is the colored matrix
+
+```@example tutorial
+A_img, B_img = show_colors(result; scale, pad, border) # hide
+A_img # hide
+```
+
+and its columnwise compression
+
+```@example tutorial
+B_img # hide
+```
+
+### Acyclic coloring
+
+Acyclic coloring is the algorithm used for symmetric matrices with decompression by substitution.
+
+```@example tutorial
+problem = ColoringProblem(; structure=:symmetric, partition=:column)
+algo = GreedyColoringAlgorithm(; decompression=:substitution)
+result = coloring(S, problem, algo)
+nothing # hide
+```
+
+```@example tutorial
+ncolors(result), column_colors(result)
+```
+
+Here is the colored matrix
+
+```@example tutorial
+A_img, B_img = show_colors(result; scale, pad, border) # hide
+A_img # hide
+```
+
+and its columnwise compression
+
+```@example tutorial
+B_img # hide
+```
+
+## Bidirectional variants
+
+We finish with bidirectional colorings, where both rows and columns are colored and the matrix is not assumed to be symmetric.
+
+Bicoloring is most relevant for matrices with dense rows and columns, which is why we consider the following test case:
+
+```@example tutorial
+S_bi = copy(S)
+S_bi[:, 1] .= true
+S_bi[1, :] .= true
+S_bi
+```
+
+With our implementations, bidirectional coloring works better using a [`RandomOrder`](@ref).
+
+### Star bicoloring
+
+```@example tutorial
+problem = ColoringProblem(; structure=:nonsymmetric, partition=:bidirectional)
+algo = GreedyColoringAlgorithm(RandomOrder(StableRNG(0), 0); decompression=:direct)
+result = coloring(S_bi, problem, algo)
+nothing # hide
+```
+
+```@example tutorial
+ncolors(result), column_colors(result)
+```
+
+Here is the colored matrix
+
+```@example tutorial
+Arc_img, _, _, Br_img, Bc_img = show_colors(result; scale, pad, border) # hide
+Arc_img # hide
+```
+
+its columnwise compression
+
+```@example tutorial
+Bc_img # hide
+```
+
+and rowwise compression
+
+```@example tutorial
+Br_img # hide
+```
+
+Both are necessary to reconstruct the original.
+
+#### Acyclic bicoloring
+
+```@example tutorial
+problem = ColoringProblem(; structure=:nonsymmetric, partition=:bidirectional)
+algo = GreedyColoringAlgorithm(RandomOrder(StableRNG(0), 0); decompression=:substitution)
+result = coloring(S_bi, problem, algo)
+nothing # hide
+```
+
+```@example tutorial
+ncolors(result), column_colors(result)
+```
+
+Here is the colored matrix
+
+```@example tutorial
+Arc_img, _, _, Br_img, Bc_img = show_colors(result; scale=20, pad, border) # hide
+Arc_img # hide
+```
+
+its columnwise compression
+
+```@example tutorial
+Bc_img # hide
+```
+
+and rowwise compression
+
+```@example tutorial
+Br_img # hide
+```
+
+Both are necessary to reconstruct the original.
diff --git a/docs/src/vis.md b/docs/src/vis.md
index 75085162..f1f9e7b3 100644
--- a/docs/src/vis.md
+++ b/docs/src/vis.md
@@ -60,7 +60,7 @@ Finally, a background color can be passed via the `background` keyword argument.
 
 We demonstrate this on a bidirectional coloring.
 
-```@example img
+```@example img; continue = true
 S = sparse([
     1 1 1 1 1 1 1 1 1
     1 0 0 0 0 0 0 0 1
diff --git a/src/decompression.jl b/src/decompression.jl
index 97aa0900..bc04394c 100644
--- a/src/decompression.jl
+++ b/src/decompression.jl
@@ -116,6 +116,9 @@ The in-place alternative is [`decompress!`](@ref).
 Compression means summing either the columns or the rows of `A` which share the same color.
 It is done by calling [`compress`](@ref).
 
+!!! warning
+    For some coloring variants, the `result` object is mutated during decompression.
+
 # Example
 
 ```jldoctest
@@ -196,6 +199,9 @@ It is done by calling [`compress`](@ref).
 For `:symmetric` coloring results (and for those only), an optional positional argument `uplo in (:U, :L, :F)` can be passed to specify which part of the matrix `A` should be updated: the Upper triangle, the Lower triangle, or the Full matrix.
 When `A isa SparseMatrixCSC`, using the `uplo` argument requires a target matrix which only stores the relevant triangle(s).
 
+!!! warning
+    For some coloring variants, the `result` object is mutated during decompression.
+
 # Example
 
 ```jldoctest
@@ -261,6 +267,9 @@ Decompress the vector `b` corresponding to color `c` in-place into `A`, given a
 For `:symmetric` coloring results (and for those only), an optional positional argument `uplo in (:U, :L, :F)` can be passed to specify which part of the matrix `A` should be updated: the Upper triangle, the Lower triangle, or the Full matrix.
 When `A isa SparseMatrixCSC`, using the `uplo` argument requires a target matrix which only stores the relevant triangle(s).
 
+!!! warning
+    For some coloring variants, the `result` object is mutated during decompression.
+
 # Example
 
 ```jldoctest