From 636aeca3feda96a5f08463e2bbcab441b5e1c418 Mon Sep 17 00:00:00 2001
From: Alexis Montoison <alexis.montoison@polymtl.ca>
Date: Mon, 12 Jan 2026 17:55:45 -0600
Subject: [PATCH 1/7] Add a function substitutable_columns

---
 docs/src/dev.md |   1 +
 src/check.jl    | 128 +++++++++++++++++++++++++++++++++++++++++++++++-
 test/check.jl   |  94 +++++++++++++++++++++++++++++++++++
 test/utils.jl   |  37 +++++++++++++-
 4 files changed, 257 insertions(+), 3 deletions(-)

diff --git a/docs/src/dev.md b/docs/src/dev.md
index 82579dd6..cfb11f56 100644
--- a/docs/src/dev.md
+++ b/docs/src/dev.md
@@ -50,6 +50,7 @@ SparseMatrixColorings.directly_recoverable_columns
 SparseMatrixColorings.symmetrically_orthogonal_columns
 SparseMatrixColorings.structurally_orthogonal_columns
 SparseMatrixColorings.structurally_biorthogonal
+SparseMatrixColorings.substitutable_columns
 SparseMatrixColorings.valid_dynamic_order
 ```
 
diff --git a/src/check.jl b/src/check.jl
index f47f0958..0e68463f 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(
@@ -261,6 +265,128 @@ function directly_recoverable_columns(
     return true
 end
 
+"""
+    substitutable_columns(
+        A::AbstractMatrix, order_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]`, `order_nonzeros[i, j]` provides its order 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, order_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, order_nonzeros; i, j, ci, cj, row_group=group, column_group=group, verbose
+        )
+        !check && return false
+    end
+    return true
+end
+
+function _substitutable_check(
+    A::AbstractMatrix,
+    order_nonzeros::AbstractMatrix;
+    i::Integer,
+    j::Integer,
+    ci::Integer,
+    cj::Integer,
+    row_group::AbstractVector,
+    column_group::AbstractVector,
+    verbose::Bool)
+    order_ij = order_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 = order_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 = order_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)
diff --git a/test/check.jl b/test/check.jl
index 4e9e5a61..2ce7521f 100644
--- a/test/check.jl
+++ b/test/check.jl
@@ -4,6 +4,7 @@ using SparseMatrixColorings:
     symmetrically_orthogonal_columns,
     structurally_biorthogonal,
     directly_recoverable_columns,
+    substitutable_columns,
     what_fig_41,
     efficient_fig_1
 using Test
@@ -184,3 +185,96 @@ 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
+
+    substitutable_columns(A1, B1, [1, 2, 2, 2, 2])
+    substitutable_columns(A2, B2, [1, 2, 3, 1, 2])
+    substitutable_columns(A3, B3, [1, 2, 3, 4, 0])
+
+    # failure
+
+    @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
diff --git a/test/utils.jl b/test/utils.jl
index 2462f8c1..159d599d 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -7,18 +7,45 @@ using SparseMatrixColorings
 using SparseMatrixColorings:
     AdjacencyGraph,
     LinearSystemColoringResult,
+    TreeSetColoringResult,
     directly_recoverable_columns,
     matrix_versions,
     respectful_similar,
     structurally_orthogonal_columns,
     symmetrically_orthogonal_columns,
-    structurally_biorthogonal
+    structurally_biorthogonal,
+    substitutable_columns
 using Test
 
 const _ALL_ORDERS = (
     NaturalOrder(), LargestFirst(), SmallestLast(), IncidenceDegree(), DynamicLargestFirst()
 )
 
