Complex dot_product and multiply(row_vector, vector) conjugate the first argument; columns_dot_product and rows_dot_product do not #3483

This is not really a bug per-se, but more like an inconsistency and a lack of convention. Before resolving this in code, a decision needs to be made for what the convention should be and document that in the relevant place.

Description

For complex arguments the dot-product functions disagree on whether the first argument is conjugated, and multiply(rv, v) disagrees with every other complex shape of multiply. The prim code behaves as follows (the rev overloads take only var scalars, e.g. stan/math/rev/fun/multiply.hpp:27, and stan/math/rev/fun/dot_product.hpp:33 excludes complex explicitly; the only fwd overload among these, multiply, takes only fvar value types, stan/math/fwd/fun/multiply.hpp:17; so complex<var> and complex<fvar<T>> use the same prim code):

The conjugation comes from Eigen: MatrixBase::dot is documented as "conjugate-linear in the first variable" (lib/eigen_5.0.1/Eigen/src/Core/Dot.h:38-40, dot_impl uses Conj = true at lib/eigen_5.0.1/Eigen/src/Core/InnerProduct.h:255). The two C++-only overloads that Stan does not expose do not conjugate: dot_product(v1, v2, length) (stan/math/prim/fun/dot_product.hpp:47) and scalar dot_self(z), which returns z * z (stan/math/prim/fun/dot_self.hpp:15).

As a result, in Stan rv * v differs from to_matrix(rv) * v, and dot_product(u, v) differs from columns_dot_product(u, v) for column vectors u and v. Also, dot_product(v, v) equals dot_self(v), but columns_dot_product(v, v) does not equal columns_dot_self(v). For real arguments all of these agree. Apart from "squared norm" in the doc comment of dot_self (stan/math/prim/fun/dot_self.hpp:27), neither the doc comments nor the Functions Reference say which convention applies; the Functions Reference says only "the dot product of x and y", and the doc comment of multiply(rv, v) says "The return is the same as the dot product" (stan/math/prim/fun/multiply.hpp:164-166). Only multiply has complex unit tests, and its complex multiply(row_vector, vector) test uses real-valued entries (test/unit/math/mix/fun/multiply_complex_test.cpp:27, :36, :43), for which both conventions agree.

The contract needs a decision: should the complex dot products conjugate their first argument (Hermitian inner product) or not (bilinear)? There is precedent for both: NumPy's vdot conjugates and its dot does not.

Example

#include <stan/math/prim.hpp>
#include <iostream>

int main() {
  using cd = std::complex<double>;
  Eigen::RowVectorXcd rv(2);
  rv << cd(2, -1), cd(0.25, 0.75);
  Eigen::VectorXcd v(2);
  v << cd(1, 2), cd(-0.5, 1.5);
  Eigen::VectorXcd u = rv.transpose();
  Eigen::MatrixXcd m = stan::math::multiply(stan::math::to_matrix(rv), v);
  Eigen::RowVectorXcd c = stan::math::columns_dot_product(u, v);
  Eigen::RowVectorXcd cvv = stan::math::columns_dot_product(v, v);
  Eigen::RowVectorXd cs = stan::math::columns_dot_self(v);
  std::cout << stan::math::multiply(rv, v) << "\n"     // (1,5.75)   conjugated
            << m(0, 0) << "\n"                          // (2.75,3)   plain
            << stan::math::dot_product(rv, v) << "\n"  // (1,5.75)   conjugated
            << c(0) << "\n"                             // (2.75,3)   plain
            << stan::math::dot_self(v) << "\n"          // (7.5,0)    conjugated
            << cs(0) << "\n"                            // 7.5        conjugated
            << cvv(0) << "\n";                          // (-5,2.5)   plain
}

For $rv = (2-i, 0.25+0.75i)$ and $v = (1+2i, -0.5+1.5i)$ the conjugated product is $(2+i)(1+2i) + (0.25-0.75i)(-0.5+1.5i) = 5i + (1+0.75i) = 1+5.75i$. The plain product is $(2-i)(1+2i) + (0.25+0.75i)(-0.5+1.5i) = (4+3i) - 1.25 = 2.75+3i$.

Expected Output

One documented convention that all these functions follow, and multiply(rv, v) equal to multiply(to_matrix(rv), v)(0, 0), as for every other shape of multiply. If the dot products are meant to be Hermitian, columns_dot_product and rows_dot_product would conjugate their first argument, and multiply(rv, v) would stop delegating to dot_product for complex input. If they are meant to be bilinear, dot_product, dot_self, columns_dot_self and rows_dot_self would stop conjugating, and multiply(rv, v) would then agree automatically.

Current Version:

v5.4.0 (develop at 5792530)