Add a function substitutable_columns

Author: amontoison

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
 ```
 

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)

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

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))