New user self-registration is disabled due to spam. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Before reporting a bug, please make sure that your Eigen version is up-to-date!
Bug 602 - [Feature request] more functions for the selfadjointView class
Summary: [Feature request] more functions for the selfadjointView class
Status: NEW
Alias: None
Product: Eigen
Classification: Unclassified
Component: Core - general (show other bugs)
Version: unspecified
Hardware: All All
: Normal Unknown
Assignee: Nobody
URL:
Whiteboard:
Keywords:
Depends on:
Blocks:
 
Reported: 2013-05-18 10:41 UTC by Elisa Friebel
Modified: 2014-05-27 19:46 UTC (History)
4 users (show)



Attachments

Description Elisa Friebel 2013-05-18 10:41:02 UTC
Dear Eigen developers,

First, I want to emphasise the benefit of a class for symmetric matrices as proposed in http://eigen.tuxfamily.org/bz/show_bug.cgi?id=42. But I need not only for dense but also for sparse matrices.

Second, as a temporary alternative it would be useful to add some features for the selfadjointView class.

For example for sparse matrices the expression

  sm2.selfadjointView<Lower>() = sm1.selfadjointView<Upper>(); //copies the upper triangular part to the lower triangular

is possible while the same code with dense matrices gives compiler errors. The line

 m.triangularView()=m2;

works. It would be nice if that was available in the selfadjointView class, too.


On the other hand the selfadjointView class combined with the operator () only works for dense matrices: The line

  m.selfadjointView<Eigen::Upper>()(0,2)=11;

compiles but with a SparseMatrix sm

 sm.selfadjointView<Eigen::Upper>()(0,2)=11

it does not.

In my opinion it makes sense that the call

 m.selfadjointView<Eigen::Upper>()(1,0)

evaluates to the the entry (0,1) in m because we're looking at its selfadjoint.In the current implementation this code leads to a range check error.
Comment 1 Gael Guennebaud 2013-05-21 23:46:48 UTC
(In reply to comment #0)
> First, I want to emphasise the benefit of a class for symmetric matrices as
> proposed in http://eigen.tuxfamily.org/bz/show_bug.cgi?id=42. But I need not
> only for dense but also for sparse matrices.

Such compact representations of triangular/selfadjoint matrices are useless for sparse storages that are already compact, and the one of SparseMatrix.
 
> For example for sparse matrices the expression
> 
>   sm2.selfadjointView<Lower>() = sm1.selfadjointView<Upper>(); //copies the
> upper triangular part to the lower triangular
> 
> is possible while the same code with dense matrices gives compiler errors. 

good point.

> On the other hand the selfadjointView class combined with the operator () only
> works for dense matrices: The line
> 
>   m.selfadjointView<Eigen::Upper>()(0,2)=11;
> 
> compiles but with a SparseMatrix sm
> 
>  sm.selfadjointView<Eigen::Upper>()(0,2)=11
> 
> it does not.

That's on purpose. A(0,2) = 11 is not even allowed with A sparse.

> In my opinion it makes sense that the call
> 
>  m.selfadjointView<Eigen::Upper>()(1,0)
> 
> evaluates to the the entry (0,1) in m because we're looking at its
> selfadjoint.In the current implementation this code leads to a range check
> error.

That's also on purpose for the following two reasons:
1- it could only be read-only because for complexes we would have to conjugate the new values (would be possible with a proxy object, but that would be overkill)
2- this would kill the performance of SelfadjointView::operator(i,j)
Comment 2 Godeffroy Valet 2014-05-27 18:59:11 UTC
Hello,

I reply to this post instead of creating a new one.

For my application, I would like to be able to make arithmetic operations like A=k*B+C, where k is a scalar and A,B,C self-adjoint matrices. Only a triangular part of the derived matrix would be referenced so that no unnecessary processing is done (like twice the same).

I would also like to be able to use .block<>() on it, and it would be great if it returned a SelfAdjointView of the block in case the requested block is on the diagonal.


What do you think ? Did I miss something ? I feel that for now, it can only be used with very few specific functions or convert back to its derived type.

Thank you.
Comment 3 Christoph Hertzberg 2014-05-27 19:17:26 UTC
(In reply to comment #2)
> For my application, I would like to be able to make arithmetic operations like
> A=k*B+C, where k is a scalar and A,B,C self-adjoint matrices. Only a triangular
> part of the derived matrix would be referenced so that no unnecessary
> processing is done (like twice the same).

You can already do
 A.triangularView<Upper>() = k*B+C;
and Eigen will only access and process the upper half of each matrix (same is possible with Lower, of course).

> I would also like to be able to use .block<>() on it, and it would be great if
> it returned a SelfAdjointView of the block in case the requested block is on
> the diagonal.

For generic blocks this will not be possible, because we would need to know at compile-time that the block is on the diagonal. It could be possible for the special case of topLeftCorner and bottomRightCorner, but still we'd have to be sure that both dimensions are the same.
If you know that the blocks you are accessing are self-adjoint, simply add an .selfAdjointView() after the .block(...).

> What do you think ? Did I miss something ? I feel that for now, it can only be
> used with very few specific functions or convert back to its derived type.

I agree that there are many cases where directly declaring a SelfAdjointMatrix would be more handy than having to add selfAdjoitView all the time.
And these types should implement basic arithmetic operations and also give compile errors, if non-selfadjoint matrices are assigned.
Comment 4 Godeffroy Valet 2014-05-27 19:46:24 UTC
(In reply to comment #3) 
> You can already do
>  A.triangularView<Upper>() = k*B+C;
> and Eigen will only access and process the upper half of each matrix (same is
> possible with Lower, of course).

Yes. I do it right now. 

The issue is that I use a Matrix4f Q, and afterwards, I access Q.block<3,3>(0,0) , Qc.block<1,3>(3,0) and Q.block<3,1>(0,3). So, if the blocks were at least mapped to the corresponding values in the upper triangular part of Q (not as trivial as with a full matrix but doable), I could use them easily. 

So, for now, I have to use Q.block<3,3>(0,0).selfAdjointView() and Qc.block<3,1>(0,3).transpose() and Q.triangularView<Upper>() = k*B+C;. And all those functions could be hidden from my code.


> I agree that there are many cases where directly declaring a SelfAdjointMatrix
> would be more handy than having to add selfAdjoitView all the time.
> And these types should implement basic arithmetic operations and also give
> compile errors, if non-selfadjoint matrices are assigned.

Exactly what I thought, like the already implemented DiagonalMatrix for diagonal matrices for example.

Is there reasons not having this SelfAdjointMatrix type ? (I agree that my blocks are easy ones and that the general case would be more arguable though)

Note You need to log in before you can comment on or make changes to this bug.