[devmotion]

Fixes #164.

to_vec(::HermOrSym) vectorizes the dense matrix, but its from_vec rebuilt the result from only one triangle. For primals that's fine, because both triangles agree. For cotangents (j′vp, grad) each off-diagonal entry collects contributions from both triangles, so the result had its off-diagonals doubled:

julia> grad(central_fdm(5, 1), sum, Symmetric([1.0 2.0; 3.0 4.0]))[1]
2×2 Symmetric{Float64, Matrix{Float64}}:
 1.0  2.0
 2.0  1.0

from_vec now rebuilds from the symmetric part (X + transpose(X)) / 2, or the hermitian part (X + X') / 2. That gives the natural cotangent, i.e. all ones above. It also matches ChainRulesCore.ProjectTo(::Symmetric), and primals still round-trip exactly.

Behaviour change: j′vp, grad and jacobian with respect to Symmetric/Hermitian inputs return different (correct) values. For example, the raw jacobian(fdm, sum, Symmetric(A)) changes from [1 0 2 1] to [1 1 1 1], with the two triangles now sharing the derivative. ChainRules.jl has rules and tests that compensate for the old double counting, for example rrule(::Type{<:Matrix}, ::HermOrSym). They will fail once this is released, so I'll open a ChainRules.jl PR that switches them to natural tangents and requires FiniteDifferences 0.12.35. Please coordinate the release with that PR. This change came up via JuliaDiff/ChainRules.jl#783.

Tests: new jvp/j′vp testset for Symmetric/Hermitian × Float64/ComplexF64 × both uplo. It fails on main (off-diagonals 2×) and passes here. Bumps the version to 0.12.35.

🤖 Generated with Claude Code

[oxinabox]

I am tossing up, is this the right way or should we be actually slicing out just one triangle.
So that the to_vec has the same dimensionality as the true degrees of freedom.

[y[i,j] for i in axes(y, 1) for j in axes(y,2) if i >= j]

(and flip i<=j for other value of uplo)
and a more complex back.

rather than

x_vec, back = to_vec(Matrix(x))

[devmotion]

That was my first thought as well (see also #164 and #132). The problem is that to_vec returns a single back, which is used for both values and cotangents. With one triangle, the two need different maps.

Example: for a 2×2 Symmetric matrix A, take the lower triangle v = (A₁₁, A₂₁, A₂₂) ∈ ℝ³ as to_vec(A). The corresponding back is the linear map φ: ℝ³ → ℝ²ˣ², φ(v) = [v₁ v₂; v₂ v₃]. Its adjoint w.r.t. the Frobenius inner product ⟨X, Y⟩ = tr(XᵀY) is φᵀ(X) = (X₁₁, X₂₁ + X₁₂, X₂₂). Let G = φᵀφ = diag(1, 2, 1).

For f = sum, f(φ(v)) = v₁ + 2v₂ + v₃. So j′vp with respect to v returns the gradient g = (1, 2, 1), and back(g) = φ(g) = [1 2; 2 1]. That's the same double counting as on main. The natural cotangent is the symmetric Ā with ⟨Ā, φ(w)⟩ = gᵀw for all w ∈ ℝ³, i.e. Ā = φ(G⁻¹g) = [1 1; 1 1]. So values would need φ and cotangents φ ∘ G⁻¹.

ProjectTo doesn't fix this either: φ(g) is already symmetric, so ProjectTo(A)(φ(g)) = φ(g). It only gives the natural cotangent when applied to the gradient w.r.t. the dense matrix, here ProjectTo(A)(ones(2, 2)) = [1 1; 1 1].

That's what this PR does. to_vec uses the dense matrix, and from_vec projects onto the symmetric part, X ↦ (X + Xᵀ)/2 = φ(G⁻¹φᵀ(X)), the same as ProjectTo(::Symmetric). On primals the projection is the identity, so one back is correct for both values and cotangents. It is also consistent with how Diagonal and triangular matrices are handled (#186). Intrinsic coordinates could still be used internally to reduce the number of function evaluations, without changing the output. If J_v is the Jacobian of v ↦ f(φ(v)) (3 columns) and J the Jacobian of X ↦ f((X + Xᵀ)/2) w.r.t. the dense matrix (4 columns), the chain rule gives J = J_v G⁻¹ φᵀ.