diff --git a/docs/src/dev.md b/docs/src/dev.md
index 82579dd6..a9e0d303 100644
--- a/docs/src/dev.md
+++ b/docs/src/dev.md
@@ -50,6 +50,9 @@ SparseMatrixColorings.directly_recoverable_columns
 SparseMatrixColorings.symmetrically_orthogonal_columns
 SparseMatrixColorings.structurally_orthogonal_columns
 SparseMatrixColorings.structurally_biorthogonal
+SparseMatrixColorings.substitutable_columns
+SparseMatrixColorings.substitutable_bidirectional
+SparseMatrixColorings.rank_nonzeros_from_trees
 SparseMatrixColorings.valid_dynamic_order
 ```
 
diff --git a/src/check.jl b/src/check.jl
index f47f0958..db5c0a34 100644
--- a/src/check.jl
+++ b/src/check.jl
@@ -86,11 +86,15 @@ A partition of the columns of a symmetric matrix `A` is _symmetrically orthogona
 1. the group containing the column `A[:, j]` has no other column with a nonzero in row `i`
 2. the group containing the column `A[:, i]` has no other column with a nonzero in row `j`
 
+It is equivalent to a __star coloring__.
+
 !!! warning
     This function is not coded with efficiency in mind, it is designed for small-scale tests.
 
 # References
 
+> [_On the Estimation of Sparse Hessian Matrices_](https://doi.org/10.1137/0716078), Powell and Toint (1979)
+> [_Estimation of sparse hessian matrices and graph coloring problems_](https://doi.org/10.1007/BF02612334), Coleman and Moré (1984)
 > [_What Color Is Your Jacobian? Graph Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005)
 """
 function symmetrically_orthogonal_columns(
@@ -102,7 +106,7 @@ function symmetrically_orthogonal_columns(
     end
     issymmetric(A) || return false
     group = group_by_color(color)
-    for i in axes(A, 2), j in axes(A, 2)
+    for i in axes(A, 1), j in axes(A, 2)
         iszero(A[i, j]) && continue
         ci, cj = color[i], color[j]
         check = _bilateral_check(
@@ -126,6 +130,8 @@ A bipartition of the rows and columns of a matrix `A` is _structurally biorthogo
 1. the group containing the column `A[:, j]` has no other column with a nonzero in row `i`
 2. the group containing the row `A[i, :]` has no other row with a nonzero in column `j`
 
+It is equivalent to a __star bicoloring__.
+
 !!! warning
     This function is not coded with efficiency in mind, it is designed for small-scale tests.
 """
@@ -138,8 +144,8 @@ function structurally_biorthogonal(
     if !proper_length_bicoloring(A, row_color, column_color; verbose)
         return false
     end
-    column_group = group_by_color(column_color)
     row_group = group_by_color(row_color)
+    column_group = group_by_color(column_color)
     for i in axes(A, 1), j in axes(A, 2)
         iszero(A[i, j]) && continue
         ci, cj = row_color[i], column_color[j]
@@ -261,6 +267,174 @@ function directly_recoverable_columns(
     return true
 end
 
+"""
+    substitutable_columns(
+        A::AbstractMatrix, rank_nonzeros::AbstractMatrix, color::AbstractVector{<:Integer};
+        verbose=false
+    )
+
+Return `true` if coloring the columns of the symmetric matrix `A` with the vector `color` results in a partition that is substitutable, and `false` otherwise.
+For all nonzeros `A[i, j]`, `rank_nonzeros[i, j]` provides its rank of recovery.
+
+A partition of the columns of a symmetric matrix `A` is _substitutable_ if, for every nonzero element `A[i, j]`, either of the following statements holds:
+
+1. the group containing the column `A[:, j]` has all nonzeros in row `i` ordered before `A[i, j]`
+2. the group containing the column `A[:, i]` has all nonzeros in row `j` ordered before `A[i, j]`
+
+It is equivalent to an __acyclic coloring__.
+
+!!! warning
+    This function is not coded with efficiency in mind, it is designed for small-scale tests.
+
+# References
+
+> [_On the Estimation of Sparse Hessian Matrices_](https://doi.org/10.1137/0716078), Powell and Toint (1979)
+> [_The Cyclic Coloring Problem and Estimation of Sparse Hessian Matrices_](https://doi.org/10.1137/0607026), Coleman and Cai (1986)
+> [_What Color Is Your Jacobian? Graph Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005)
+"""
+function substitutable_columns(
+    A::AbstractMatrix,
+    rank_nonzeros::AbstractMatrix,
+    color::AbstractVector{<:Integer};
+    verbose::Bool=false,
+)
+    checksquare(A)
+    if !proper_length_coloring(A, color; verbose)
+        return false
+    end
+    issymmetric(A) || return false
+    group = group_by_color(color)
+    for i in axes(A, 1), j in axes(A, 2)
+        iszero(A[i, j]) && continue
+        ci, cj = color[i], color[j]
+        check = _substitutable_check(
+            A, rank_nonzeros; i, j, ci, cj, row_group=group, column_group=group, verbose
+        )
+        !check && return false
+    end
+    return true
+end
+
+"""
+    substitutable_bidirectional(
+        A::AbstractMatrix, rank_nonzeros::AbstractMatrix, row_color::AbstractVector{<:Integer}, column_color::AbstractVector{<:Integer};
+        verbose=false
+    )
+
+Return `true` if bicoloring of the matrix `A` with the vectors `row_color` and `column_color` results in a bipartition that is substitutable, and `false` otherwise.
+For all nonzeros `A[i, j]`, `rank_nonzeros[i, j]` provides its rank of recovery.
+
+A bipartition of the rows and columns of a matrix `A` is _substitutable_ if, for every nonzero element `A[i, j]`, either of the following statements holds:
+
+1. the group containing the column `A[:, j]` has all nonzeros in row `i` ordered before `A[i, j]`
+2. the group containing the row `A[i, :]` has all nonzeros in column `j` ordered before `A[i, j]`
+
+It is equivalent to an __acyclic bicoloring__.
+
+!!! warning
+    This function is not coded with efficiency in mind, it is designed for small-scale tests.
+"""
+function substitutable_bidirectional(
+    A::AbstractMatrix,
+    rank_nonzeros::AbstractMatrix,
+    row_color::AbstractVector{<:Integer},
+    column_color::AbstractVector{<:Integer};
+    verbose::Bool=false,
+)
+    if !proper_length_bicoloring(A, row_color, column_color; verbose)
+        return false
+    end
+    row_group = group_by_color(row_color)
+    column_group = group_by_color(column_color)
+    for i in axes(A, 1), j in axes(A, 2)
+        iszero(A[i, j]) && continue
+        ci, cj = row_color[i], column_color[j]
+        check = _substitutable_check(
+            A, rank_nonzeros; i, j, ci, cj, row_group, column_group, verbose
+        )
+        !check && return false
+    end
+    return true
+end
+
+function _substitutable_check(
+    A::AbstractMatrix,
+    rank_nonzeros::AbstractMatrix;
+    i::Integer,
+    j::Integer,
+    ci::Integer,
+    cj::Integer,
+    row_group::AbstractVector,
+    column_group::AbstractVector,
+    verbose::Bool,
+)
+    order_ij = rank_nonzeros[i, j]
+    k_row = 0
+    k_column = 0
+    if ci != 0
+        for k in row_group[ci]
+            (k == i) && continue
+            if !iszero(A[k, j])
+                order_kj = rank_nonzeros[k, j]
+                @assert !iszero(order_kj)
+                if order_kj > order_ij
+                    k_row = k
+                end
+            end
+        end
+    end
+    if cj != 0
+        for k in column_group[cj]
+            (k == j) && continue
+            if !iszero(A[i, k])
+                order_ik = rank_nonzeros[i, k]
+                @assert !iszero(order_ik)
+                if order_ik > order_ij
+                    k_column = k
+                end
+            end
+        end
+    end
+    if ci == 0 && cj == 0
+        if verbose
+            @warn """
+            For coefficient (i=$i, j=$j) with colors (ci=$ci, cj=$cj):
+            - Row color ci=$ci is neutral.
+            - Column color cj=$cj is neutral.
+            """
+        end
+        return false
+    elseif ci == 0 && !iszero(k_column)
+        if verbose
+            @warn """
+            For coefficient (i=$i, j=$j) with colors (ci=$ci, cj=$cj):
+            - Row color ci=$ci is neutral.
+            - For the column $k_column in column color cj=$cj, A[$i, $k_column] is ordered after A[$i, $j].
+            """
+        end
+        return false
+    elseif cj == 0 && !iszero(k_row)
+        if verbose
+            @warn """
+            For coefficient (i=$i, j=$j) with colors (ci=$ci, cj=$cj):
+            - For the row $k_row in row color ci=$ci, A[$k_row, $j] is ordered after A[$i, $j].
+            - Column color cj=$cj is neutral.
+            """
+        end
+        return false
+    elseif !iszero(k_row) && !iszero(k_column)
+        if verbose
+            @warn """
+            For coefficient (i=$i, j=$j) with colors (ci=$ci, cj=$cj):
+            - For the row $k_row in row color ci=$ci, A[$k_row, $j] is ordered after A[$i, $j].
+            - For the column $k_column in column color cj=$cj, A[$i, $k_column] is ordered after A[$i, $j].
+            """
+        end
+        return false
+    end
+    return true
+end
+
 """
     valid_dynamic_order(g::AdjacencyGraph, π::AbstractVector{<:Integer}, order::DynamicDegreeBasedOrder)
     valid_dynamic_order(bg::AdjacencyGraph, ::Val{side}, π::AbstractVector{<:Integer}, order::DynamicDegreeBasedOrder)
@@ -326,3 +500,75 @@ function valid_dynamic_order(
     end
     return true
 end
+
+"""
+    rank_nonzeros_from_trees(result::TreeSetColoringResult)
+    rank_nonzeros_from_trees(result::BicoloringResult)
+
+Construct a sparse matrix `rank_nonzeros` that assigns a unique recovery rank
+to each nonzero coefficient associated with an acyclic coloring or bicoloring.
+
+For every nonzero entry `result.A[i, j]`, `rank_nonzeros[i, j]` stores a strictly positive
+integer representing the order in which this coefficient is recovered during the decompression.
+A larger value means the coefficient is recovered later.
+
+This ranking is used to test substitutability (acyclicity) of colorings:
+for a given nonzero `result.A[i, j]`, the ranks allow one to check whether all competing
+nonzeros in the same row or column (within a color group) are recovered before it.
+"""
+function rank_nonzeros_from_trees end
+
+function rank_nonzeros_from_trees(result::TreeSetColoringResult)
+    (; A, ag, reverse_bfs_orders, diagonal_indices, tree_edge_indices, nt) = result
+    (; S) = ag
+    m, n = size(A)
+    nnzS = nnz(S)
+    nzval = zeros(Int, nnzS)
+    rank_nonzeros = SparseMatrixCSC(n, n, S.colptr, S.rowval, nzval)
+    counter = 0
+    for i in diagonal_indices
+        counter += 1
+        rank_nonzeros[i, i] = counter
+    end
+    for k in 1:nt
+        first = tree_edge_indices[k]
+        last = tree_edge_indices[k + 1] - 1
+        for pos in first:last
+            (i, j) = reverse_bfs_orders[pos]
+            counter += 1
+            rank_nonzeros[i, j] = counter
+            rank_nonzeros[j, i] = counter
+        end
+    end
+    return rank_nonzeros
+end
+
+function rank_nonzeros_from_trees(result::BicoloringResult)
+    (; A, abg, row_color, column_color, symmetric_result, large_colptr, large_rowval) =
+        result
+    @assert symmetric_result isa TreeSetColoringResult
+    (; ag, reverse_bfs_orders, tree_edge_indices, nt) = symmetric_result
+    (; S) = ag
+    m, n = size(A)
+    nnzA = nnz(S) ÷ 2
+    nzval = zeros(Int, nnzA)
+    colptr = large_colptr[1:(n + 1)]
+    rowval = large_rowval[1:nnzA]
+    rowval .-= n
+    rank_nonzeros = SparseMatrixCSC(m, n, colptr, rowval, nzval)
+    counter = 0
+    for k in 1:nt
+        first = tree_edge_indices[k]
+        last = tree_edge_indices[k + 1] - 1
+        for pos in first:last
+            (i, j) = reverse_bfs_orders[pos]
+            counter += 1
+            if i > j
+                rank_nonzeros[i - n, j] = counter
+            else
+                rank_nonzeros[j - n, i] = counter
+            end
+        end
+    end
+    return rank_nonzeros
+end
diff --git a/test/check.jl b/test/check.jl
index 4e9e5a61..b2701095 100644
--- a/test/check.jl
+++ b/test/check.jl
@@ -4,6 +4,8 @@ using SparseMatrixColorings:
     symmetrically_orthogonal_columns,
     structurally_biorthogonal,
     directly_recoverable_columns,
+    substitutable_columns,
+    substitutable_bidirectional,
     what_fig_41,
     efficient_fig_1
 using Test
@@ -184,3 +186,125 @@ For coefficient (i=1, j=1) with colors (ci=0, cj=0):
 """,
     ) structurally_biorthogonal(A, [0, 2, 2, 3], [0, 2, 2, 2, 3], verbose=true)
 end
+
+@testset "Substitutable columns" begin
+    A1 = [
+        1 1 1 1 1
+        1 1 0 0 0
+        1 0 1 0 0
+        1 0 0 1 0
+        1 0 0 0 1
+    ]
+    B1 = [
+        1 6 7 8 9
+        6 2 0 0 0
+        7 0 3 0 0
+        8 0 0 4 0
+        9 0 0 0 5
+    ]
+    A2 = [
+        1 1 0 0 0
+        1 1 1 0 0
+        0 1 1 1 0
+        0 0 1 1 1
+        0 0 0 1 1
+    ]
+    B2 = [
+        5 1 0 0 0
+        1 6 2 0 0
+        0 2 7 3 0
+        0 0 3 8 4
+        0 0 0 4 9
+    ]
+    A3 = [
+        0 1 1 1 1
+        1 0 1 1 1
+        1 1 0 1 1
+        1 1 1 0 1
+        1 1 1 1 0
+    ]
+    B3 = [
+        0 1 2 3 4
+        1 0 5 6 7
+        2 5 0 8 9
+        3 6 8 0 10
+        4 7 9 10 0
+    ]
+
+    # success
+
+    @test substitutable_columns(A1, B1, [1, 2, 2, 2, 2])
+    @test substitutable_columns(A2, B2, [1, 2, 3, 1, 2])
+    @test substitutable_columns(A3, B3, [1, 2, 3, 4, 0])
+
+    # failure
+
+    @test !substitutable_columns(A1, B1, [1, 1, 1, 1])
+    log = (:warn, "4 colors provided for 5 columns.")
+    @test_logs log substitutable_columns(A1, B1, [1, 1, 1, 1]; verbose=true)
+
+    @test !substitutable_columns(A1, B1, [1, 1, 1, 1, 1])
+    @test_logs (
+        :warn,
+        """
+For coefficient (i=1, j=1) with colors (ci=1, cj=1):
+- For the row 5 in row color ci=1, A[5, 1] is ordered after A[1, 1].
+- For the column 5 in column color cj=1, A[1, 5] is ordered after A[1, 1].
+""",
+    ) substitutable_columns(A1, B1, [1, 1, 1, 1, 1]; verbose=true)
+
+    @test !substitutable_columns(A2, B2, [1, 2, 0, 1, 2])
+    @test_logs (
+        :warn,
+        """
+For coefficient (i=3, j=3) with colors (ci=0, cj=0):
+- Row color ci=0 is neutral.
+- Column color cj=0 is neutral.
+""",
+    ) substitutable_columns(A2, B2, [1, 2, 0, 1, 2]; verbose=true)
+
+    @test !substitutable_columns(A3, B3, [0, 1, 2, 3, 3])
+    @test_logs (
+        :warn,
+        """
+For coefficient (i=1, j=4) with colors (ci=0, cj=3):
+- Row color ci=0 is neutral.
+- For the column 5 in column color cj=3, A[1, 5] is ordered after A[1, 4].
+""",
+    ) substitutable_columns(A3, B3, [0, 1, 2, 3, 3]; verbose=true)
+
+    @test !substitutable_columns(A3, B3, [1, 2, 3, 3, 0])
+    @test_logs (
+        :warn,
+        """
+For coefficient (i=3, j=5) with colors (ci=3, cj=0):
+- For the row 4 in row color ci=3, A[4, 5] is ordered after A[3, 5].
+- Column color cj=0 is neutral.
+""",
+    ) substitutable_columns(A3, B3, [1, 2, 3, 3, 0]; verbose=true)
+end
+
+@testset "Substitutable bidirectional" begin
+    A = [
+        1 0 0
+        0 1 0
+        0 0 1
+    ]
+    B = [
+        1 0 0
+        0 2 0
+        0 0 3
+    ]
+
+    # success
+
+    @test substitutable_bidirectional(A, B, [1, 0, 0], [0, 1, 1])
+
+    # failure
+
+    log = (:warn, "2 colors provided for 3 columns.")
+    @test_logs log !substitutable_bidirectional(A, B, [1, 0, 0], [0, 1]; verbose=true)
+
+    log = (:warn, "4 colors provided for 3 rows.")
+    @test_logs log !substitutable_bidirectional(A, B, [1, 0, 0, 1], [0, 1, 1]; verbose=true)
+end
diff --git a/test/utils.jl b/test/utils.jl
index 2462f8c1..77c88610 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -12,7 +12,10 @@ using SparseMatrixColorings:
     respectful_similar,
     structurally_orthogonal_columns,
     symmetrically_orthogonal_columns,
-    structurally_biorthogonal
+    structurally_biorthogonal,
+    substitutable_columns,
+    substitutable_bidirectional,
+    rank_nonzeros_from_trees
 using Test
 
 const _ALL_ORDERS = (
@@ -77,7 +80,6 @@ function test_coloring_decompression(
         end
 
         @testset "Recoverability" begin
-            # TODO: find tests for recoverability for substitution decompression
             if decompression == :direct
                 if structure == :nonsymmetric
                     if partition == :column
@@ -97,6 +99,15 @@ function test_coloring_decompression(
             end
         end
 
+        @testset "Substitutable" begin
+            if decompression == :substitution
+                if structure == :symmetric
+                    rank_nonzeros = rank_nonzeros_from_trees(result)
+                    @test substitutable_columns(A0, rank_nonzeros, color)
+                end
+            end
+        end
+
         @testset "Single-color decompression" begin
             if decompression == :direct  # TODO: implement for :substitution too
                 A2 = respectful_similar(A, eltype(B))
@@ -235,6 +246,15 @@ function test_bicoloring_decompression(
                 @test structurally_biorthogonal(A0, row_color, column_color)
             end
         end
+
+        if decompression == :substitution
+            @testset "Substitutable" begin
+                rank_nonzeros = rank_nonzeros_from_trees(result)
+                @test substitutable_bidirectional(
+                    A0, rank_nonzeros, row_color, column_color
+                )
+            end
+        end
     end
 
     @testset "More orders is better" begin