`Multifrontal.selinv(A, F)` returns wrong derivatives when `A`'s pattern is not chordal

## Summary

The two-argument differentiable form `Multifrontal.selinv(A, F)` computes a **correct primal**, but both its reverse and forward rules return **wrong derivatives** whenever `A`'s sparsity pattern is a strict subset of the Cholesky factor's fill pattern — equivalently, whenever the graph of `A` is not already chordal.

The error is silent and not small: **≈4.7% on a 5×5 lattice**, and it grows with the amount of fill. The smallest failing case is the 4-cycle (`n = 4`, one fill entry, 0.63% error).

`logdet(A, F)`'s rule is correct, for a reason explained below — so the two rules disagree in a way that makes this easy to miss.

## Reproduction

Julia 1.11.9, CliqueTrees 1.19.4, Mooncake 0.5.40, macOS/aarch64. Fresh environment with only `CliqueTrees` and `Mooncake` added.

The objective is `f(θ) = dot(w, diag(Σ(θ)))` with `Q(θ) = Q₀ + θI` and `Σ = Q⁻¹`. Since `dQ/dθ = I`, we have `dΣ/dθ = -ΣΣ`, so `df/dθ = -dot(w, vec(sum(abs2, Σ, dims = 2)))` exactly — evaluated densely as ground truth.