+function order_from_trees(result::TreeSetColoringResult)
+    (; ag, color, reverse_bfs_orders, diagonal_indices, tree_edge_indices, nt) = result
+    (; S) = ag
+    n = length(color)
+    nnzS = nnz(S)
+    nzval = zeros(Int, nnzS)
+    order_nonzeros = SparseMatrixCSC(n, n, S.colptr, S.rowval, nzval)
+    counter = 0
+    for i in diagonal_indices
+        counter += 1
+        order_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
+            order_nonzeros[i, j] = counter
+            order_nonzeros[j, i] = counter
+        end
+    end
+    return order_nonzeros
+end
+
 function test_coloring_decompression(
     A0::AbstractMatrix,
     problem::ColoringProblem{structure,partition},
@@ -77,7 +104,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 +123,13 @@ function test_coloring_decompression(
             end
         end
 
+        if decompression == :substitution
+            if structure == :symmetric
+                order_nonzeros = order_from_trees(result)
+                @test substitutable_columns(A0, order_nonzeros, color)
+            end
+        end
+
         @testset "Single-color decompression" begin
             if decompression == :direct  # TODO: implement for :substitution too
                 A2 = respectful_similar(A, eltype(B))

From b3bf3c05c66f146c8283847f60cfb82b0f0b5057 Mon Sep 17 00:00:00 2001
From: Alexis Montoison <alexis.montoison@polymtl.ca>
Date: Mon, 12 Jan 2026 18:05:30 -0600
Subject: [PATCH 2/7] Use the formatter

---
 src/check.jl  | 10 +++++++---
 test/check.jl |  6 +++---
 2 files changed, 10 insertions(+), 6 deletions(-)

diff --git a/src/check.jl b/src/check.jl
index 0e68463f..76ba4dd4 100644
--- a/src/check.jl
+++ b/src/check.jl
@@ -291,7 +291,10 @@ It is equivalent to an __acyclic coloring__.
 > [_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, order_nonzeros::AbstractMatrix, color::AbstractVector{<:Integer}; verbose::Bool=false
+    A::AbstractMatrix,
+    order_nonzeros::AbstractMatrix,
+    color::AbstractVector{<:Integer};
+    verbose::Bool=false,
 )
     checksquare(A)
     if !proper_length_coloring(A, color; verbose)
@@ -319,8 +322,9 @@ function _substitutable_check(
     cj::Integer,
     row_group::AbstractVector,
     column_group::AbstractVector,
-    verbose::Bool)
-    order_ij = order_nonzeros[i,j]
+    verbose::Bool,
+)
+    order_ij = order_nonzeros[i, j]
     k_row = 0
     k_column = 0
     if ci != 0
diff --git a/test/check.jl b/test/check.jl
index 2ce7521f..172cf418 100644
--- a/test/check.jl
+++ b/test/check.jl
@@ -223,9 +223,9 @@ end
         1 1 1 1 0
     ]
     B3 = [
-        0 1 2 3  4
-        1 0 5 6  7
-        2 5 0 8  9
+        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
     ]

From ee27403a86a79f7dbf15d11009d56cf21395d067 Mon Sep 17 00:00:00 2001
From: Alexis Montoison <alexis.montoison@polymtl.ca>
Date: Mon, 12 Jan 2026 18:36:46 -0600
Subject: [PATCH 3/7] Fix the code coverage

---
 test/check.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/test/check.jl b/test/check.jl
index 172cf418..589e6f97 100644
--- a/test/check.jl
+++ b/test/check.jl
@@ -238,6 +238,10 @@ end
 
     # 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,

From 19623f838963d7792d1e2428d248a5ed0934dc76 Mon Sep 17 00:00:00 2001
From: Alexis Montoison <alexis.montoison@polymtl.ca>
Date: Tue, 13 Jan 2026 13:11:33 -0600
Subject: [PATCH 4/7] Add a function substitutable_bidirectional

---
 docs/src/dev.md |  1 +
 src/check.jl    | 46 ++++++++++++++++++++++++++++++++++++++-
 test/check.jl   | 25 ++++++++++++++++++++++
 test/utils.jl   | 57 +++++++++++++++++++++++++++++++++++++++++++------
 4 files changed, 121 insertions(+), 8 deletions(-)

diff --git a/docs/src/dev.md b/docs/src/dev.md
index cfb11f56..55df94d5 100644
--- a/docs/src/dev.md
+++ b/docs/src/dev.md
@@ -51,6 +51,7 @@ SparseMatrixColorings.symmetrically_orthogonal_columns
 SparseMatrixColorings.structurally_orthogonal_columns
 SparseMatrixColorings.structurally_biorthogonal
 SparseMatrixColorings.substitutable_columns
