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

Offen
  • Stimmt etwas nicht?
#3,483 0 Kommentare 0 Reaktionen 0 zugewiesene Personen Auf GitHub ansehen

Maintainer antworten meist innerhalb von 1 Tag

Dieses Issue hat noch niemand übernommen.

Bewertung

Schwierigkeit
5/5
Geschätzter Aufwand
Über eine Woche
Anfängerfreundlichkeit
35/100
Issue-Typ
Feature
Klarheit
Klar beschrieben
Aktivitätsstatus
Aktiv
Tech-Stack
cpp
Bereich
backend

Rechercherichtung

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.

Vom Indexierungsmodell aus dem Issue-Text verfasst.

Beschreibung

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)

Vorherrschende Sprache
C++
Sterne
845
Forks
225
Ø Merge
4 T. 21 Std.
Gemergte PRs (30 T.)
13

Entwicklungsumgebung

Erste Schritte

  1. Lesen Sie das ganze Issue und danach den Beitragsleitfaden des Projekts.
  2. Schreiben Sie ins Issue, dass Sie es übernehmen — das erspart doppelte Arbeit.
  3. Forken Sie das Repository und arbeiten Sie in einem Branch.
  4. Öffnen Sie einen Pull Request, der die Issue-Nummer nennt.

Mehr aus stan-dev/math

Alle Issues in stan-dev/math

Ähnliche Issues

Weitere Issues zu C++

Neue Issues direkt in Ihr Postfach

Eine kurze Übersicht über anfängerfreundliche GitHub-Issues.