Fix grad symmetry and release metadata gaps

.github/workflows/release.yml

@@ -61,6 +61,10 @@ jobs:
 
   build-linux:
     runs-on: ubuntu-latest
+    strategy:
+      fail-fast: false
+      matrix:
+        python-version: ["3.10", "3.11", "3.12"]
     env:
       ELLPHI_USE_EIGEN: "1"
       ELLPHI_EIGEN_INCLUDE: "/usr/include/eigen3"
@@ -70,7 +74,7 @@ jobs:
     - name: Set up Python
       uses: actions/setup-python@v5
       with:
-        python-version: "3.12"
+        python-version: ${{ matrix.python-version }}
 
     - name: Install Eigen headers
       run: sudo apt-get update && sudo apt-get install -y libeigen3-dev
@@ -84,7 +88,7 @@ jobs:
     - name: Upload wheel
       uses: actions/upload-artifact@v4
       with:
-        name: dist-linux
+        name: dist-linux-py${{ matrix.python-version }}
         path: dist/*.whl
 
   build-windows:

src/ellphi/_tangency_cpp_impl.cpp

@@ -27,8 +27,19 @@
 #endif
 #endif
 
-extern "C" const char* tangency_backend_version() { return TANGENCY_VERSION; }
-extern "C" const char* tangency_linalg_kind() { return TANGENCY_LINALG_KIND; }
+#if defined(_WIN32) || defined(__CYGWIN__)
+#define ELLPHI_EXPORT __declspec(dllexport)
+#else
+#define ELLPHI_EXPORT
+#endif
+
+ELLPHI_EXPORT extern "C" const char* tangency_backend_version() {
+    return TANGENCY_VERSION;
+}
+
+ELLPHI_EXPORT extern "C" const char* tangency_linalg_kind() {
+    return TANGENCY_LINALG_KIND;
+}
 
 namespace {
 
@@ -1100,12 +1111,6 @@ void copy_error(char* buffer, std::size_t size, const std::string& message) {
 
 }  // namespace
 
-#if defined(_WIN32) || defined(__CYGWIN__)
-#define ELLPHI_EXPORT __declspec(dllexport)
-#else
-#define ELLPHI_EXPORT
-#endif
-
 ELLPHI_EXPORT extern "C" int tangency_solve(
     const double* pcoef,
     const double* qcoef,

src/ellphi/grad.py

@@ -248,6 +248,11 @@ def vjp(
         # dL/dA[i,j] += grad_c * center[i] * center[j]
         grad_A += grad_c[:, None, None] * np.einsum("ni,nj->nij", centers, centers)
 
+        # A is parameterized by symmetric covariances, so project the
+        # unconstrained matrix gradient back onto the symmetric subspace before
+        # exposing grad_cov or differentiating through inv(cov).
+        grad_A = 0.5 * (grad_A + np.swapaxes(grad_A, -2, -1))
+
         # --- grad w.r.t. centers ---
         # From b = -A @ center: grad_center += -A^T @ grad_b
         grad_centers = -np.einsum("nji,nj->ni", A, grad_b)
@@ -260,13 +265,11 @@ def vjp(
         # dA = -inv(cov) @ dCov @ inv(cov) / s^2
         # <grad_A, dA> = tr(grad_A^T (-inv(cov) dCov inv(cov) / s^2))
         #              = -1/s^2 * tr(inv(cov)^T grad_A^T inv(cov)^T dCov)
-        # Since inv(cov) is symmetric:
-        #   grad_cov = -inv(cov) @ grad_A^T @ inv(cov) / s^2
-        #            = -s^2 A @ grad_A^T @ s^2 A / s^2
-        #            = -s^2 * A @ grad_A^T @ A
-        # A @ grad_A^T @ A: note grad_A^T[i,j] = grad_A[j,i]
-        grad_A_T = np.swapaxes(grad_A, -2, -1)
-        grad_cov = -s2 * np.einsum("nij,njk,nkl->nil", A, grad_A_T, A)
+        # Since inv(cov) and grad_A are symmetric:
+        #   grad_cov = -inv(cov) @ grad_A @ inv(cov) / s^2
+        #            = -s^2 A @ grad_A @ s^2 A / s^2
+        #            = -s^2 * A @ grad_A @ A
+        grad_cov = -s2 * np.einsum("nij,njk,nkl->nil", A, grad_A, A)
 
         return grad_centers, grad_cov
 

tests/test_coef_from_cov_grad.py

@@ -39,7 +39,8 @@ def _numerical_jacobian_cov(centers, covs, scale, h=1e-6):
 
     Perturbs each entry of the symmetric covariance matrix while maintaining
     symmetry. For the symmetric parameterization, we perturb both (j,k) and
-    (k,j) together and assign the derivative to both entries.
+    (k,j) together and split the off-diagonal derivative evenly across the
+    returned matrix entries.
     """
     n, d = centers.shape
     coefs0 = coef_from_cov(centers, covs, scale=scale)
@@ -50,15 +51,20 @@ def _numerical_jacobian_cov(centers, covs, scale, h=1e-6):
             for k in range(j, d):
                 cov_plus = covs.copy()
                 cov_plus[i, j, k] += h
-                cov_plus[i, k, j] += h  # keep symmetric (no-op when j==k)
+                if j != k:
+                    cov_plus[i, k, j] += h
                 cov_minus = covs.copy()
                 cov_minus[i, j, k] -= h
-                cov_minus[i, k, j] -= h
+                if j != k:
+                    cov_minus[i, k, j] -= h
                 c_plus = coef_from_cov(centers, cov_plus, scale=scale)
                 c_minus = coef_from_cov(centers, cov_minus, scale=scale)
                 deriv = (c_plus - c_minus) / (2 * h)
-                jac[:, :, i, j, k] = deriv
-                jac[:, :, i, k, j] = deriv
+                if j == k:
+                    jac[:, :, i, j, k] = deriv
+                else:
+                    jac[:, :, i, j, k] = 0.5 * deriv
+                    jac[:, :, i, k, j] = 0.5 * deriv
     return jac
 
 
@@ -84,13 +90,10 @@ def _check_vjp_against_fd(centers, covs, scale, rng):
     np.testing.assert_allclose(grad_X, expected_grad_X, rtol=1e-5, atol=1e-8)
 
     # --- Check grad_cov via FD ---
-    # The FD perturbs both (j,k) and (k,j) together (symmetric perturbation),
-    # so the FD derivative is grad_cov[j,k] + grad_cov[k,j]. We symmetrize
-    # our VJP output before comparing.
+    np.testing.assert_allclose(grad_cov, np.swapaxes(grad_cov, -2, -1), atol=1e-12)
     jac_cov = _numerical_jacobian_cov(centers, covs, scale)
     expected_grad_cov = np.einsum("ab,abijk->ijk", g_coefs, jac_cov)
-    grad_cov_sym = grad_cov + np.swapaxes(grad_cov, -2, -1)
-    np.testing.assert_allclose(grad_cov_sym, expected_grad_cov, rtol=1e-5, atol=1e-8)
+    np.testing.assert_allclose(grad_cov, expected_grad_cov, rtol=1e-5, atol=1e-8)
 
 
 # ---------------------------------------------------------------------------