Start symmetric decompression (#10)

* Start symmetric decompression

* Tests on more matrix formats

* Add missing tests

Author: gdalle

Project.toml

@@ -1,7 +1,7 @@
 name = "SparseMatrixColorings"
 uuid = "0a514795-09f3-496d-8182-132a7b665d35"
 authors = ["Guillaume Dalle <22795598+gdalle@users.noreply.github.com>"]
-version = "0.3.1"
+version = "0.3.2"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"

docs/src/api.md

@@ -26,15 +26,25 @@ LargestFirst
 ### Decompression
 
 ```@docs
-decompress_columns!
+color_groups
 decompress_columns
-decompress_rows!
+decompress_columns!
 decompress_rows
-color_groups
+decompress_rows!
+decompress_symmetric
+decompress_symmetric!
 ```
 
 ## Private
 
+### Matrices
+
+```@docs
+matrix_versions
+respectful_similar
+same_sparsity_pattern
+```
+
 ### Graphs
 
 ```@docs

src/SparseMatrixColorings.jl

@@ -8,7 +8,16 @@ module SparseMatrixColorings
 using ADTypes: ADTypes, AbstractColoringAlgorithm
 using Compat: @compat
 using DocStringExtensions: README
-using LinearAlgebra: Diagonal, Transpose, checksquare, parent, transpose
+using LinearAlgebra:
+    Adjoint,
+    Diagonal,
+    Symmetric,
+    Transpose,
+    adjoint,
+    checksquare,
+    issymmetric,
+    parent,
+    transpose
 using Random: AbstractRNG, default_rng, randperm
 using SparseArrays:
     SparseArrays,
@@ -25,15 +34,18 @@ using SparseArrays:
 include("graph.jl")
 include("order.jl")
 include("coloring.jl")
+include("groups.jl")
 include("adtypes.jl")
-include("check.jl")
+include("matrices.jl")
 include("decompression.jl")
+include("check.jl")
 
 @compat public GreedyColoringAlgorithm
 @compat public NaturalOrder, RandomOrder, LargestFirst
-@compat public decompress_columns!, decompress_columns
-@compat public decompress_rows!, decompress_rows
 @compat public color_groups
+@compat public decompress_columns, decompress_columns!
+@compat public decompress_rows, decompress_rows!
+@compat public decompress_symmetric, decompress_symmetric!
 
 export GreedyColoringAlgorithm
 

src/check.jl

@@ -12,14 +12,15 @@ A partition of the columns of a matrix `A` is _structurally orthogonal_ if, for
     This function is not coded with efficiency in mind, it is designed for small-scale tests.
 """
 function check_structurally_orthogonal_columns(
-    A::AbstractMatrix, colors::AbstractVector{<:Integer}; verbose::Bool=false
+    A::AbstractMatrix, colors::AbstractVector{<:Integer}; verbose::Bool=true
 )
-    for c in unique(colors)
-        js = filter(j -> colors[j] == c, axes(A, 2))
-        Ajs = @view A[:, js]
-        nonzeros_per_row = count(!iszero, Ajs; dims=2)
-        if maximum(nonzeros_per_row) > 1
-            verbose && @warn "Color $c has columns $js sharing nonzeros"
+    groups = color_groups(colors)
+    for (c, g) in enumerate(groups)
+        Ag = @view A[:, g]
+        nonzeros_per_row = dropdims(count(!iszero, Ag; dims=2); dims=2)
+        max_nonzeros_per_row, i = findmax(nonzeros_per_row)
+        if max_nonzeros_per_row > 1
+            verbose && @warn "Columns $g (with color $c) share nonzeros in row $i"
             return false
         end
     end
@@ -40,14 +41,15 @@ A partition of the rows of a matrix `A` is _structurally orthogonal_ if, for eve
     This function is not coded with efficiency in mind, it is designed for small-scale tests.
 """
 function check_structurally_orthogonal_rows(
-    A::AbstractMatrix, colors::AbstractVector{<:Integer}; verbose::Bool=false
+    A::AbstractMatrix, colors::AbstractVector{<:Integer}; verbose::Bool=true
 )
-    for c in unique(colors)
-        is = filter(i -> colors[i] == c, axes(A, 1))
-        Ais = @view A[is, :]
-        nonzeros_per_column = count(!iszero, Ais; dims=1)
-        if maximum(nonzeros_per_column) > 1
-            verbose && @warn "Color $c has rows $is sharing nonzeros"
+    groups = color_groups(colors)
+    for (c, g) in enumerate(groups)
+        Ag = @view A[g, :]
+        nonzeros_per_col = dropdims(count(!iszero, Ag; dims=1); dims=1)
+        max_nonzeros_per_col, j = findmax(nonzeros_per_col)
+        if max_nonzeros_per_col > 1
+            verbose && @warn "Rows $g (with color $c) share nonzeros in column $j"
             return false
         end
     end
@@ -71,24 +73,26 @@ A partition of the columns of a symmetrix matrix `A` is _symmetrically orthogona
     This function is not coded with efficiency in mind, it is designed for small-scale tests.
 """
 function check_symmetrically_orthogonal(
-    A::AbstractMatrix, colors::AbstractVector{<:Integer}; verbose::Bool=false
+    A::AbstractMatrix, colors::AbstractVector{<:Integer}; verbose::Bool=true
 )
+    checksquare(A)
+    issymmetric(A) || return false
+    groups = color_groups(colors)
     for i in axes(A, 2), j in axes(A, 2)
-        if !iszero(A[i, j])
-            group_i = filter(i2 -> (i2 != i) && (colors[i2] == colors[i]), axes(A, 2))
-            group_j = filter(j2 -> (j2 != j) && (colors[j2] == colors[j]), axes(A, 2))
-            A_group_i_column_j = @view A[group_i, j]
-            A_group_j_column_i = @view A[group_j, i]
-            nonzeros_group_i_column_j = count(!iszero, A_group_i_column_j)
-            nonzeros_group_j_column_i = count(!iszero, A_group_j_column_i)
-            if nonzeros_group_i_column_j > 0 && nonzeros_group_j_column_i > 0
-                verbose && @warn """
-                For coefficient $((i, j)), both of the following have confounding zeros:
-                - color $(colors[j]) with group $group_j
-                - color $(colors[i]) with group $group_i
-                """
-                return false
-            end
+        iszero(A[i, j]) && continue
+        ki, kj = colors[i], colors[j]
+        gi, gj = groups[ki], groups[kj]
+        A_gj_rowi = view(A, i, gj)
+        A_gi_rowj = view(A, j, gi)
+        nonzeros_gj_rowi = count(!iszero, A_gj_rowi)
+        nonzeros_gi_rowj = count(!iszero, A_gi_rowj)
+        if nonzeros_gj_rowi > 1 && nonzeros_gi_rowj > 1
+            verbose && @warn """
+            For coefficient $((i, j)):
+            - columns $gj (with color $kj) share nonzeros in row $i
+            - columns $gi (with color $ki) share nonzeros in row $j
+            """
+            return false
         end
     end
     return true

src/decompression.jl

@@ -1,50 +1,3 @@
-transpose_respecting_similar(A::AbstractMatrix, ::Type{T}) where {T} = similar(A, T)
-
-function transpose_respecting_similar(A::Transpose, ::Type{T}) where {T}
-    return transpose(similar(parent(A), T))
-end
-
-function same_sparsity_pattern(A::SparseMatrixCSC, B::SparseMatrixCSC)
-    if size(A) != size(B)
-        return false
-    elseif nnz(A) != nnz(B)
-        return false
-    else
-        for j in axes(A, 2)
-            rA = nzrange(A, j)
-            rB = nzrange(B, j)
-            if rA != rB
-                return false
-            end
-            # TODO: check rowvals?
-        end
-        return true
-    end
-end
-
-function same_sparsity_pattern(
-    A::Transpose{<:Any,<:SparseMatrixCSC}, B::Transpose{<:Any,<:SparseMatrixCSC}
-)
-    return same_sparsity_pattern(parent(A), parent(B))
-end
-
-"""
-    color_groups(colors)
-
-Return `groups::Vector{Vector{Int}}` such that `i ∈ groups[c]` iff `colors[i] == c`.
-
-Assumes the colors are contiguously numbered from `1` to some `cmax`.
-"""
-function color_groups(colors::AbstractVector{<:Integer})
-    cmin, cmax = extrema(colors)
-    @assert cmin == 1
-    groups = [Int[] for c in 1:cmax]
-    for (k, c) in enumerate(colors)
-        push!(groups[c], k)
-    end
-    return groups
-end
-
 ## Column decompression
 
 """
@@ -67,6 +20,9 @@ function decompress_columns!(
     C::AbstractMatrix{R},
     colors::AbstractVector{<:Integer},
 ) where {R<:Real}
+    if !same_sparsity_pattern(A, S)
+        throw(DimensionMismatch("`A` and `S` must have the same sparsity pattern."))
+    end
     A .= zero(R)
     for j in axes(A, 2)
         k = colors[j]
@@ -114,7 +70,7 @@ Here, `colors` is a column coloring of `S`, while `C` is a compressed representa
 function decompress_columns(
     S::AbstractMatrix{Bool}, C::AbstractMatrix{R}, colors::AbstractVector{<:Integer}
 ) where {R<:Real}
-    A = transpose_respecting_similar(S, R)
+    A = respectful_similar(S, R)
     return decompress_columns!(A, S, C, colors)
 end
 
@@ -140,6 +96,9 @@ function decompress_rows!(
     C::AbstractMatrix{R},
     colors::AbstractVector{<:Integer},
 ) where {R<:Real}
+    if !same_sparsity_pattern(A, S)
+        throw(DimensionMismatch("`A` and `S` must have the same sparsity pattern."))
+    end
     A .= zero(R)
     for i in axes(A, 1)
         k = colors[i]
@@ -152,8 +111,8 @@ function decompress_rows!(
 end
 
 function decompress_rows!(
-    A::Transpose{R,<:SparseMatrixCSC{R}},
-    S::Transpose{Bool,<:SparseMatrixCSC{Bool}},
+    A::TransposeOrAdjoint{R,<:SparseMatrixCSC{R}},
+    S::TransposeOrAdjoint{Bool,<:SparseMatrixCSC{Bool}},
     C::AbstractMatrix{R},
     colors::AbstractVector{<:Integer},
 ) where {R<:Real}
@@ -188,6 +147,75 @@ Here, `colors` is a row coloring of `S`, while `C` is a compressed representatio
 function decompress_rows(
     S::AbstractMatrix{Bool}, C::AbstractMatrix{R}, colors::AbstractVector{<:Integer}
 ) where {R<:Real}
-    A = transpose_respecting_similar(S, R)
+    A = respectful_similar(S, R)
     return decompress_rows!(A, S, C, colors)
 end
+
+## Symmetric decompression
+
+"""
+    decompress_symmetric!(
+        A::AbstractMatrix{R},
+        S::AbstractMatrix{Bool},
+        C::AbstractMatrix{R},
+        colors::AbstractVector{<:Integer}
+    ) where {R<:Real}
+
+Decompress the thin matrix `C` into the symmetric matrix `A` which must have the same sparsity pattern as `S`.
+
+Here, `colors` is a symmetric coloring of `S`, while `C` is a compressed representation of matrix `A` obtained by summing the columns that share the same color.
+"""
+function decompress_symmetric! end
+
+function decompress_symmetric!(
+    A::AbstractMatrix{R},
+    S::AbstractMatrix{Bool},
+    C::AbstractMatrix{R},
+    colors::AbstractVector{<:Integer},
+) where {R<:Real}
+    A .= zero(R)
+    groups = color_groups(colors)
+    checksquare(A)
+    for i in axes(A, 1), j in axes(A, 2)
+        iszero(S[i, j]) && continue
+        ki, kj = colors[i], colors[j]
+        gi, gj = groups[ki], groups[kj]
+        if sum(!iszero, view(S, i, gj)) == 1
+            A[i, j] = C[i, kj]
+        elseif sum(!iszero, view(S, j, gi)) == 1
+            A[i, j] = C[j, ki]
+        else
+            error("Symmetric coloring is not valid")
+        end
+    end
+    return A
+end
+
+function decompress_symmetric!(
+    A::Symmetric{R},
+    S::AbstractMatrix{Bool},
+    C::AbstractMatrix{R},
+    colors::AbstractVector{<:Integer},
+) where {R<:Real}
+    # requires parent decompression to handle both upper and lower triangles
+    decompress_symmetric!(parent(A), S, C, colors)
+    return A
+end
+
+"""
+    decompress_symmetric(
+        S::AbstractMatrix{Bool},
+        C::AbstractMatrix{R},
+        colors::AbstractVector{<:Integer}
+    ) where {R<:Real}
+
+Decompress the thin matrix `C` into a new symmetric matrix `A` with the same sparsity pattern as `S`.
+
+Here, `colors` is a symmetric coloring of `S`, while `C` is a compressed representation of matrix `A` obtained by summing the columns that share the same color.
+"""
+function decompress_symmetric(
+    S::AbstractMatrix{Bool}, C::AbstractMatrix{R}, colors::AbstractVector{<:Integer}
+) where {R<:Real}
+    A = respectful_similar(S, R)
+    return decompress_symmetric!(A, S, C, colors)
+end

src/groups.jl

@@ -0,0 +1,16 @@
+"""
+    color_groups(colors)
+
+Return `groups::Vector{Vector{Int}}` such that `i ∈ groups[c]` iff `colors[i] == c`.
+
+Assumes the colors are contiguously numbered from `1` to some `cmax`.
+"""
+function color_groups(colors::AbstractVector{<:Integer})
+    cmin, cmax = extrema(colors)
+    @assert cmin == 1
+    groups = [Int[] for c in 1:cmax]
+    for (k, c) in enumerate(colors)
+        push!(groups[c], k)
+    end
+    return groups
+end