Discussion:
[eigen] using views for symmetric matricies
Smith, Louis
2018-03-05 18:39:43 UTC
Permalink
Hello,


I have a very basic question about views. Thanks in advance for being understanding, I'm sure this a totally elementary problem that just shows my lack of understanding, but this is the appropriate forum for getting help with this library. Right?*


I'm trying to use views the way I perceive that they are being used in the quick reference <https://eigen.tuxfamily.org/dox/group__QuickRefPage.html> to, for example, only compute and store half the entries in a symmetric matrix, then have later operations view the matrix produced in this way as symmetric. I'm working with things this way as a substitute for writing nested for-loops where the inner loop's index is constrained to rolling over only one of the triangles of the array I am either getting compile time or runtime errors for all of it, with code that works if the triangularView and selfadjointView statements I included are removed. The question is how can I use views to only compute and store half of a symmetric matrix, while still using the visitors and reductions that are part of the Matrix class? Following are some examples I tried that do not work.


I use the following function to get a matrix of pairwise distances:

// takes a nxd data matrix, returns an nxn matrix containing pairwise distances
// use formula (a - b)^2 = a^2 + b^2 -2a*b.
MatrixXd pairwise_dists(const Ref<const MatrixXd>& data){
const VectorXd data_sq = data.rowwise().squaredNorm();
MatrixXd distances;
distances = data_sq.rowwise().replicate(data.rows())
+ data_sq.transpose().colwise().replicate(data.rows())
- 2.*data*data.transpose();
distances.diagonal().setZero(); // prevents nans from occurring along diag.
distances = distances.cwiseSqrt();
return distances;
}

I would like to do something like:

// takes a nxd data matrix, returns an nxn matrix containing pairwise distances
// use formula (a - b)^2 = a^2 + b^2 -2a*b.
MatrixXd pairwise_dists(const Ref<const MatrixXd>& data){
const VectorXd data_sq = data.rowwise().squaredNorm();
MatrixXd distances;
distances.triangularView<StrictlyLower>() = data_sq.rowwise().replicate(data.rows())
+ data_sq.transpose().colwise().replicate(data.rows())
- 2.*data*data.transpose();
distances = distances.cwiseSqrt();
return distances;
}


And then when I use this function later, for example in a reduction:

VectorXd rho = (-0.5*pairwise_dists(scdata)).array().exp().rowwise().sum();

Would become something like:

VectorXd rho = (-0.5*pairwise_dists(scdata).selfadjointView<StrictlyLower,ZeroDiag>()).array().exp().rowwise().sum();


Adding the triangular view statement to the function 'pairwise_dists' compiles but gives a runtime error whenever the function is called by my code, whereas using the 'selfadjointView' statement I give causes a compile error. Changing the template parameters to only <StrictlyLower> or changing both template parameters to <Lower> do not fix the issue.


Even in a toy example I seem to be having issues:

Matrix3f m0, m1, m2, m3, m4;
m2 << 1,2,3,
4,5,6,
7,8,9;
m3 << 1,0,1,
3,0,0,
0,4,0;

m0 = m2 + m3;
m1.triangularView<StrictlyLower>() = m2+m3;
cout << m0 << "\nNo triangular view\n";
cout << m1 << "\nS. Lower triangular view\n";
m4 = m1.selfadjointView<StrictlyLower>(); //+ m3;
cout << m4 << "\nself-adjoint lower view\n";

The above code will compile and run, but m4 is not a dense symmetric matrix like I'd expect it to be. It has uninitialised values on the diagonal and upper triangle, as though selfadjointView didn't do anything.


If I uncomment the addition operation on the line defining m4, the code no longer compiles.


I hope this is enough information to germinate a response. Thanks in advance for taking the time to read it.


Regards,
Louis


*Sorry if this is not true. It really seems like many of the posts to this mailing list are higher level discussions, comments on use, or bug reports, so if there's a more correct place to mail basic use questions I'd happily direct subsequent questions there.
Gael Guennebaud
2018-03-08 16:38:18 UTC
Permalink
Hi,

let me answer a few shortcomings.

1 - In:

m1.triangularView<StrictlyLower>() = m2+m3;

triangularView on the left-hand-side of an assignment works as a writing
mask, so the diagonal and upper-left triangular part of m1 is left
unchanged. Either write:

m1.setZero();
m1.triangularView<StrictlyLower>() = m2+m3;

or

m1 = (m2+m3).triangularView<StrictlyLower>();

2 - The line:

m4 = m1.selfadjointView<StrictlyLower>();

is invalid and should give you a compiler error (a static assert seems to
be missing in Eigen) because selfadjointView only accept Lower of Upper
parameters, no Strictly* nor Unit* stuff. If the diagonal of m1 is already
zero, then using Lower is fine.

I hope these two explanations answer your questions.

