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

Aberta
  • Algo errado?
#3,483 0 comentários 0 reações 0 responsáveis Ver no GitHub

Mantenedores costumam responder em até 1 dia

Ninguém assumiu esta issue ainda.

Avaliação

Dificuldade
5/5
Tempo estimado
Mais de uma semana
Facilidade para iniciantes
35/100
Tipo de issue
Funcionalidade
Clareza
Claramente especificada
Status de atividade
Ativa
Stack de tecnologia
cpp
Domínio
backend

Direção de pesquisa

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.

Escrita pelo modelo de indexação a partir do texto da issue.

Descrição

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)

Linguagem predominante
C++
Estrelas
845
Forks
225
Merge médio
4d 21h
PRs com merge (30d)
13

Preparar o ambiente

Primeiros passos

  1. Leia a issue inteira e depois o guia de contribuição do projeto.
  2. Comente na issue dizendo que vai assumir — evita que duas pessoas façam o mesmo trabalho.
  3. Faça um fork do repositório e trabalhe em uma branch.
  4. Abra um pull request que referencie o número da issue.

Mais de stan-dev/math

Todas as issues de stan-dev/math

Issues semelhantes

Mais issues de C++

Receba novas issues na sua caixa de entrada

Um resumo curto de issues do GitHub para quem está começando.