```julia
using CliqueTrees, LinearAlgebra, SparseArrays, Random, Mooncake
using CliqueTrees.Multifrontal: ChordalCholesky, selinv, selinv!, fisher!, project, selaxpy!
using CliqueTrees.Multifrontal.Differential: selinv_rrule_impl!, selinv_frule_impl, scldia!

const SymSp = Symmetric{Float64, SparseMatrixCSC{Float64, Int}}

# Mooncake cannot trace the BLAS calls inside cholesky!, and ChordalCholesky has
# tangent_type NoTangent anyway, so the factorization is an opaque constant.
# All of the θ-dependence is supposed to be supplied by the selinv rule.
factorize(H) = cholesky!(ChordalCholesky(H))
Mooncake.@zero_adjoint Mooncake.MinimalCtx Tuple{typeof(factorize), SymSp}

# ---------------------------------------------------------------- patterns
# each returns (tril(Q₀), tril(dQ/dθ)) sharing one structural pattern

function path(n)
    return tril(spdiagm(-1 => -ones(n - 1), 0 => fill(2.0, n), 1 => -ones(n - 1))),
        tril(spdiagm(-1 => zeros(n - 1), 0 => ones(n), 1 => zeros(n - 1)))
end

function cycle(k)
    I_, J_, V_, D_ = Int[], Int[], Float64[], Float64[]
    for i in 1:k
        push!(I_, i); push!(J_, i); push!(V_, 4.0); push!(D_, 1.0)
        for j in (mod1(i - 1, k), mod1(i + 1, k))
            push!(I_, i); push!(J_, j); push!(V_, -1.0); push!(D_, 0.0)
        end
    end
    return tril(sparse(I_, J_, V_, k, k)), tril(sparse(I_, J_, D_, k, k))
end

function lattice(m)
    n = m * m; idx(i, j) = (j - 1) * m + i
    I_, J_, V_, D_ = Int[], Int[], Float64[], Float64[]
    for i in 1:m, j in 1:m
        push!(I_, idx(i, j)); push!(J_, idx(i, j)); push!(V_, 4.0); push!(D_, 1.0)
        for (di, dj) in ((1, 0), (-1, 0), (0, 1), (0, -1))
            ii, jj = i + di, j + dj
            if 1 <= ii <= m && 1 <= jj <= m
                push!(I_, idx(i, j)); push!(J_, idx(ii, jj)); push!(V_, -1.0); push!(D_, 0.0)
            end
        end
    end
    return tril(sparse(I_, J_, V_, n, n)), tril(sparse(I_, J_, D_, n, n))
end

function complete(n)
    Random.seed!(7); M = randn(n, n); Q0 = M * M' + n * I
    I_, J_, V_, D_ = Int[], Int[], Float64[], Float64[]
    for j in 1:n, i in j:n
        push!(I_, i); push!(J_, j); push!(V_, Q0[i, j]); push!(D_, i == j ? 1.0 : 0.0)
    end
    return sparse(I_, J_, V_, n, n), sparse(I_, J_, D_, n, n)
end

# ---------------------------------------------------------------- helpers

zeropat(P) = SparseMatrixCSC(P.m, P.n, copy(P.colptr), copy(P.rowval), zeros(nnz(P)))

# structural fill of the Cholesky factor under CliqueTrees' own ordering;
# zero exactly when A's pattern is already chordal
function structural_fill(P, p)
    n = size(P, 1); A = (Matrix(Symmetric(Matrix(P), :L)) .!= 0)[p, p]; f = 0
    for k in 1:n
        nb = [i for i in (k + 1):n if A[i, k]]
        for a in nb, b in nb
            a != b && !A[a, b] && (A[a, b] = true; f += 1)
        end
    end
    return f ÷ 2
end

# the proposed fix: hand fisher! the UNPROJECTED selected inverse
function rrule_fixed!(SA, A, F, DB)
    dA = scldia!(project(A, fisher!(scldia!(copyto!(similar(F), DB), 2), F,
                selinv!(copy(F)); inv = false)), 1 / 2)
    selaxpy!(-1, dA, SA)
    return SA
end

function check(name, L0, dL; theta = 0.7)
    n = size(L0, 1); Random.seed!(42); w = randn(n)
    mkH(t) = Symmetric(SparseMatrixCSC(n, n, L0.colptr, L0.rowval,
            nonzeros(L0) .+ t .* nonzeros(dL)), :L)

    H = mkH(theta)
    S = inv(Matrix(Symmetric(Matrix(parent(H)), :L)))
    exact = -dot(w, vec(sum(abs2, S, dims = 2)))       # d/dθ dot(w, diag Σ)
    exact_ld = tr(S)                                    # d/dθ logdet Q

    F = factorize(H)
    fill = structural_fill(parent(H), collect(F.perm))

    f(t)  = (G = mkH(t); dot(w, diag(selinv(G, factorize(G)))))
    fld(t) = (G = mkH(t); logdet(G, factorize(G)))
    mc(g) = (r = Mooncake.build_rrule(g, theta);
        Mooncake.value_and_gradient!!(r, g, theta)[2][2])

    B = selinv(H, F)
    DBz = zeropat(parent(B)); for i in 1:n; DBz[i, i] = w[i]; end
    DB = Symmetric(DBz, :L)
    SA = Symmetric(zeropat(parent(H)), :L); selinv_rrule_impl!(SA, H, F, B, DB)
    SA2 = Symmetric(zeropat(parent(H)), :L); rrule_fixed!(SA2, H, F, DB)
    _, dB = selinv_frule_impl(H, F,
        Symmetric(SparseMatrixCSC(n, n, L0.colptr, L0.rowval, copy(nonzeros(dL))), :L))

    e(g, x) = abs(g - x) / abs(x)
    println(rpad(name, 15), lpad(n, 4), lpad(fill, 7), "   ",
        rpad(round(exact, digits = 12), 16), rpad(round(mc(f), digits = 12), 16),
        rpad(round(e(mc(f), exact), sigdigits = 2), 10),
        rpad(round(e(sum(diag(SA)), exact), sigdigits = 2), 10),
        rpad(round(e(dot(w, diag(dB)), exact), sigdigits = 2), 10),
        rpad(round(e(sum(diag(SA2)), exact), sigdigits = 2), 10),
        round(e(mc(fld), exact_ld), sigdigits = 2))
    return nothing
end

println("julia ", VERSION, " | CliqueTrees ", pkgversion(CliqueTrees),
    " | Mooncake ", pkgversion(Mooncake), "\n")
println(rpad("pattern", 15), lpad("n", 4), lpad("fill", 7), "   ",
    rpad("exact df/dθ", 16), rpad("Mooncake", 16),
    rpad("mc err", 10), rpad("rrule err", 10), rpad("frule err", 10),
    rpad("FIXED", 10), "logdet err")
println("-" ^ 118)
for (nm, (a, b)) in ("path n=25" => path(25), "complete n=6" => complete(6),
        "triangle" => cycle(3), "4-cycle" => cycle(4),
        "lattice 3x3" => lattice(3), "lattice 5x5" => lattice(5))
    check(nm, a, b)
end
```