+SparseMatrixColorings.substitutable_bidirectional
 SparseMatrixColorings.valid_dynamic_order
 ```
 
diff --git a/src/check.jl b/src/check.jl
index 76ba4dd4..35bdfe42 100644
--- a/src/check.jl
+++ b/src/check.jl
@@ -130,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 an __star bicoloring__.
+
 !!! warning
     This function is not coded with efficiency in mind, it is designed for small-scale tests.
 """
@@ -142,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]
@@ -313,6 +315,48 @@ function substitutable_columns(
     return true
 end
 
+"""
+    substitutable_bidirectional(
+        A::AbstractMatrix, order_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]`, `order_nonzeros[i, j]` provides its order 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,
+    order_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, order_nonzeros; i, j, ci, cj, row_group, column_group, verbose
+        )
+        !check && return false
+    end
+    return true
+end
+
 function _substitutable_check(
     A::AbstractMatrix,
     order_nonzeros::AbstractMatrix;
diff --git a/test/check.jl b/test/check.jl
index 589e6f97..45cb33b4 100644
--- a/test/check.jl
+++ b/test/check.jl
@@ -282,3 +282,28 @@ For coefficient (i=3, j=5) with colors (ci=3, cj=0):
 """,
     ) 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
+
+    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 159d599d..2595d0c3 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -8,13 +8,15 @@ using SparseMatrixColorings:
     AdjacencyGraph,
     LinearSystemColoringResult,
     TreeSetColoringResult,
+    BicoloringResult,
     directly_recoverable_columns,
     matrix_versions,
     respectful_similar,
     structurally_orthogonal_columns,
     symmetrically_orthogonal_columns,
     structurally_biorthogonal,
-    substitutable_columns
+    substitutable_columns,
+    substitutable_bidirectional
 using Test
 
 const _ALL_ORDERS = (
@@ -22,9 +24,9 @@ const _ALL_ORDERS = (
 )
 
 function order_from_trees(result::TreeSetColoringResult)
-    (; ag, color, reverse_bfs_orders, diagonal_indices, tree_edge_indices, nt) = result
+    (; A, ag, reverse_bfs_orders, diagonal_indices, tree_edge_indices, nt) = result
     (; S) = ag
-    n = length(color)
+    m, n = size(A)
     nnzS = nnz(S)
     nzval = zeros(Int, nnzS)
     order_nonzeros = SparseMatrixCSC(n, n, S.colptr, S.rowval, nzval)
@@ -46,6 +48,36 @@ function order_from_trees(result::TreeSetColoringResult)
     return order_nonzeros
 end
 
+function order_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
+    order_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
+                order_nonzeros[i - n, j] = counter
+            else
+                order_nonzeros[j - n, i] = counter
+            end
+        end
+    end
+    return order_nonzeros
+end
+
 function test_coloring_decompression(
     A0::AbstractMatrix,
     problem::ColoringProblem{structure,partition},
@@ -123,10 +155,12 @@ function test_coloring_decompression(
             end
         end
 
-        if decompression == :substitution
-            if structure == :symmetric
-                order_nonzeros = order_from_trees(result)
-                @test substitutable_columns(A0, order_nonzeros, color)
+        @testset "Substitutable" begin
+            if decompression == :substitution
+                if structure == :symmetric
+                    order_nonzeros = order_from_trees(result)
+                    @test substitutable_columns(A0, order_nonzeros, color)
+                end
             end
         end
 
@@ -268,6 +302,15 @@ function test_bicoloring_decompression(
                 @test structurally_biorthogonal(A0, row_color, column_color)
             end
         end
+
+        if decompression == :substitution
+            @testset "Substitutable" begin
+                order_nonzeros = order_from_trees(result)
+                @test substitutable_bidirectional(
+                    A0, order_nonzeros, row_color, column_color
+                )
+            end
+        end
     end
 
     @testset "More orders is better" begin

From de6bdde471ce8b256b52ef9ca621de6ea64963c2 Mon Sep 17 00:00:00 2001
From: Alexis Montoison <alexis.montoison@polymtl.ca>
Date: Tue, 13 Jan 2026 13:15:23 -0600
Subject: [PATCH 5/7] Update test/check.jl

---
 test/check.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/check.jl b/test/check.jl
index 45cb33b4..857752fe 100644
--- a/test/check.jl
+++ b/test/check.jl
@@ -5,6 +5,7 @@ using SparseMatrixColorings:
     structurally_biorthogonal,
     directly_recoverable_columns,
     substitutable_columns,
+    substitutable_bidirectional,
     what_fig_41,
     efficient_fig_1
 using Test

From 572329d3ab118d493949809f752eabede6cd6b3f Mon Sep 17 00:00:00 2001
From: Alexis Montoison <35051714+amontoison@users.noreply.github.com>
Date: Wed, 14 Jan 2026 07:20:48 -0600
Subject: [PATCH 6/7] Apply suggestions from code review

Co-authored-by: Guillaume Dalle <22795598+gdalle@users.noreply.github.com>
---
 src/check.jl  | 2 +-
 test/check.jl | 8 ++++----
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/check.jl b/src/check.jl
index 35bdfe42..55c5b35a 100644
--- a/src/check.jl
+++ b/src/check.jl
@@ -130,7 +130,7 @@ 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 an __star bicoloring__.
+It is equivalent to a __star bicoloring__.
 
 !!! warning
     This function is not coded with efficiency in mind, it is designed for small-scale tests.
diff --git a/test/check.jl b/test/check.jl
index 857752fe..b2701095 100644
--- a/test/check.jl
+++ b/test/check.jl
@@ -233,9 +233,9 @@ end
 
     # success
 
-    substitutable_columns(A1, B1, [1, 2, 2, 2, 2])
-    substitutable_columns(A2, B2, [1, 2, 3, 1, 2])
-    substitutable_columns(A3, B3, [1, 2, 3, 4, 0])
+    @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
 
@@ -298,7 +298,7 @@ end
 
     # success
 
-    substitutable_bidirectional(A, B, [1, 0, 0], [0, 1, 1])
+    @test substitutable_bidirectional(A, B, [1, 0, 0], [0, 1, 1])
 
     # failure
 

From 2a36a59f13e6857bc3739f2644c2b229a80cf045 Mon Sep 17 00:00:00 2001
From: Alexis Montoison <alexis.montoison@polymtl.ca>
Date: Wed, 14 Jan 2026 08:46:52 -0600
Subject: [PATCH 7/7] Take into account the comments of Guillaume

---
 docs/src/dev.md |  1 +
 src/check.jl    | 96 ++++++++++++++++++++++++++++++++++++++++++-------
 test/utils.jl   | 68 ++++-------------------------------
 3 files changed, 91 insertions(+), 74 deletions(-)

diff --git a/docs/src/dev.md b/docs/src/dev.md
index 55df94d5..a9e0d303 100644
--- a/docs/src/dev.md
+++ b/docs/src/dev.md
@@ -52,6 +52,7 @@ 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 55c5b35a..db5c0a34 100644
--- a/src/check.jl
+++ b/src/check.jl
@@ -269,12 +269,12 @@ end
 
 """
     substitutable_columns(
-        A::AbstractMatrix, order_nonzeros::AbstractMatrix, color::AbstractVector{<:Integer};
+        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]`, `order_nonzeros[i, j]` provides its order of recovery.
+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:
 
@@ -294,7 +294,7 @@ It is equivalent to an __acyclic coloring__.
 """
 function substitutable_columns(
     A::AbstractMatrix,
-    order_nonzeros::AbstractMatrix,
+    rank_nonzeros::AbstractMatrix,
     color::AbstractVector{<:Integer};
     verbose::Bool=false,
 )
@@ -308,7 +308,7 @@ function substitutable_columns(
         iszero(A[i, j]) && continue
         ci, cj = color[i], color[j]
         check = _substitutable_check(
-            A, order_nonzeros; i, j, ci, cj, row_group=group, column_group=group, verbose
+            A, rank_nonzeros; i, j, ci, cj, row_group=group, column_group=group, verbose
         )
         !check && return false
     end
@@ -317,12 +317,12 @@ end
 
 """
     substitutable_bidirectional(
-        A::AbstractMatrix, order_nonzeros::AbstractMatrix, row_color::AbstractVector{<:Integer}, column_color::AbstractVector{<:Integer};
+        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]`, `order_nonzeros[i, j]` provides its order of recovery.
+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:
 
@@ -336,7 +336,7 @@ It is equivalent to an __acyclic bicoloring__.
 """
 function substitutable_bidirectional(
     A::AbstractMatrix,
-    order_nonzeros::AbstractMatrix,
+    rank_nonzeros::AbstractMatrix,
     row_color::AbstractVector{<:Integer},
     column_color::AbstractVector{<:Integer};
     verbose::Bool=false,
@@ -350,7 +350,7 @@ function substitutable_bidirectional(
         iszero(A[i, j]) && continue
         ci, cj = row_color[i], column_color[j]
         check = _substitutable_check(
-            A, order_nonzeros; i, j, ci, cj, row_group, column_group, verbose
+            A, rank_nonzeros; i, j, ci, cj, row_group, column_group, verbose
         )
         !check && return false
     end
@@ -359,7 +359,7 @@ end
 
 function _substitutable_check(
     A::AbstractMatrix,
-    order_nonzeros::AbstractMatrix;
+    rank_nonzeros::AbstractMatrix;
     i::Integer,
     j::Integer,
     ci::Integer,
@@ -368,14 +368,14 @@ function _substitutable_check(
     column_group::AbstractVector,
     verbose::Bool,
 )
-    order_ij = order_nonzeros[i, j]
+    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 = order_nonzeros[k, j]
+                order_kj = rank_nonzeros[k, j]
                 @assert !iszero(order_kj)
                 if order_kj > order_ij
                     k_row = k
@@ -387,7 +387,7 @@ function _substitutable_check(
         for k in column_group[cj]
             (k == j) && continue
             if !iszero(A[i, k])
-                order_ik = order_nonzeros[i, k]
+                order_ik = rank_nonzeros[i, k]
                 @assert !iszero(order_ik)
                 if order_ik > order_ij
                     k_column = k
@@ -500,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/utils.jl b/test/utils.jl
index 2595d0c3..77c88610 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -7,8 +7,6 @@ using SparseMatrixColorings
 using SparseMatrixColorings:
     AdjacencyGraph,
     LinearSystemColoringResult,
-    TreeSetColoringResult,
-    BicoloringResult,
     directly_recoverable_columns,
     matrix_versions,
     respectful_similar,
@@ -16,68 +14,14 @@ using SparseMatrixColorings:
     symmetrically_orthogonal_columns,
     structurally_biorthogonal,
     substitutable_columns,
-    substitutable_bidirectional
+    substitutable_bidirectional,
+    rank_nonzeros_from_trees
 using Test
 
 const _ALL_ORDERS = (
     NaturalOrder(), LargestFirst(), SmallestLast(), IncidenceDegree(), DynamicLargestFirst()
 )
 
-function order_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)
-    order_nonzeros = SparseMatrixCSC(n, n, S.colptr, S.rowval, nzval)
-    counter = 0
-    for i in diagonal_indices
-        counter += 1
-        order_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
-            order_nonzeros[i, j] = counter
-            order_nonzeros[j, i] = counter
-        end
-    end
-    return order_nonzeros
-end
-
-function order_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
-    order_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
-                order_nonzeros[i - n, j] = counter
-            else
-                order_nonzeros[j - n, i] = counter
-            end
-        end
-    end
-    return order_nonzeros
-end
-
 function test_coloring_decompression(
     A0::AbstractMatrix,
     problem::ColoringProblem{structure,partition},
@@ -158,8 +102,8 @@ function test_coloring_decompression(
         @testset "Substitutable" begin
             if decompression == :substitution
                 if structure == :symmetric
-                    order_nonzeros = order_from_trees(result)
-                    @test substitutable_columns(A0, order_nonzeros, color)
+                    rank_nonzeros = rank_nonzeros_from_trees(result)
+                    @test substitutable_columns(A0, rank_nonzeros, color)
                 end
             end
         end
@@ -305,9 +249,9 @@ function test_bicoloring_decompression(
 
         if decompression == :substitution
             @testset "Substitutable" begin
-                order_nonzeros = order_from_trees(result)
+                rank_nonzeros = rank_nonzeros_from_trees(result)
                 @test substitutable_bidirectional(
-                    A0, order_nonzeros, row_color, column_color
+                    A0, rank_nonzeros, row_color, column_color
                 )
             end
         end