Fix double counting of Symmetric/Hermitian off-diagonals

Project.toml

@@ -1,6 +1,6 @@
 name = "FiniteDifferences"
 uuid = "26cc04aa-876d-5657-8c51-4c34ba976000"
-version = "0.12.34"
+version = "0.12.35"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

src/to_vec.jl

@@ -106,7 +106,11 @@ end
 function to_vec(x::T) where {T<:LinearAlgebra.HermOrSym}
     x_vec, back = to_vec(Matrix(x))
     function HermOrSym_from_vec(x_vec)
-        return T(back(x_vec), x.uplo)
+        # Project onto symmetric/hermitian matrices instead of reading only one triangle:
+        # that would double the off-diagonals of cotangents, see issue #164
+        X = back(x_vec)
+        X_sym = (X .+ (T <: Symmetric ? transpose(X) : adjoint(X))) ./ 2
+        return T(X_sym, x.uplo)
     end
     return x_vec, HermOrSym_from_vec
 end

test/grad.jl

@@ -218,6 +218,24 @@ using FiniteDifferences: grad, jacobian, _jvp, jvp, j′vp, _j′vp, to_vec
     end
 end
 
+# https://github.com/JuliaDiff/FiniteDifferences.jl/issues/164
+@testset "jvp / j′vp(::$Op{$T})" for Op in (Symmetric, Hermitian), T in (Float64, ComplexF64)
+    rng, fdm = MersenneTwister(123456), central_fdm(5, 1)
+    @testset "uplo = $uplo" for uplo in (:U, :L)
+        X = Op(randn(rng, T, 3, 3), uplo)
+
+        # The cotangent of `Matrix(X)` is the projection of the output cotangent onto
+        # symmetric/hermitian matrices, without doubling the off-diagonals
+        Ȳ = randn(rng, T, 3, 3)
+        X̄ = j′vp(fdm, Matrix, Ȳ, X)[1]
+        @test X̄ isa Op
+        @test X̄ ≈ (Ȳ .+ (Op === Symmetric ? transpose(Ȳ) : adjoint(Ȳ))) ./ 2
+
+        Ẋ = Op(randn(rng, T, 3, 3), uplo)
+        @test jvp(fdm, Matrix, (X, Ẋ)) ≈ Ẋ
+    end
+end
+
 using LinearAlgebra
 
 function partial_nan_returning(x)