### Output

```
julia 1.11.9 | CliqueTrees 1.19.4 | Mooncake 0.5.40

pattern           n   fill   exact df/dθ     Mooncake        mc err    rrule err frule err FIXED     logdet err
----------------------------------------------------------------------------------------------------------------------
path n=25        25      0   1.3251629127    1.3251629127    1.7e-16   3.4e-16   1.3e-15   3.4e-16   1.3e-16
complete n=6      6      0   0.03440002267   0.03440002267   4.0e-16   2.0e-16   0.0       2.0e-16   1.9e-16
triangle          3      0   0.063944975439  0.063944975439  2.2e-16   2.2e-16   2.2e-16   2.2e-16   0.0
4-cycle           4      1   0.106155045798  0.105490515173  0.0063    0.0063    0.0063    1.3e-16   2.3e-16
lattice 3x3       9      5   0.117557009541  0.115443261633  0.018     0.018     0.018     2.4e-16   0.0
lattice 5x5      25     38   0.354031933357  0.337345534952  0.047     0.047     0.047     6.3e-16   1.4e-16
```

`fill` is the structural fill of the Cholesky factor under CliqueTrees' own elimination ordering. The correlation is exact over every case tried:

> **zero fill (chordal pattern) ⟺ derivative exact to machine precision; any fill ⟺ derivative wrong.**

Two things worth noting about the table:

- The `rrule err` column calls `selinv_rrule_impl!` **directly, with no AD framework involved at all**, and matches the Mooncake column to every digit. This is not a Mooncake integration problem.
- The `frule err` column shows `selinv_frule_impl` has the **identical** defect, matching the reverse rule to ~15 digits.

Central differences on the same `f(θ)` (`h = 1e-5`) agree with the exact analytic column to ~10 digits in every row, so the reference values are not in question — e.g. for `lattice 5x5`, exact `0.35403193335739863` vs. central difference `0.35403193341387856`.

The `@zero_adjoint` on `factorize` is not hiding a derivative path. `Mooncake.tangent_type(::Type{<:ChordalCholesky})` is already `NoTangent` upstream, and the primal `selinv(A, F)` reads only `A`'s *pattern*, never its values (`diag(selinv(H1, F1)) == diag(selinv(H2, F1))` is `true` for two matrices sharing a pattern but with different values). By design, then, the entire θ-dependence of the objective has to be reconstructed by the `selinv` rule — which is exactly what the direct, AD-free `rrule err` column measures.

Note that `nnz(F)` is *not* a usable proxy for fill here — it counts padded dense supernodal blocks, so it exceeds `nnz(tril(A))` even for chordal patterns that produce exact gradients (e.g. `path n=25` has `nnz(tril A) = 49`, `nnz(F) = 50`, but zero structural fill and an exact gradient).

## Cause

In `src/Multifrontal.jl/src/Differential.jl/src/selinv.jl`:

```julia
function selinv(A::HermOrSymSparse, F::ChordalCholesky)
    Y = copy(F)
    return project(A, selinv!(Y))       # Σ, then truncated to A's pattern; Y discarded
end

function selinv_rrule_impl!(ΣA, A, F, B, ΔB)
     cB = copyto!(similar(F),  B)       # B is the PROJECTED Σ
    cΔB = copyto!(similar(F), ΔB)
     dA = scldia!(project(A, fisher!(scldia!(cΔB, 2), F, cB; inv=false)), 1 / 2)
    selaxpy!(-1, dA, ΣA)
    return ΣA
end
```

`selinv` computes the selected inverse `Y` on the **factor's** pattern and then discards everything outside `A`'s pattern:

- `project(A, B, P)` (`chordal_triangular.jl:1069`) builds its result by `sympermute`-ing `parent(A)`, so the result carries **`A`'s** pattern, not the factor's.
- `copyto!(::ChordalTriangular, ::SparseMatrixCSC)` (`chordal_triangular.jl:852`) begins with `zerorec!(A.Dval); zerorec!(A.Lval)` and then writes only the entries present in the source.

