Re: [eigen] discrepancy in (triangular view + transpose) and self-adjointness

[Manoj Rajagopalan <rmanoj@xxxxxxxxx>]

Attached compilable C++ program.

thanks,
Manoj


On Wednesday 14 July 2010 07:14:02 pm Benoit Jacob wrote:
> Hi,
>
> It would help if you could make a compilable test program, so we can
> quickly run stuff.
>
> Benoit
>
> 2010/7/14 Manoj Rajagopalan <rmanoj@xxxxxxxxx>:
> > Hi all,
> >
> >   Could someone please help me figure out the source of the following
> > discrepancy?
> >
> >   Given:
> >        int const N = 5;
> >        MatrixXi A(N,N);
> >        A.setRandom();
> >
> >        TriangularView<MatrixXi,Lower> const A_l(A);
> >        TriangularView<MatrixXi,Lower|ZeroDiag> const A_zl(A);
> >
> > the following must be true:
> >
> >   A.selfadjointView<Lower>() == A_l + A_zl.transpose()
> >
> > Basically, I am adding the lower triangle of A to its own transpose,
> > without adding the diagonal a second time. So the following is calculated
> > correctly by Eigen:
> >
> >    MatrixXi C1(N,N), C2(N,N);
> >    C1 = A.selfadjointView<Lower>();
> >    C2 = A_l.toDenseMatrix() + A_zl.transpose().toDenseMatrix();
> >    assert(C1 == C2);
> >
> > But, the following code produces a discrepancy:
> >
> >    MatrixXi B(N,N);
> >    B.setRandom();
> >    C1 = A.selfadjointView<Lower>() * B;
> >    C2 = (A_l * B)   +   (A_zl.transpose() * B);
> >    assert(C1 == C2); // <--------  FALSE!
> >
> >     After a long day my brain is unable to figure out if I am doing
> > something extremely silly - it seems too trivial to be a bug. Can someone
> > offer comments?
> >
> > thanks,
> > Manoj

#include <Eigen/Core>
#include <iostream>

using namespace Eigen;
using namespace std;

int main(void)
{
	int const N = 5;
	MatrixXi A(N,N), B(N,N), C1(N,N), C2(N,N), C_diff(N,N);
	A.setRandom();
	B.setRandom();
	C1 = A.selfadjointView<Lower>();
	// C1 = A.selfadjointView<Lower>() * B;

	TriangularView<MatrixXi,Lower> const A_l(A);
	TriangularView<MatrixXi,Lower|ZeroDiag> const A_zl(A);
	// C2 = A_l * B;
	// C2 += A_zl.transpose() * B;
	C2 = A_l.toDenseMatrix() + A_zl.transpose().toDenseMatrix();
	C_diff = C1 - C2;

	cout << "A =\n" << A << endl;
	cout << "\nA.triangularView<Lower>() =\n" << A_l.toDenseMatrix() << endl;
	cout << "\nA.triangularView<Lower|ZeroDiag>() =\n" << A_zl.toDenseMatrix() << endl;
	cout << "C1 =\n" << C1 << endl;
	cout << "C2 =\n" << C2 << endl;
	cout << "C_diff =\n" << C_diff << endl;

	assert(C1 == C2);
	return 0;
}