Gael.
Post by Smith, Louis
Hello,
I have a very basic question about views. Thanks in advance for being
understanding, I'm sure this a totally elementary problem that just shows
my lack of understanding, but this is the appropriate forum for getting
help with this library. Right?*
I'm trying to use views the way I perceive that they are being used in the quick
reference <https://eigen.tuxfamily.org/dox/group__QuickRefPage.html>to,
for example, only compute and store half the entries in a symmetric matrix,
then have later operations view the matrix produced in this way as
symmetric. I'm working with things this way as a substitute for writing
nested for-loops where the inner loop's index is constrained to rolling
over only one of the triangles of the array I am either getting compile
time or runtime errors for all of it, with code that works if the
triangularView and selfadjointView statements I included are removed. The
question is how can I use views to only compute and store half of a
symmetric matrix, while still using the visitors and reductions that are
part of the Matrix class? Following are some examples I tried that do not
work.
// takes a nxd data matrix, returns an nxn matrix containing pairwise distances
// use formula (a - b)^2 = a^2 + b^2 -2a*b.
MatrixXd pairwise_dists(const Ref<const MatrixXd>& data){
const VectorXd data_sq = data.rowwise().squaredNorm();
MatrixXd distances;
distances = data_sq.rowwise().replicate(data.rows())
+ data_sq.transpose().colwise().replicate(data.rows())
- 2.*data*data.transpose();
distances.diagonal().setZero(); // prevents nans from occurring along diag.
distances = distances.cwiseSqrt();
return distances;
}
// takes a nxd data matrix, returns an nxn matrix containing pairwise distances
// use formula (a - b)^2 = a^2 + b^2 -2a*b.
MatrixXd pairwise_dists(const Ref<const MatrixXd>& data){
const VectorXd data_sq = data.rowwise().squaredNorm();
MatrixXd distances;
distances.triangularView<StrictlyLower>() = data_sq.rowwise().replicate(
data.rows())
+ data_sq.transpose().colwise().replicate(data.rows())
- 2.*data*data.transpose();
distances = distances.cwiseSqrt();
return distances;
}
VectorXd rho = (-0.5*pairwise_dists(scdata)).array().exp().rowwise().sum
();
VectorXd rho = (-0.5*pairwise_dists(scdata).selfadjointView<StrictlyLower,
ZeroDiag>()).array().exp().rowwise().sum();
Adding the triangular view statement to the function 'pairwise_dists'
compiles but gives a runtime error whenever the function is called by my
code, whereas using the 'selfadjointView' statement I give causes a compile
error. Changing the template parameters to only <StrictlyLower> or changing
both template parameters to <Lower> do not fix the issue.
Matrix3f m0, m1, m2, m3, m4;
m2 << 1,2,3,
4,5,6,
7,8,9;
m3 << 1,0,1,
3,0,0,
0,4,0;
m0 = m2 + m3;
m1.triangularView<StrictlyLower>() = m2+m3;
cout << m0 << "\nNo triangular view\n";
cout << m1 << "\nS. Lower triangular view\n";
m4 = m1.selfadjointView<StrictlyLower>(); //+ m3;
cout << m4 << "\nself-adjoint lower view\n";
The above code will compile and run, but m4 is not a dense symmetric
matrix like I'd expect it to be. It has uninitialised values on the
diagonal and upper triangle, as though selfadjointView didn't do anything.
If I uncomment the addition operation on the line defining m4, the code no longer compiles.
I hope this is enough information to germinate a response. Thanks in
advance for taking the time to read it.
Regards,
Louis
*Sorry if this is not true. It really seems like many of the posts to this
mailing list are higher level discussions, comments on use, or bug reports,
so if there's a more correct place to mail basic use questions I'd happily
direct subsequent questions there.
Christoph Hertzberg
2018-03-08 16:52:26 UTC
Permalink
Hi!
Post by Smith, Louis
m4 = m1.selfadjointView<StrictlyLower>();
is invalid and should give you a compiler error (a static assert seems to
be missing in Eigen) because selfadjointView only accept Lower of Upper
parameters, no Strictly* nor Unit* stuff. If the diagonal of m1 is already
zero, then using Lower is fine.
I'm not sure if there are many applications for that, but having a
selfadjointView with a zero diagonal might have some use cases.

What also could have some use cases is to have a skew-symmetric or
skew-hermitian view on a matrix (i.e., "the other half" equals
-A.transpose() or -A.adjoint())


