gdalle / SparseMatrixColorings.jl

Coloring algorithms for sparse Jacobian and Hessian matrices
https://gdalle.github.io/SparseMatrixColorings.jl/
MIT License
17 stars 4 forks source link

Postprocessing for symmetric coloring in bicoloring #122

Open gdalle opened 1 month ago

gdalle commented 1 month ago

We need some way to detect useless colors

amontoison commented 1 week ago

@gdalle

using SparseMatrixColorings

###################################################################
# Rectangle matrix
###################################################################

m = 20
n = 30
p = 2
J = zeros(m,n)
for i = 1:p
  J[i,:] .= rand()
end
for i = m-p+1:m
  J[i,:] .= rand()
end
for j = 1:p
  J[:,j] .= rand()
end
for j = n-p+1:n
  J[:,j] .= rand()
end

J = sparse(J)
H = [0*I J'; J 0*I]

###################################################################
# rectangle -- star bicoloring
###################################################################

problem = ColoringProblem(; structure=:nonsymmetric,
                            partition=:bidirectional)

order = NaturalOrder()

algo = GreedyColoringAlgorithm(order;
                               decompression=:direct)

result = coloring(J, problem, algo)

rcolors = row_colors(result)
ccolors = column_colors(result)
num_colors = ncolors(result)

Br, Bc = compress(J, result)

# H2[i,j] is the color of the column j of H.
#
# Because of the shape of H:
# - If the column j of H represents a column of J (1 ≤ j ≤ n),
# H2[i,j] stores the color of the row `i-n` of J such that J[i-n,j] != 0.
# - If the column j of H represents a row of J (n+1 ≤ j ≤ n+m),
# H2[i,j] stores the color of the column `i` of J such that J[j-n,i] != 0.
H2 = copy(H)
for j in 1:n+m
  for index in H2.colptr[j] : H2.colptr[j+1]-1
    i = H2.rowval[index]
    if i > n
      H2.nzval[index] = rcolors[i-n]
    else
      H2.nzval[index] = ccolors[i]
    end
  end
end
H2 = Int.(H2)

# `checker` is a vector of size `num_colors` that will store, for a given column `j` of H,
# the number of times each color is used in `rcolors` or `ccolors`.
#
# `checker[c]` stores the number of times that color `c` appears in column `j` of H2.
#
# If each non-zero entry in column `j` of H2 has a different color
# (checker[c] <= 1 for all c), it means we can use a neutral color 0.
# This implies that all non-zero entries in the column `j` of H can be recovered from the rows `i` of H where H[i, j] != 0.
checker = Vector{Int}(undef, num_colors)

# List of columns of H that be colored with the neutral colors 0.
neutral_columns = Int[]

# Here, we need to loop through the list of columns that don't have a diagonal entry H[k,k].
# This means all columns for bicoloring, but it could also be used on a subset of columns
# in more traditional star coloring for symmetric matrices!
for j in 1:n+m
  fill!(checker, 0)
  nnz_col = H2.colptr[j+1] - H2.colptr[j]
  neutral = true
  index = 1
  pos = H2.colptr[j] - 1
  while neutral && index ≤ nnz_col
    color = H2.nzval[pos+index]
    if checker[color] != 0
      neutral = false
    else
      checker[color] += 1
      index += 1
    end
  end
  if neutral
    push!(neutral_columns, j)
  end
end

# Use `neutral_columns` to set the neutral color 0 in `rcolors` and `colors`.
for val in neutral_columns
  if val ≤ n
    ccolors[val] = 0
  else
    rcolors[val-n] = 0
  end
end

# Remark: We can't avoid storing `H2` in a more efficient implementation.
amontoison commented 1 week ago

@gdalle I did some plots and it seems to do (almost) what we want: rectangle_star_decompressed_H rectangle_star_pp_decompressed_H

rectangle_star_decompressed_J rectangle_star_pp_decompressed_J

[rectangle] The number of colors before preprocessing is 10.
[rectangle] The number of colors after preprocessing is 8.

arrowhead_star_decompressed_H arrowhead_star_pp_decompressed_H arrowhead_star_decompressed_J arrowhead_star_pp_decompressed_J

[arrowhead] The number of colors before preprocessing is 22.
[arrowhead] The number of colors after preprocessing is 20.

We can see in the last figure that some non-zeros are not colored, probably because we need to consider the previous neutral colors.

If we have a square matrix with an unsmart coloring that assigns a distinct color to each row and each column, the post-processing will set the neutral color for everyone, because each column can be recovered from all rows and vice versa. Each row can be also recovered from all columns.