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

Đang mở
  • Có gì không ổn?
#3,483 0 bình luận 0 reaction 0 người được giao Xem trên GitHub

Maintainer thường phản hồi trong vòng 1 ngày

Chưa có ai nhận issue này.

Đánh giá

Độ khó
5/5
Thời gian dự kiến
Hơn một tuần
Mức phù hợp với người mới
35/100
Loại issue
Tính năng
Độ rõ ràng
Đặc tả rõ ràng
Mức độ hoạt động
Sôi nổi
Công nghệ
cpp
Lĩnh vực
backend

Hướng nghiên cứu

Read the listed implementations in stan/math/prim/fun/dot_product.hpp, multiply.hpp, dot_self.hpp, columns_dot_product.hpp, and rows_dot_product.hpp, along with the Eigen dot documentation cited in the issue. First decide which complex dot product convention should apply consistently, then document that contract in the relevant comments and Functions Reference. Done means the affected functions follow the documented convention and the row-vector multiply agrees with the equivalent matrix multiply.

Do mô hình lập chỉ mục viết ra từ nội dung của issue.

Mô tả

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):

Function Implementation Complex result
dot_product(u, v) (vectors, row vectors, arrays) u.dot(v) (stan/math/prim/fun/dot_product.hpp:28; the array overload delegates to it at :70) $\sum_i \bar{u}_i v_i$
multiply(rv, v) forwards to dot_product (stan/math/prim/fun/multiply.hpp:181) $\sum_i \bar{rv}_i v_i$
dot_self(v) squaredNorm() (stan/math/prim/fun/dot_self.hpp:36) $\sum_i \lvert v_i \rvert^2$
columns_dot_self(x), rows_dot_self(x) colwise() / rowwise().squaredNorm() (stan/math/prim/fun/columns_dot_self.hpp:24, stan/math/prim/fun/rows_dot_self.hpp:22) $\sum \lvert x_{ij} \rvert^2$
columns_dot_product(a, b), rows_dot_product(a, b) cwiseProduct(...).colwise() / rowwise().sum() (stan/math/prim/fun/columns_dot_product.hpp:31, stan/math/prim/fun/rows_dot_product.hpp:31) $\sum a_{ij} b_{ij}$, no conjugate
multiply(A, B), every other complex matrix or vector shape lazyProduct (stan/math/prim/fun/multiply.hpp:156) matrix product, no conjugate

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 579253050c)

Ngôn ngữ chính
C++
Star
845
Fork
225
Merge trung bình
3 ngày 9 giờ
Pull request đã merge (30 ngày)
12

Chuẩn bị môi trường

Bắt đầu từ đâu

  1. Đọc hết issue, rồi đọc hướng dẫn đóng góp của dự án.
  2. Bình luận trên issue rằng bạn sẽ nhận — tránh hai người làm cùng một việc.
  3. Fork repository và làm thay đổi trên một nhánh.
  4. Mở pull request có tham chiếu số hiệu của issue.

Issue khác của stan-dev/math

Tất cả issue của stan-dev/math

Issue tương tự

Thêm issue về C++

Nhận issue mới trong hộp thư của bạn

Bản tóm tắt ngắn những issue GitHub phù hợp với người mới.