Christoph
Post by Smith, Louis
I hope these two explanations answer your questions.
Gael.
Post by Smith, Louis
Hello,
I have a very basic question about views. Thanks in advance for being
understanding, I'm sure this a totally elementary problem that just shows
my lack of understanding, but this is the appropriate forum for getting
help with this library. Right?*
I'm trying to use views the way I perceive that they are being used in the quick
reference <https://eigen.tuxfamily.org/dox/group__QuickRefPage.html>to,
for example, only compute and store half the entries in a symmetric matrix,
then have later operations view the matrix produced in this way as
symmetric. I'm working with things this way as a substitute for writing
nested for-loops where the inner loop's index is constrained to rolling
over only one of the triangles of the array I am either getting compile
time or runtime errors for all of it, with code that works if the
triangularView and selfadjointView statements I included are removed. The
question is how can I use views to only compute and store half of a
symmetric matrix, while still using the visitors and reductions that are
part of the Matrix class? Following are some examples I tried that do not
work.
// takes a nxd data matrix, returns an nxn matrix containing pairwise distances
// use formula (a - b)^2 = a^2 + b^2 -2a*b.
MatrixXd pairwise_dists(const Ref<const MatrixXd>& data){
const VectorXd data_sq = data.rowwise().squaredNorm();
MatrixXd distances;
distances = data_sq.rowwise().replicate(data.rows())
+ data_sq.transpose().colwise().replicate(data.rows())
- 2.*data*data.transpose();
distances.diagonal().setZero(); // prevents nans from occurring along diag.
distances = distances.cwiseSqrt();
return distances;
}
// takes a nxd data matrix, returns an nxn matrix containing pairwise distances
// use formula (a - b)^2 = a^2 + b^2 -2a*b.
MatrixXd pairwise_dists(const Ref<const MatrixXd>& data){
const VectorXd data_sq = data.rowwise().squaredNorm();
MatrixXd distances;
distances.triangularView<StrictlyLower>() = data_sq.rowwise().replicate(
data.rows())
+ data_sq.transpose().colwise().replicate(data.rows())
- 2.*data*data.transpose();
distances = distances.cwiseSqrt();
return distances;
}
VectorXd rho = (-0.5*pairwise_dists(scdata)).array().exp().rowwise().sum
();
VectorXd rho = (-0.5*pairwise_dists(scdata).selfadjointView<StrictlyLower,
ZeroDiag>()).array().exp().rowwise().sum();
Adding the triangular view statement to the function 'pairwise_dists'
compiles but gives a runtime error whenever the function is called by my
code, whereas using the 'selfadjointView' statement I give causes a compile
error. Changing the template parameters to only <StrictlyLower> or changing
both template parameters to <Lower> do not fix the issue.
Matrix3f m0, m1, m2, m3, m4;
m2 << 1,2,3,
4,5,6,
7,8,9;
m3 << 1,0,1,
3,0,0,
0,4,0;
m0 = m2 + m3;
m1.triangularView<StrictlyLower>() = m2+m3;
cout << m0 << "\nNo triangular view\n";
cout << m1 << "\nS. Lower triangular view\n";
m4 = m1.selfadjointView<StrictlyLower>(); //+ m3;
cout << m4 << "\nself-adjoint lower view\n";
The above code will compile and run, but m4 is not a dense symmetric
matrix like I'd expect it to be. It has uninitialised values on the
diagonal and upper triangle, as though selfadjointView didn't do anything.
If I uncomment the addition operation on the line defining m4, the code no
longer compiles.
I hope this is enough information to germinate a response. Thanks in
advance for taking the time to read it.
Regards,
Louis
*Sorry if this is not true. It really seems like many of the posts to this
mailing list are higher level discussions, comments on use, or bug reports,
so if there's a more correct place to mail basic use questions I'd happily
direct subsequent questions there.
--
Dr.-Ing. Christoph Hertzberg

Besuchsadresse der Nebengeschäftsstelle:
DFKI GmbH
Robotics Innovation Center
Robert-Hooke-Straße 5
28359 Bremen, Germany

Postadresse der Hauptgeschäftsstelle Standort Bremen:
DFKI GmbH
Robotics Innovation Center
Robert-Hooke-Straße 1
28359 Bremen, Germany

Tel.: +49 421 178 45-4021
Zentrale: +49 421 178 45-0
E-Mail: ***@dfki.de

Weitere Informationen: http://www.dfki.de/robotik
-----------------------------------------------------------------------
Deutsches Forschungszentrum fuer Kuenstliche Intelligenz GmbH
Firmensitz: Trippstadter Straße 122, D-67663 Kaiserslautern
Geschaeftsfuehrung: Prof. Dr. Dr. h.c. mult. Wolfgang Wahlster
(Vorsitzender) Dr. Walter Olthoff
Vorsitzender des Aufsichtsrats: Prof. Dr. h.c. Hans A. Aukes
Amtsgericht Kaiserslautern, HRB 2313
Sitz der Gesellschaft: Kaiserslautern (HRB 2313)
USt-Id.Nr.: DE 148646973
Steuernummer: 19/672/50006
-----------------------------------------------------------------------
Loading...