From which HVP do we determine the off-diagonal entries ?
For the edges 34 and 64, respectively between the yellow nodes 3, 6 and the blue nodes 4, they are part of the star centered at the blue node 4 so they are obtained from the entries 3 and 6 blue HVP.
For the edges 12, 32, 52 and 62, respectively between the yellow nodes 1, 3, 5, 6 and the pink node 2, they are part of the star centered at the pink node 2 so they are obtained from the entries 1, 3, 5 and 6 of the pink HVP.
Acyclic coloring
Definition 6.3 A mapping $\phi: V \to \{1, \ldots, p\}$ is an acyclic coloring if
$\phi$ is a distance-1 coloring
every cycle uses at least 3 colors
Theorem 4.6 Let $A$ be a symmetric matrix with nonzero diagonal elements, $G$ be its adjacency matrix. A mapping $\phi$ is an acyclic coloring of the adjacency graph iff it induces a substitutable partition of the columns of $A$.
Name acyclic coloring comes from the fact that the subgraph induced by any pair of colors is a forest.
How is the subgraph induced by a pair of colors ?
It does not contain cycles so it is a union of disjoint trees or in other words a forest.
The concept of tree generalizes the concept of stars. Indeed, a star is a tree and, when rooted at the center of the star, the depth of the tree is 1. Since the depth is one, every edge is incident to a leaf so no substitutions are needed. With acyclic coloring, the edges incident to a leaf will be obtained directly but the edges that are not incident to any tree leaves will need the entries corresponding to the edges deeper in their tree to be determined first.
Sparse AD
Benoît Legat
$(PlutoTeachingTools.ChooseDisplayMode())
$(PlutoUI.TableOfContents(depth=1))
""")metadatashow_logsèdisabled®skip_as_script«code_folded$76162106-fb90-4883-b2d2-b87e1419c072cell_id$76162106-fb90-4883-b2d2-b87e1419c072code,md"""
## Formulation as coloring problem
"""metadatashow_logsèdisabled®skip_as_script«code_folded$d23e7238-5edb-42d3-96d9-0dffa9e7f326cell_id$d23e7238-5edb-42d3-96d9-0dffa9e7f326codeٲhbox([
md"scale = $(@bind(scale, Slider(1:15, default = 10)))",
md"pad = $(@bind(pad, Slider(1:10, default = 3)))",
md"border = $(@bind(border, Slider(1:5, default = 2)))",
])metadatashow_logsèdisabled®skip_as_script«code_folded$a7b85908-404b-475b-9875-c1190de5ff71cell_id$a7b85908-404b-475b-9875-c1190de5ff71code}md"""
[**Definition 4.2**](https://epubs.siam.org/doi/10.1137/S0036144504444711) A coloring of the columns of a symmetric matrix ``A`` is *symmetrically orthogonal* if for every nonzero element ``a_{uv}``, either (1) there is no other column of color ``\phi(u)`` with a nonzero in row ``r_u``, or (2) there is no other column of color ``\phi(v)`` with a nonzero in row ``r_v``.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$15cf368c-cada-4a72-b59b-a0f9dd0115b4cell_id$15cf368c-cada-4a72-b59b-a0f9dd0115b4codeqa(
md"How to compute the non-bold entries ?",
md"""
For star graphs, the entry ``a_{uv}`` was the edge of a star with either center ``u`` and leaf ``v`` or center ``v`` and leaf ``u``.
It was then obtained from the HVP of the color of the center.
We need to generalize this for arbitrary trees.
For any edge from a leaf ``u`` to a node ``v``, we can compute
``a_{uv}`` from the ``u``th entry of the HVP of color ``\phi(v)``.
Then, we subtract ``a_{uv}`` from the ``v``th row of the HVP of color ``\phi(u)``.
So we can do
* ``a_{71} \gets (h_\text{red})_7`` then ``(h_\text{blue})_1 \gets (h_\text{blue})_1 - a_{71}``
* ``a_{43} \gets (h_\text{red})_4`` then ``(h_\text{blue})_3 \gets (h_\text{blue})_3 - a_{43}``
* ``a_{52} \gets (h_\text{blue})_5`` then ``(h_\text{red})_2 \gets (h_\text{red})_2 - a_{52}``
Then, we can remove these leaves and edges. Doing so, some of their parents in the tree become leaves and we can apply the same procedure recursively. So we get
* ``a_{12} \gets (h_\text{blue})_5`` then ``(h_\text{red})_2 \gets (h_\text{red})_2 - a_{12}``
* ``a_{32} \gets (h_\text{blue})_5`` then ``(h_\text{red})_2 \gets (h_\text{red})_2 - a_{32}``
""",
)metadatashow_logsèdisabled®skip_as_script«code_folded$45f4edeb-d934-413f-9f45-568465111b6ecell_id$45f4edeb-d934-413f-9f45-568465111b6ecode/md"## Sparsity detection with multiple outputs"metadatashow_logsèdisabled®skip_as_script«code_folded$b18a3bf6-102e-49f5-af8c-f1c35d2e2a8bcell_id$b18a3bf6-102e-49f5-af8c-f1c35d2e2a8bcodeemd"The sign function is constant over ``x > 0`` and ``x < 0`` so its output carries **no** variables"metadatashow_logsèdisabled®skip_as_script«code_folded$e2eb1894-06ad-47a0-9ee0-c65b24a669b5cell_id$e2eb1894-06ad-47a0-9ee0-c65b24a669b5code\md"""
* [What color is your Jacobian ? Graph Coloring for Computing Derivatives](https://epubs.siam.org/doi/10.1137/S0036144504444711)
* [An Illustrated Guide to Automatic Sparse Differentiation](https://iclr-blogposts.github.io/2025/blog/sparse-autodiff/)
* [Revisiting Sparse Matrix Coloring and Bicoloring](https://arxiv.org/abs/2505.07308)
"""metadatashow_logsèdisabled®skip_as_script«code_folded$44a6b427-bd4b-4fa1-95d8-a8ea55719192cell_id$44a6b427-bd4b-4fa1-95d8-a8ea55719192codeqa(md"How can the diagonal entries be determined ?", md"""
The value of ``(h_{\phi(i)})_i`` is the sum of the entries of the ``i``th row of ``B_{\phi(i)}``.
Since the only nonzero entry of this row is ``A_{ii}``, we have ``A_{ii} = (h_{\phi(i)})_i``.
""")metadatashow_logsèdisabled®skip_as_script«code_folded$dfe3dffd-926b-4083-9e66-8536e5df198dcell_id$dfe3dffd-926b-4083-9e66-8536e5df198dcodebegin
function qa(question, answer)
return @htl("$question
$answer` if `question.content[]` is `Markdown.Paragraph`
# This will print the question on a new line and we don't want that:
h = HTML(_inline_html(question.content[]))
return qa(h, answer)
end
endmetadatashow_logsèdisabled®skip_as_script«code_folded$874b6593-7310-4788-beb0-2f18c0ca09a7cell_id$874b6593-7310-4788-beb0-2f18c0ca09a7code$md"## Three columns in just one JVP"metadatashow_logsèdisabled®skip_as_script«code_folded$9e24f674-42af-4fa8-bf3e-42e14d855a58cell_id$9e24f674-42af-4fa8-bf3e-42e14d855a58codetwo_columns(
md"""
* The columns 1, 2 and 5 do not share any rows with nonzero entries.
* So their entries can be recovered unambiguously with just one JVP!
""",
img("https://iclr-blogposts.github.io/2025/assets/img/2025-04-28-sparse-autodiff/sparse_ad.svg", :width => 400),
)metadatashow_logsèdisabled®skip_as_script«code_folded$6bd04f2f-ed5f-4ecf-91dd-8868ec5b83e1cell_id$6bd04f2f-ed5f-4ecf-91dd-8868ec5b83e1code1using SparseArrays, Images, SparseMatrixColoringsmetadatashow_logsèdisabled®skip_as_script«code_folded$29d9f355-26aa-4714-9f9e-3d020016662ecell_id$29d9f355-26aa-4714-9f9e-3d020016662ecodemd"## Small Example"metadatashow_logsèdisabled®skip_as_script«code_folded$e3dec4ca-69c0-4955-9303-5f8a239546f2cell_id$e3dec4ca-69c0-4955-9303-5f8a239546f2code_md"Let ``h_c`` be the sum of the columns of ``B_c`` hence the result of the corresponding HVP."metadatashow_logsèdisabled®skip_as_script«code_folded$5a6aec48-54fb-4d88-a8a1-b06a2bb99108cell_id$5a6aec48-54fb-4d88-a8a1-b06a2bb99108codeOS = sparse([
0 0 1 1 0 1
1 0 0 0 1 0
0 1 0 0 1 0
0 1 1 0 0 0
])metadatashow_logsèdisabled®skip_as_script«code_folded$75597493-3e69-437f-9408-b43a89b35559cell_id$75597493-3e69-437f-9408-b43a89b35559codeOusing PlutoUI, PlutoUI.ExperimentalLayout, HypertextLiteral, PlutoTeachingToolsmetadatashow_logsèdisabled®skip_as_script«code_folded$1efb0f40-fb9e-4e5b-be59-84910afbf376cell_id$1efb0f40-fb9e-4e5b-be59-84910afbf376codemd"""
## Comparison of chromatic numbers
[Theorem 7.1](https://epubs.siam.org/doi/10.1137/S0036144504444711) For every graph ``G``,
```math
\xi_1(G) \le \xi_\text{acyclic}(G) \le \xi_\text{star}(G) \le \xi_2(G) = \xi_1(G^2)
```
"""metadatashow_logsèdisabled®skip_as_script«code_folded$e882ebee-36ad-43d9-8ab7-20f7e8b2d9fecell_id$e882ebee-36ad-43d9-8ab7-20f7e8b2d9fecodeUmd"""
## Example
Consider the following sparsity pattern of the Jacobian matrix:
"""metadatashow_logsèdisabled®skip_as_script«code_folded$cabeb42c-f40b-4549-b9a2-3bda7a988672cell_id$cabeb42c-f40b-4549-b9a2-3bda7a988672codeimg("Gebremedhin_Figure_6_1")metadatashow_logsèdisabled®skip_as_script«code_folded$e9135060-75d4-4f5a-9315-2e5fbb493cdfcell_id$e9135060-75d4-4f5a-9315-2e5fbb493cdfcodemfunction colored_plots(A; decompression, structure, partition = :column, kws...)
problem = ColoringProblem(; structure, partition)
algo = GreedyColoringAlgorithm(; decompression)
result = coloring(A, problem, algo)
background_color = RGBA(0, 0, 0, 0)
border_color = RGB(0, 0, 0)
colorscheme = distinguishable_colors(
ncolors(result),
[convert(RGB, background_color), convert(RGB, border_color)];
dropseed=true,
)
A_img, B_img = SparseMatrixColorings.show_colors(result; colorscheme, scale, pad, border, background_color, border_color)
if structure == :symmetric
S = A
else
if partition == :column
S = A' * A
else
S = A * A'
end
end
adj = SparseMatrixColorings.AdjacencyGraph(A)
G = Graphs.SimpleDiGraph(S - Diagonal(diag(S)))
gp = gplot(G; nodefillc = colorscheme[result.color], nodelabel=1:Graphs.nv(G), kws...)
A_img, gp, B_img
endmetadatashow_logsèdisabled®skip_as_script«code_folded$d46dcdb3-ac4c-4ba9-ac48-6d6a693630d3cell_id$d46dcdb3-ac4c-4ba9-ac48-6d6a693630d3code$f(x) = x[1] + x[2]*x[3] * sign(x[4])metadatashow_logsèdisabled®skip_as_script«code_folded$a546edb2-6a72-4c5f-95f7-cc7709bbb505cell_id$a546edb2-6a72-4c5f-95f7-cc7709bbb505codeF(xtracer)metadatashow_logsèdisabled®skip_as_script«code_folded$4ab46ef8-0f6e-4e6d-8f6f-ac466977b476cell_id$4ab46ef8-0f6e-4e6d-8f6f-ac466977b476codejmd"Do we need to aim for distance-2 coloring or does the symmetry of the matrix allows us to relax this ?"metadatashow_logsèdisabled®skip_as_script«code_folded$18594f1f-6036-4a95-8682-cf415c998311cell_id$18594f1f-6036-4a95-8682-cf415c998311code$md"# Sparsity detection of Jacobian"metadatashow_logsèdisabled®skip_as_script«code_folded$2fe7f4eb-ca33-4338-ba8e-29b769b0dccecell_id$2fe7f4eb-ca33-4338-ba8e-29b769b0dccecodeElarge = sparse(Symmetric(sprand(StableRNG(0), Bool, n, n, 0.2) + I));metadatashow_logsèdisabled®skip_as_script«code_folded$97c94c5a-fab5-4226-b803-10c85a298f2acell_id$97c94c5a-fab5-4226-b803-10c85a298f2acode2struct MyGradientTracer
indexset::Set{Int}
endmetadatashow_logsèdisabled®skip_as_script«code_folded$60070a94-d9ca-474b-9c6f-37468e3f030dcell_id$60070a94-d9ca-474b-9c6f-37468e3f030dcodeqa(md"How do we get the value of ``a_4`` ?",
md"""
Two ways:
1) Get the value of ``a_2`` from ``(AD)_{12}`` then get ``a_4`` from ``(AD)_{21} - a_2``
2) Get the value of ``a_6`` from ``(AD)_{41}`` then get ``a_4`` from ``(AD)_{32} - a_6``
``a_4`` cannot be read directly from a single entry, it needs **substitutions**. This example shows that enabling substitutions allows us to use less colors to find all entries of the Hessian!
""")metadatashow_logsèdisabled®skip_as_script«code_folded$8e2045f3-ee09-4439-b504-d2ca9f58a3ddcell_id$8e2045f3-ee09-4439-b504-d2ca9f58a3ddcodecBase.:+(a::MyGradientTracer, b::MyGradientTracer) = MyGradientTracer(union(a.indexset, b.indexset))metadatashow_logsèdisabled®skip_as_script«code_folded$f2418632-3261-4acb-8cc6-df333b5480a5cell_id$f2418632-3261-4acb-8cc6-df333b5480a5codemd"## Illustrative example"metadatashow_logsèdisabled®skip_as_script«code_folded$592510fd-a0e4-42e5-ae5b-cac301e74f0acell_id$592510fd-a0e4-42e5-ae5b-cac301e74f0acodeEqa(md"How is the subgraph induced by a pair of colors ?", md"""
It does not contain paths of more than 3 vertices and adjacent vertices have opposite colors.
So it is a union of disjoint connected components and each component is a stars with one root node of a given color connected with leaf nodes of the other color.
""")metadatashow_logsèdisabled®skip_as_script«code_folded$5892dfbb-b375-407b-a0a7-da66adb71b67cell_id$5892dfbb-b375-407b-a0a7-da66adb71b67code2md"""
[**Theorem 3.5**](https://epubs.siam.org/doi/10.1137/S0036144504444711) A coloring of the columns is distance-2 in the bipartite graph iff columns of the same color are structurally orthogonal.
Let the ``k``-th power ``G^k`` of a graph ``G`` be the graph with same nodes as ``G`` and an edge in ``G^k`` for each distance ``k`` neighbors in ``G``.
[**Lemma 3.7**](https://epubs.siam.org/doi/10.1137/S0036144504444711) The column intersection graph is the square of the biparted graph restricted to its column vertices.
[**Lemma 2.1**](https://epubs.siam.org/doi/10.1137/S0036144504444711) A coloring is distance-``k`` in ``G`` iff it is distance-1 in ``G^k``.
**Corollary** A coloring of the columns is distance-1 in the column intersection graph iff columns of the same color are structurally orthogonal.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$726208c6-e856-42a2-b027-555388bc505ecell_id$726208c6-e856-42a2-b027-555388bc505ecodetwo_columns(
img("https://iclr-blogposts.github.io/2025/assets/img/2025-04-28-sparse-autodiff/sparse_ad_reverse_full.svg"),
img("https://iclr-blogposts.github.io/2025/assets/img/2025-04-28-sparse-autodiff/sparse_ad_reverse_decompression.svg"),
)metadatashow_logsèdisabled®skip_as_script«code_folded$12ffe7c2-2bf6-4949-ae28-f3d8e03a5da5cell_id$12ffe7c2-2bf6-4949-ae28-f3d8e03a5da5codeZtwo_columns(
md"Color's columns are structurally orthogonal",
md"Jacobian from 3 JVP",
)metadatashow_logsèdisabled®skip_as_script«code_folded$9df9c76a-48af-4ca6-83f1-b8798d6d4e34cell_id$9df9c76a-48af-4ca6-83f1-b8798d6d4e34code+Base.:/(a::MyGradientTracer, b::Number) = ametadatashow_logsèdisabled®skip_as_script«code_folded$760ef943-e9d9-4cdc-8977-932cf6fcb78acell_id$760ef943-e9d9-4cdc-8977-932cf6fcb78acodemd"## Coloring of a star"metadatashow_logsèdisabled®skip_as_script«code_folded$a01688ec-8f16-4d39-b39d-6d1acd24d7adcell_id$a01688ec-8f16-4d39-b39d-6d1acd24d7adcodetwo_columns(
img("https://iclr-blogposts.github.io/2025/assets/img/2025-04-28-sparse-autodiff/sparse_ad_forward_full.svg", :width => 300),
img("https://iclr-blogposts.github.io/2025/assets/img/2025-04-28-sparse-autodiff/sparse_ad_forward_decompression.svg", :width => 300),
)metadatashow_logsèdisabled®skip_as_script«code_folded$a3bf5961-323b-4adc-99a0-9a456a0bfd48cell_id$a3bf5961-323b-4adc-99a0-9a456a0bfd48code2F(x) = [x[1]*x[2]+sign(x[3]), sign(x[3]) * x[4]/2]metadatashow_logsèdisabled®skip_as_script«code_folded$70296870-49af-4564-bdb7-bd2d49f9114bcell_id$70296870-49af-4564-bdb7-bd2d49f9114bcodeRviz(star, decompression = :direct, structure = :symmetric, plot_size = (5cm, 5cm))metadatashow_logsèdisabled®skip_as_script«code_folded$dcf92bc9-5e0f-4716-bb07-ca0d9b73b166cell_id$dcf92bc9-5e0f-4716-bb07-ca0d9b73b166codeGviz(args...; kws...) = three_columns(colored_plots(args...; kws...)...)metadatashow_logsèdisabled®skip_as_script«code_folded$9b427766-5290-4dee-8dd5-5aab20080d0acell_id$9b427766-5290-4dee-8dd5-5aab20080d0acodeYviz(large, decompression = :substitution, structure = :symmetric, plot_size = (6cm, 6cm))metadatashow_logsèdisabled®skip_as_script«code_folded$4f8de444-dcb4-411a-b087-6948699614e0cell_id$4f8de444-dcb4-411a-b087-6948699614e0code2using Graphs, GraphPlot, LinearAlgebra, StableRNGsmetadatashow_logsèdisabled®skip_as_script«code_folded$e820e135-2df2-4fa3-a09b-8ce7e5a1f745cell_id$e820e135-2df2-4fa3-a09b-8ce7e5a1f745codedtwo_columns(left, right) = hbox([
left,
Div(html" ", style = Dict("flex-grow" => "1")),
right,
])metadatashow_logsèdisabled®skip_as_script«code_folded$26191290-2b2d-43fe-b2bf-d690fd0bbe15cell_id$26191290-2b2d-43fe-b2bf-d690fd0bbe15codeqa(md"From which HVP do we determine the off-diagonal entries ?",
md"""
* For the edges 34 and 64, respectively between the yellow nodes 3, 6 and the blue nodes 4, they are part of the star centered at the blue node 4 so they are obtained from the entries 3 and 6 blue HVP.
* For the edges 12, 32, 52 and 62, respectively between the yellow nodes 1, 3, 5, 6 and the pink node 2, they are part of the star centered at the pink node 2 so they are obtained from the entries 1, 3, 5 and 6 of the pink HVP.
""")metadatashow_logsèdisabled®skip_as_script«code_folded$e3ae5697-6b0b-4d5f-99a0-d579b6eff655cell_id$e3ae5697-6b0b-4d5f-99a0-d579b6eff655codemd"""## Acyclic coloring
[Definition 6.3](https://epubs.siam.org/doi/10.1137/S0036144504444711) A mapping ``\phi: V \to \{1, \ldots, p\}`` is an *acyclic coloring* if
1. ``\phi`` is a distance-1 coloring
2. every **cycle** uses at least 3 colors
[Theorem 4.6](https://epubs.siam.org/doi/10.1137/S0036144504444711)
Let ``A`` be a symmetric matrix with **nonzero diagonal elements**,
``G`` be its adjacency matrix. A mapping ``\phi`` is an acyclic coloring of the adjacency graph iff it induces a substitutable partition of the columns of ``A``.
Name acyclic coloring comes from the fact that the subgraph induced by any pair of colors is a forest.
"""metadatashow_logsèdisabled®skip_as_script«code_folded$95d7e9cd-b4a1-4205-b7d4-8c10916be8a2cell_id$95d7e9cd-b4a1-4205-b7d4-8c10916be8a2codeٱstar = let
I = collect(1:star_n)
J = ones(Int, star_n)
for i in 2:star_n
push!(I, 1)
push!(J, i)
push!(I, i)
push!(J, i)
end
sparse(I, J, ones(Int, length(I)))
end;metadatashow_logsèdisabled®skip_as_script«code_folded$04ea8b79-63ba-4eb4-938c-98a21668924acell_id$04ea8b79-63ba-4eb4-938c-98a21668924acoderqa(md"How is the subgraph induced by a pair of colors ?", md"""
It does not contain cycles so it is a union of disjoint trees or in other words a forest.
The concept of tree generalizes the concept of stars. Indeed, a star is a tree and, when rooted at the center of the star, the depth of the tree is 1.
Since the depth is one, every edge is incident to a leaf so no substitutions are needed. With acyclic coloring, the edges incident to a leaf will be obtained directly but the edges that are not incident to any tree leaves will need the entries corresponding to the edges deeper in their tree to be determined first.
""")metadatashow_logsèdisabled®skip_as_script«code_folded$dbef0e7d-b731-4d83-98e3-104b23f26042cell_id$dbef0e7d-b731-4d83-98e3-104b23f26042codeCmd"`n` = $(@bind n Slider(10:20, show_value = true, default = 16))"metadatashow_logsèdisabled®skip_as_script«code_foldedënotebook_id$5c8963ce-9fcd-11f1-b08d-e5d05c640adcin_temp_dir¨metadata
mimetext/htmlrootassigneelast_run_timestampAڣNE$persist_js_state·has_pluto_hook_features§cell_id$22c1030b-bef5-4af5-9d5f-6afb4a7ae699depends_on_disabled_cells§runtime 2published_object_keysdepends_on_skipped_cells§errored$dc23180e-2b0f-41ae-aebf-0159a8f8677aqueued¤logsrunning¦outputbodyattributesstyle#display: flex; flex-direction: row;tagdivchildrenG
mimetext/htmlrootassigneelast_run_timestampAڣMpersist_js_state·has_pluto_hook_features§cell_id$cabeb42c-f40b-4549-b9a2-3bda7a988672depends_on_disabled_cells§runtime \~published_object_keysdepends_on_skipped_cells§errored$e9135060-75d4-4f5a-9315-2e5fbb493cdfqueued¤logsrunning¦outputbody.colored_plots (generic function with 1 method)mimetext/plainrootassigneelast_run_timestampAڣyBpersist_js_state·has_pluto_hook_features§cell_id$e9135060-75d4-4f5a-9315-2e5fbb493cdfdepends_on_disabled_cells§runtime ;
