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

Aperta
  • Qualcosa non va?
#3,483 0 commenti 0 reazioni 0 assegnatari Vedi su GitHub

I maintainer di solito rispondono entro 1 giorno

Nessuno ha ancora preso questa issue.

Valutazione

Difficoltà
5/5
Tempo stimato
Più di una settimana
Idoneità per principianti
35/100
Tipo di issue
Funzionalità
Chiarezza
Specificata chiaramente
Stato di attività
Attiva
Stack tecnologico
cpp
Ambito
backend

Direzione di ricerca

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.

Scritto dal modello di indicizzazione a partire dal testo della issue.

Descrizione

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)

Lingua principale
C++
Stelle
845
Fork
225
Merge medio
3g 9h
PR unite (30g)
12

Preparare l'ambiente

Come iniziare

  1. Leggi tutta la issue e poi la guida ai contributi del progetto.
  2. Commenta sulla issue per dire che te ne occupi tu — evita che due persone facciano lo stesso lavoro.
  3. Fai un fork del repository e lavora su un branch.
  4. Apri una pull request che faccia riferimento al numero della issue.

Altre issue di stan-dev/math

Tutte le issue di stan-dev/math

Issue simili

Altre issue su C++

Ricevi le nuove issue nella tua casella

Un breve riepilogo di issue GitHub adatte ai principianti.