So `cB` is a factor-shaped container in which **every Σ entry on the factor's fill pattern but outside `A`'s pattern is exactly zero**. But

```
∂Σ_ij/∂A_kl = -Σ_ik Σ_lj
```

genuinely depends on those entries — `Σ` is dense, and the fill positions carry real values. `fisher!` is therefore fed a truncated `Σ` and returns a correspondingly truncated (too small) derivative. `selinv_frule_impl` does `cB = copyto!(similar(F), B)` with the same already-projected `B` and is wrong in exactly the same way.

## Why `logdet(A, F)` is unaffected

`logdet`'s rule (sibling `logdet.jl`) only ever needs

```
∂logdet(A)/∂A_kl = Σ_lk
```

restricted to the positions of `A` it is contracted against — so projecting Σ onto `A`'s pattern loses nothing there. The last column of the table confirms `logdet` is exact to machine precision on the very same fill-in matrices. This asymmetry is why the bug survives: any test suite that checks `logpdf`/`logdet` gradients will pass while marginal-variance gradients are quietly wrong.

## Suggested fix

Feed `fisher!` the **unprojected** selected inverse and project only at the end. Substituting `selinv!(copy(F))` for `cB` makes every failing case exact to machine precision — that is the verified `FIXED` column above.

Minimally:

```diff
 function selinv_rrule_impl!(ΣA, A, F, B, ΔB)
-     cB = copyto!(similar(F),  B)
+     cB = selinv!(copy(F))
     cΔB = copyto!(similar(F), ΔB)
      dA = scldia!(project(A, fisher!(scldia!(cΔB, 2), F, cB; inv=false)), 1 / 2)
```

and the same substitution in `selinv_frule_impl`.

That recomputes the selected inversion, which is wasteful, since `selinv(A, F)` has already computed exactly this object and thrown it away. The zero-extra-cost version is to thread the factor-shaped `Y` through instead of the projected `B` — e.g. in `ext/MooncakeExt/selinv.jl`:

```julia
function Mooncake.rrule!!(::CoDual{typeof(selinv)}, cdA, cdF)
    A, dA = primaltangent(cdA)
    Y = selinv!(copy(primal(cdF)))     # keep the unprojected Σ
    B = project(A, Y)                  # what selinv(A, F) returns
    dB = zero(B)
    function pullback!!(::NoRData)
        selinv_rrule_impl!(dA, A, primal(cdF), Y, dB)   # pass Y, not B
        return NoRData(), NoRData(), NoRData()
    end
    return CoDual(B, tofdata(B, dB)), pullback!!
end
```

with `selinv_rrule_impl!` taking the factor-shaped `Y` directly and dropping its `copyto!(similar(F), B)` line. `ΔB` still needs the `copyto!(similar(F), ΔB)` treatment — that one is genuinely a cotangent living on `A`'s pattern, so zero-filling it is correct.

## Downstream impact

This silently corrupts marginal-variance hyperparameter gradients for every [GaussianMarkovRandomFields.jl](https://github.com/timweiland/GaussianMarkovRandomFields.jl) user on the Mooncake path whose precision matrix has Cholesky fill-in — which is most non-trivial spatial models. `logpdf` gradients are unaffected, so the failure does not show up in likelihood-only workflows.


Navigation Menu

Uh oh!

`Multifrontal.selinv(A, F)` returns wrong derivatives when `A`'s pattern is not chordal #19

Summary

Reproduction

Output

Cause

Why `logdet(A, F)` is unaffected

Suggested fix

Downstream impact

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Uh oh!

Multifrontal.selinv(A, F) returns wrong derivatives when A's pattern is not chordal #19

Description

Summary

Reproduction

Output

Cause

Why logdet(A, F) is unaffected

Suggested fix

Downstream impact

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions

`Multifrontal.selinv(A, F)` returns wrong derivatives when `A`'s pattern is not chordal #19

Why `logdet(A, F)` is unaffected