diff --git a/c++/cppdlr/dlr_imtime.hpp b/c++/cppdlr/dlr_imtime.hpp
index 0bd94a9..1ca7729 100644
--- a/c++/cppdlr/dlr_imtime.hpp
+++ b/c++/cppdlr/dlr_imtime.hpp
@@ -534,6 +534,51 @@ namespace cppdlr {
       }
     }
 
+    /** 
+    * @brief Compute inner product of two imaginary time Green's functions
+    *
+    * We define the inner product of complex matrix-valued f and g as
+    *
+    * (f,g) = 1/beta * sum_{ij} int_0^beta dt conj(f_ij(t)) g_ij(t), 
+    *
+    * where conj refers to complex conjugation. This method takes
+    * the DLR coefficients of f and g as input and returns the inner product.
+    *
+    * We use the numerically stable method described in Appendix B of 
+    *
+    * H. LaBollita, J. Kaye, A. Hampel, "Stabilizing the calculation of the
+    * self-energy in dynamical mean-field theory using constrained residual
+    * minimization," arXiv:2310.01266 (2023).
+    *
+    * @param[in] fc DLR coefficients of f
+    * @param[in] gc DLR coefficients of g
+    *
+    * @return Inner product of f and g
+    * */
+    template <nda::MemoryArray T, nda::Scalar S = nda::get_value_t<T>> std::complex<double> innerprod(T const &fc, T const &gc) const {
+
+      if (r != fc.shape(0) || r != gc.shape(0)) throw std::runtime_error("First dim of input arrays must be equal to DLR rank r.");
+      if (fc.shape() != gc.shape()) throw std::runtime_error("Input arrays must have the same shape.");
+
+      // Initialize inner product matrix, if it hasn't been done already
+      if (ipmat.empty()) { innerprod_init(); }
+
+      std::complex<double> ip = 0;
+      if constexpr (T::rank == 1) { // Scalar-valued Green's function
+        ip = nda::blas::dotc(fc, matvecmul(ipmat, gc));
+      } else if (T::rank == 3) { // Matrix-valued Green's function
+        int norb = fc.shape(1);
+        // Compute inner product
+        for (int i = 0; i < norb; ++i) {
+          for (int j = 0; j < norb; ++j) { ip += nda::blas::dotc(fc(_, i, j), matvecmul(ipmat, gc(_, i, j))); }
+        }
+      } else {
+        throw std::runtime_error("Input arrays must be rank 1 (scalar-valued Green's function) or 3 (matrix-valued Green's function).");
+      }
+
+      return ip;
+    }
+
     /** 
     * @brief Get DLR imaginary time nodes
     *
@@ -668,6 +713,33 @@ namespace cppdlr {
       }
     }
 
+    /**
+    * @brief Initialization for inner product method
+    *
+    * This method is called automatically the first time the innerprod method is
+    * called, but it may also be called manually to avoid the additional
+    * overhead in the first inner product call.
+    */
+    void innerprod_init() const {
+
+      ipmat = nda::matrix<double>(r, r);
+
+      // Matrix of inner product of two DLR expansions
+      double ssum = 0;
+      for (int k = 0; k < r; ++k) {
+        for (int l = 0; l < r; ++l) {
+          ssum = dlr_rf(k) + dlr_rf(l);
+          if (ssum == 0) {
+            ipmat(k, l) = k_it(0.0, dlr_rf(k)) * k_it(0.0, dlr_rf(l));
+          } else if (abs(ssum) < 1) {
+            ipmat(k, l) = -k_it(0.0, dlr_rf(k)) * k_it(0.0, dlr_rf(l)) * std::expm1(-ssum) / ssum;
+          } else {
+            ipmat(k, l) = (k_it(0.0, dlr_rf(k)) * k_it(0.0, dlr_rf(l)) - k_it(1.0, dlr_rf(k)) * k_it(1.0, dlr_rf(l))) / ssum;
+          }
+        }
+      }
+    }
+
     /**
     * @brief Initialization for reflection method
     *
@@ -715,7 +787,10 @@ namespace cppdlr {
     mutable nda::matrix<double> thilb;   ///< "Discrete Hilbert transform" matrix, modified for time-ordered convolution
     mutable nda::matrix<double> ttcf2it; ///< A matrix required for time-ordered convolution
 
-    // Arrays used for dlr_imtime::reflect
+    // Array used for dlr_imtime::innerprod
+    mutable nda::matrix<double> ipmat; ///< Inner product matrix
+
+    // Array used for dlr_imtime::reflect
     mutable nda::matrix<double> refl; ///< Matrix of reflection
 
     // -------------------- hdf5 -------------------
diff --git a/doc/examples.rst b/doc/examples.rst
index 9d425d8..64e3bd3 100644
--- a/doc/examples.rst
+++ b/doc/examples.rst
@@ -192,6 +192,8 @@ list several such use cases below.
 - Perform a "reflection" operation :math:`G(\tau) \mapsto G(\beta-\tau)` on a
   Green's function: see the test ``imtime_ops.refl_matrix`` in the file
   ``test/imtime_ops.cpp``.
+- Compute the inner product of two DLR expansions: see the test
+  ``imtime_ops.innerprod`` in the file ``test/imtime_ops.cpp``.
 - Obtain a DLR expansion by interpolation on the DLR Matsubara frequency nodes:
   see the tests ``imfreq_ops.interp_scalar`` and ``imfreq_ops.interp_matrix`` in
   the file ``test/imfreq_ops.cpp``.
diff --git a/test/c++/imtime_ops.cpp b/test/c++/imtime_ops.cpp
index 2312103..71ef461 100644
--- a/test/c++/imtime_ops.cpp
+++ b/test/c++/imtime_ops.cpp
@@ -1034,6 +1034,78 @@ TEST(imtime_ops, interp_matrix_sym_bos) {
   std::cout << fmt::format("Imag freq: l^2 err = {:e}, L^inf err = {:e}\n", errl2, errlinf);
 }
 
+/**
+* @brief Test inner product of two DLR expansions
+*/
+TEST(imtime_ops, innerprod) {
+
+  double lambda = 200;   // DLR cutoff
+  double eps    = 1e-10; // DLR tolerance
+
+  double beta = 100; // Inverse temperature
+  int norb    = 2;   // Orbital dimensions
+
+  std::cout << fmt::format("eps = {:e}, Lambda = {:e}\n", eps, lambda);
+
+  // Get DLR frequencies
+  auto dlr_rf = build_dlr_rf(lambda, eps);
+  int r       = dlr_rf.size();
+
+  // Get DLR imaginary time object
+  auto itops  = imtime_ops(lambda, dlr_rf);
+  auto dlr_it = itops.get_itnodes();
+
+  // We take the inner product of two matrix-valued functions, each of the form
+  //
+  // c_{1,2} * exp(-tau * om_{1,2})/(1+exp(-beta * om_{1,2}))
+  //
+  // where c_1, c_2 are 2x2 complex-valued matrices, and om_1, om_2 are real
+  // frequencies. The inner product is given by
+  //
+  // (sum_ij conj(c1_ij) * c2_ij) * (1+exp(-beta * om_1))^{-1} * (1+exp(-beta *
+  // om_2))^{-1}  * (1 - exp(-beta * (om_1 + om_2))) / (om_1 + om_2)
+  //
+  // If om_1 + om_2 is positive and not too small, this formula is numerically
+  // stable, and can be used for an analytical reference.
+
+  auto c1    = nda::array<dcomplex, 2>(norb, norb);
+  auto c2    = nda::array<dcomplex, 2>(norb, norb);
+  c1(0, 0)   = 0.12 + 0.21i;
+  c1(0, 1)   = 0.34 + 0.43i;
+  c1(1, 1)   = 0.56 + 0.65i;
+  c1(1, 0)   = conj(c1(0, 1));
+  c2(0, 0)   = -0.22 + 0.11i;
+  c2(0, 1)   = -0.44 + 0.33i;
+  c2(1, 1)   = 0.66 - 0.55i;
+  c2(1, 0)   = conj(c2(0, 1));
+  double om1 = 0.0789;
+  double om2 = 0.456;
+
+  std::complex<double> iptrue = conj(c1(0, 0)) * c2(0, 0) + conj(c1(0, 1)) * c2(0, 1) + conj(c1(1, 0)) * c2(1, 0) + conj(c1(1, 1)) * c2(1, 1);
+  iptrue *= (1.0 - exp(-beta * (om1 + om2))) / (beta * (om1 + om2) * (1.0 + exp(-beta * om1)) * (1.0 + exp(-beta * om2)));
+
+  // Sample Green's functions at DLR nodes
+  auto f = nda::array<dcomplex, 3>(r, norb, norb);
+  auto g = nda::array<dcomplex, 3>(r, norb, norb);
+  for (int i = 0; i < r; ++i) {
+    f(i, _, _) = k_it(dlr_it(i), om1, beta) * c1;
+    g(i, _, _) = k_it(dlr_it(i), om2, beta) * c2;
+  }
+
+  // DLR coefficients of f and g
+  auto fc = itops.vals2coefs(f);
+  auto gc = itops.vals2coefs(g);
+
+  std::complex<double> iptest = itops.innerprod(fc, gc);
+
+  double err = abs(iptest - iptrue);
+  std::cout << fmt::format("True inner product: {:.16e}+1i*{:.16e}\n", iptrue.real(), iptrue.imag());
+  std::cout << fmt::format("Test inner product: {:.16e}+1i*{:.16e}\n", iptest.real(), iptest.imag());
+  std::cout << fmt::format("Inner product error: {:e}\n", err);
+
+  EXPECT_LT(err, eps);
+}
+
 TEST(dlr_imtime, h5_rw) {
 
   double lambda = 1000;  // DLR cutoff