Fix Hessian coloring: use lower triangular CSC format

- Build lower triangular sparsity pattern (row >= col) instead of full
  symmetric matrix, so decompress_csc! with :L works correctly
- Always add diagonal entries explicitly (old Coloring module did this
  implicitly; SMC requires them to be in the sparsity pattern)
- Get CSC structure directly from the sparse matrix instead of going
  through compress/decompress
- Switch decompress_csc! from :U to :L to match lower triangular storage
- Update test expectations to CSC lower triangular order:
  J=[1,1,2] (was [1,2,1]), V=[3.4,2.1,1.3] (was [3.4,1.3,2.1])

https://claude.ai/code/session_01WBu9hZukriWDSSybN9gfBq

Author: claude

src/coloring.jl

@@ -40,15 +40,12 @@ function _hessian_color_preprocess(
 )
     resize!(seen_idx, num_total_var)
     I, J = Int[], Int[]
-    for (i, j) in edgelist
-        push!(seen_idx, i)
-        push!(seen_idx, j)
-        push!(I, i)
-        push!(J, j)
-        if i != j
-            push!(I, j)
-            push!(J, i)
-        end
+    for (ei, ej) in edgelist
+        push!(seen_idx, ei)
+        push!(seen_idx, ej)
+        # Store in lower triangular format: row >= col
+        push!(I, max(ei, ej))
+        push!(J, min(ei, ej))
     end
     local_indices = sort!(collect(seen_idx))
     empty!(seen_idx)
@@ -63,11 +60,16 @@ function _hessian_color_preprocess(
         J[k] = global_to_local_idx[J[k]]
     end
 
-    # Create sparsity pattern matrix
     n = length(local_indices)
-    S = SMC.SparsityPatternCSC(
-        SparseArrays.sparse(I, J, trues(length(I)), n, n, &),
-    )
+    # Always include diagonal entries (needed for Hessian recovery)
+    for k in 1:n
+        push!(I, k)
+        push!(J, k)
+    end
+
+    # Create lower triangular sparsity pattern (including diagonal)
+    mat = SparseArrays.sparse(I, J, trues(length(I)), n, n, |)
+    S = SMC.SparsityPatternCSC(mat)
 
     # Perform coloring using SMC
     problem = SMC.ColoringProblem(; structure = :symmetric, partition = :column)
@@ -76,12 +78,10 @@ function _hessian_color_preprocess(
     # Wrap result with local_indices
     result = ColoringResult(tree_result, local_indices)
 
-    # SparseMatrixColorings assumes that `I` and `J` are CSC-ordered
-    B = SMC.compress(S, tree_result)
-    C = SMC.decompress(B, tree_result)
-    I_sorted, J_sorted = SparseArrays.findnz(C)
+    # Get CSC-ordered indices directly from the sparse matrix
+    I_sorted, J_sorted, _ = SparseArrays.findnz(mat)
 
-    return C.colptr, I_sorted, J_sorted, result
+    return copy(mat.colptr), I_sorted, J_sorted, result
 end
 
 """
@@ -134,6 +134,6 @@ function _recover_from_matmat!(
     result::ColoringResult,
     stored_values::AbstractVector{T},
 ) where {T}
-    SMC.decompress_csc!(V, colptr, R, result.result, :U)
+    SMC.decompress_csc!(V, colptr, R, result.result, :L)
     return
 end

test/ReverseAD.jl

@@ -426,13 +426,13 @@ function test_coloring_end_to_end_hessian_coloring_and_recovery()
     R = ArrayDiff._seed_matrix(rinfo)
     ArrayDiff._prepare_seed_matrix!(R, rinfo)
     @test I == [1, 2, 2]
-    @test J == [1, 2, 1]
+    @test J == [1, 1, 2]
     @test R == [1.0 0.0; 0.0 1.0]
     hess = [3.4 2.1; 2.1 1.3]
     matmat = hess * R
     V = zeros(3)
     ArrayDiff._recover_from_matmat!(colptr, V, matmat, rinfo, zeros(3))
-    @test V == [3.4, 1.3, 2.1]
+    @test V == [3.4, 2.1, 1.3]
     return
 end