Newbie question - complex matrix algebra using Eigen

user_10449731 0
Viewed 322

I am new to C++ and am trying to translate the following MATLAB code into C++:

x = [1,2,3];
y = [4,5,7];
px = 0.5;
py = 0.5;

m = -1:1:1;

v = [1+i,2+i,3+i;1+2i,2+2i,3+2i;1+3i,2+3i,3+3i];

Tr5 = real(exp(2i*pi/px*x(1)*m)*v*exp(2i*pi/py*m'*y(1)));

fprintf('Tr5 = %f\n',Tr5);

which returns Tr5=18.

Here is my attempt in C++, where I have done as little as possible on each line to try and hunt down the bug:

#include <Eigen/Dense>
#include <unsupported/Eigen/MatrixFunctions>
#include <math.h>
#include <iterator>
#include <vector>
#include <algorithm>
#include <complex>
#include <iostream>

int main()
{

double x [3] = {1, 2, 3};
double y [3] = {4, 5, 7};
const std::complex<double> im(0.0, 1.0);
double px = 0.5;
double py = 0.5;

Eigen::MatrixXcd Tr1;
Eigen::MatrixXcd Tr2;
Eigen::MatrixXcd Tr3;
Eigen::MatrixXcd Tr4;

Eigen::MatrixXd m = Eigen::VectorXd::LinSpaced(2*1+1, -1, 1);
Eigen::MatrixXcd v4(3,3);
v4(0,0) = 1.0d + 1.0 * im;
v4(0,1) = 2.0d + 1.0 * im;
v4(0,2) = 3.0d + 1.0 * im;
v4(1,0) = 1.0d + 2.0 * im;
v4(1,1) = 2.0d + 2.0 * im;
v4(1,2) = 3.0d + 2.0 * im;
v4(2,0) = 1.0d + 3.0 * im;
v4(2,1) = 2.0d + 3.0 * im;
v4(2,2) = 3.0d + 3.0 * im;

Tr1 = 2.0*im*M_PI/px*x[1]*m;
Tr1 = Tr1.exp();

Tr2 = 2.0*im*M_PI/py*y[1]*m;
Tr2 = Tr2.exp();

Tr3 = Tr2.transpose();
Tr4 = Tr1 * v4 * Tr3;

Eigen::MatrixXcd Tr5 = Tr4.real(); // This should now be "18"

std::cout << "Tr5 = " << Tr5 << std::endl;

}

However, when I try to compile and run, I get the error:

/usr/local/include/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h:436: const Eigen::MatrixExponentialReturnValue<Derived> Eigen::MatrixBase<Derived>::exp() const [with Derived = Eigen::Matrix<std::complex<double>, -1, -1>]: Assertion `rows() == cols()' failed.

Could somebody please help me with the code?

Many thanks, Tom Waits

eigen c++
2 Answers

It seems that you are not using coefficient-wise exponential but matrix exponential (which is included from <unsupported/Eigen/MatrixFunctions>) Matrix exponential (see https://en.wikipedia.org/wiki/Matrix_exponential) is a matrix function on square matrices. You see that assertion error because your matrix 'm' is not square.

You should use coefficient-wise exponential. In Eigen, coefficient-wise functions are defined for arrays but they are also available for matrices by first casting it as an array : m.array() (see https://eigen.tuxfamily.org/dox/group__CoeffwiseMathFunctions.html) However, you'll need to cast your arrays back to matrix before performing matrix multiplication: leftSide.array().exp().matrix() * v * rightSide.array().exp().matrix()

Finally you can fill your matrix v using comma initializer syntax (which I believe ise easier to read). I will show that just as an example and of course, your initialization is also fine.

Example (with a syntax similar to MATLAB):

#define _USE_MATH_DEFINES

#include "Eigen/Core"
#include <iostream>
#include <math.h>

std::complex<double> operator"" _i(long double x ) 
{ 
    return std::complex<double>(0.0, x);
} 

int main()
{
    double x [3] = {1, 2, 3};
    double y [3] = {4, 5, 7};

    double px = 0.5;
    double py = 0.5;

    Eigen::RowVectorXd m = Eigen::RowVectorXd::LinSpaced(3, -1, 1);

    Eigen::MatrixXcd v(3,3);
    v << (1.0 + 1.0_i), (2.0 + 1.0_i), (3.0 + 1.0_i),
         (1.0 + 2.0_i), (2.0 + 2.0_i), (3.0 + 2.0_i),
         (1.0 + 3.0_i), (2.0 + 3.0_i), (3.0 + 3.0_i);

    double Tr5 = ((2.0_i * M_PI / px * x[0] * m).array().exp().matrix() * v * (2.0_i * M_PI / py * y[0] * m).array().exp().matrix().transpose()).real()(0);

    std::cout << "Tr5 = " << Tr5 << std::endl;

    return 0;
}

This will yield your desired output:

Tr5 = 18
user_141014400
  1. You forgot to #include <iostream>
  2. Please use std::cout and std::endl rather than cout, endl
  3. You incorrectly assumed that .real() member function returns a real (double) scalar. It returns a matrix, see https://eigen.tuxfamily.org/dox/namespaceEigen.html#title49 . Similarly, in Matlab (Actually I checked in Octave, https://octave-online.net/ ) real(matrix) returns the "real part" of a matrix, element-wise, that is, a matrix, in full agreement with Eigen. So replace
double Tr5;

with

Eigen::MatrixXcd Tr5;

Actually, it is a good practice in C++ to postpone object definitions till they're really needed, so you might consider combining declaration with initialization:

Eigen::MatrixXcd Tr5 = Tr4.real(); // This should now be a real double
user_124512430
Related