[eigen] Unexpected behavior of applyHouseholderOnTheRight()

[Yixuan Qiu <yixuanq@xxxxxxxxx>]

Hi All,

I have a question on the MatrixBase::applyHouseholderOnTheRight() function that applies a Householder transformation to a matrix on the right.

If the Householder transformation is H, then I expect that this function will calculate X * H or X * H^*, depending on the definition. But in fact I find that it calculates X * H^T, which is fine with real-valued matrices but is quite confusing with complex-valued matrices.

Attached is a minimal example to reproduce the issue. I wonder if this is a "bug" or my own misunderstanding of the function.

Thanks.

Best,

Yixuan

--

Yixuan Qiu <yixuanq@xxxxxxxxx>
Department of Statistics,
Purdue University

#include <Eigen/Eigenvalues>
#include <iostream>

using Eigen::VectorXcd;
using Eigen::MatrixXcd;

int main()
{
    const int n = 3;
    std::srand(0);
    VectorXcd x = VectorXcd::Random(n);
    VectorXcd essential(n - 1);
    std::complex<double> tau;
    double beta;
    // Create a random Householder transformation
    x.makeHouseholder(essential, tau, beta);
    
    std::cout << "essential = \n" << essential << std::endl << std::endl;
    std::cout << "tau = " << tau << std::endl << std::endl;
    std::cout << "beta = " << beta << std::endl << std::endl;
    
    // Compute the transformation matrix by definition
    // H = 1 - tau * v * v^*, v = [1 essential^T]^T
    VectorXcd h(n);
    h[0] = std::complex<double>(1, 0);
    h.tail(n - 1) = essential;
    MatrixXcd H = MatrixXcd::Identity(n, n) - tau * h * h.adjoint();
    std::cout << "H = \n" << H << std::endl << std::endl;
    
    // Compute the results of applying the transformation to an identity matrix
    MatrixXcd H_left = MatrixXcd::Identity(n, n);
    MatrixXcd H_right = MatrixXcd::Identity(n, n);
    VectorXcd tmp(n);
    
    H_left.applyHouseholderOnTheLeft(essential, tau, tmp.data());
    H_right.applyHouseholderOnTheRight(essential, tau, tmp.data());
    
    // It shows that H_left = H, but H_right = H^T (not H or H^*)
    std::cout << "H_left = \n" << H_left << std::endl << std::endl;
    std::cout << "H_right = \n" << H_right << std::endl << std::endl;
    
    return 0;
}