c++ Eigen making copy?

Question

Let A be a symmetric matrix and let v be a vector. I extract from A a block of n columns starting from j and multiply it by v using

VectorXd a;
a = A.middleCols(j,n).selfadjointView<Lower>() * v // does not compile

As this does not compile, whereas this

a = MatrixXd(A.middleCols(j,n).selfadjointView<Lower>()) * v

does, I am wondering whether the second version makes a copy of

A.middleCols(j,n).selfadjointView<Lower>()

or performs the computation directly?

Thanks for any hint.

EDIT: I suspect the problem has something to do with argument types, as I get the error :

invalid argument type 'typename ConstSelfAdjointViewReturnType.... to unary expression'

Indeed, A is the argument of a function passed by const reference using either of

const MatrixXd& A 
const Ref<const MatrixXd>& A

Here is an example:

// this version doesn't compile
MatrixXd returnSymmetricMatrix(const MatrixXd& A, const VectorXd& v, const MatrixXd& B){
// B is a symmetric matrix

VectorXd a;
a = A.middleCols(3, 4).selfadjointView<Lower>() * v;
MatrixXd M(code_fill(){...}); 
// code_fill is the function filling the lower triangular part of a symmetric matrix
M.block(1, 2, 3, 4).triangularView<Lower>() += B.selfadjointView<Lower>();

return M;
}

// this version compiles
MatrixXd returnSymmetricMatrix(const MatrixXd& A, const VectorXd& v, const MatrixXd& B){
// B is a symmetric matrix

VectorXd a;
a = MatrixXd(A.middleCols(3, 4).selfadjointView<Lower>()) * v;
MatrixXd M(code_fill(){...}); 
// code_fill is the function filling the lower triangular part of a symmetric matrix
Matrix(M.block(1, 2, 3, 4).triangularView<Lower>()) += B.selfadjointView<Lower>();

return M;
}

EDIT2 Regarding my initial question and the example I've added at the Edit section, I am a bit confused regarding the copying. As I understand the difference between the working and the non-working versions, the line

Matrix(M.block(1, 2, 3, 4).triangularView<Lower>()) += B.selfadjointView<Lower>();

works because its lhs tells to Eigen that M.block(1, 2, 3, 4).triangularView() is actually a matrix rather than a reference to a matrix. Otherwise, the operator += would through an error that this operator is not overloaded for .block(). So my original question is whether Matrix(...) only tells that it's a Matrix to enable the computation, or rather copy the ... into a Matrix? Thanks!

The expression compiles (and runs) fine for me. Can you double check and add any additional information needed? — Avi Ginsburg, Jan 30 '17 at 12:54
Yes, the first one should work fine, I see it's using `M` instead of `A`... — ggael, Jan 30 '17 at 13:06

ggael · Accepted Answer · 2017-01-30T21:36:32.510

1

The following expression:

A.middleCols(j,n).selfadjointView<Lower>()

does not create any copy.

On the other hand, to avoid a temporary for the result of the product, you can add .noalias():

a.noalias() = M.middleCols(j,n).selfadjointView<Lower>() * v;

This is only needed for the immediate assignment of dense product to allow code like:

a = M * a;

to behave as expected.

EDIT:

regarding your compilation issue, the following compiles fine:

#include <Eigen/Dense>
using namespace Eigen;
int main()
{
  int n = 10;
  MatrixXd M = MatrixXd::Random(n,2*n);
  VectorXd v = VectorXd::Random(n);
  VectorXd a;
  a.noalias() = M.middleCols(2,n).selfadjointView<Upper>() * v;

  const Ref<const MatrixXd>& A = M;
  a.noalias() = A.middleCols(2,n).selfadjointView<Upper>() * v;
}

EDIT2

The following line:

MatrixXd(M.block(1, 2, 3, 4).triangularView<Lower>()) += B.selfadjointView<Lower>();

makes no sense as you are creating a temporary copy that you are assigning to. Recall that here MatrixXd(whatever) calls the Eigen::Matrix constructor. Also assigning a selfadjoint-view to a triangular-view makes no sense. I cannot even think about what could be a reasonable behavior for this. If you only want to update the lower triangular part, then do:

M.block(1, 2, 3, 4).triangularView<Lower>() += B;

In this case, triangularView behaves as a writing mask.

edited Jan 30 '17 at 21:36

answered Jan 30 '17 at 13:09

ggael

28,425
2
65
71

Many thanks for the .noalias() reminder, but in this case, I do not really understand why the .noalias() is needed since I compute a=M*v (rather than a=M*a). I am not trying to avoid the storage of a as I really need it for future computations. However, I am trying to avoid the copy of A.middleCols(j,n).selfadjointView() and I am still puzzled why the first version works for you, as it doesn't for me (and yes the M instead of A is a typo here, not in my code). – itQ Jan 30 '17 at 13:45
This is because we cannot know at compile-time whether you are doing a=Mv or a=Ma, and even at runtime this is very difficult in the general case. So by default, a=Mv is evaluated as if there is aliasing between a and M or v, so: `tmp = M*v; a=tmp;`. – ggael Jan 30 '17 at 13:59
Thank you for your suggestion. It's almost at it, but a little thing remains, that I added as the EDIT2 section. – itQ Jan 30 '17 at 15:22
Many thanks! I was confused with the use of .selfadjointView<>(). I thought it deals with the storage as well. So my goal was indeed to fill the lower triangular part of A with that of B, and I thought I was warning Eigen that B is symmetric. – itQ Jan 31 '17 at 09:39

c++ Eigen making copy?

1 Answers1