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 885 - Extend API of TriangularView and SelfAdjointView
Summary: Extend API of TriangularView and SelfAdjointView
Status: DECISIONNEEDED
Alias: None
Product: Eigen
Classification: Unclassified
Component: Core - general (show other bugs)
Version: 3.3 (current stable)
Hardware: All All
: Normal Feature Request
Assignee: Nobody
URL:
Whiteboard:
Keywords:
Depends on:
Blocks: 3.x
  Show dependency treegraph
 
Reported: 2014-09-24 16:29 UTC by Christoph Hertzberg
Modified: 2019-02-20 14:56 UTC (History)
4 users (show)



Attachments

Description Christoph Hertzberg 2014-09-24 16:29:01 UTC
This causes things like the following not to compile:

template<class Derived>
void foo(const Eigen::MatrixBase<Derived>& ) {}

int main() {
	Eigen::MatrixXd A;
	foo(A.triangularView<Eigen::Upper>());
}

Also, these classes miss some useful functions, like isZero(), or sum().
Comment 1 Gael Guennebaud 2014-09-29 13:48:22 UTC
I'm not convinced that Triangular/SelfAdjoint-View should inherit MatrixBase as MatrixBase is really for "dense with full storage" objects, but I agree that we should support more features in TriangularBase.
Comment 2 Christoph Hertzberg 2014-10-01 18:05:38 UTC
Hm, then it does not make much sense to me that, e.g., the following lines work:

  Eigen::MatrixXd A(3,3);
  foo(Eigen::MatrixXd::Identity(3,3));
  foo(A.triangularView<Eigen::Upper>()*Eigen::MatrixXd::Identity(3,3));
  //foo(A.triangularView<Eigen::Upper>()); // does not work

So far, I always considered MatrixBase as any kind of dense matrix, including Matrix expressions.
Comment 3 Gael Guennebaud 2014-10-06 12:02:12 UTC
(In reply to Christoph Hertzberg from comment #2)
> Hm, then it does not make much sense to me that, e.g., the following lines
> work:
> 
>   Eigen::MatrixXd A(3,3);
>   foo(Eigen::MatrixXd::Identity(3,3));
>   foo(A.triangularView<Eigen::Upper>()*Eigen::MatrixXd::Identity(3,3));
>   //foo(A.triangularView<Eigen::Upper>()); // does not work

Here, MatrixXd::Identity is explicitly a dense MatrixXd, but you build it a DiagonalMatrix, then none of the calls to foo will work.

> So far, I always considered MatrixBase as any kind of dense matrix,
> including Matrix expressions.

True, but can we consider a TriangularView as a dense matrix? In the extreme case we also have SparseView, for which it is pretty clear that it should inherit SparseMatrixBase, and not MatrixBase. With that reasoning, I agree that SelfAdjointView is more debatable.
Comment 4 Christoph Hertzberg 2014-10-07 07:33:12 UTC
The return type of SparseMatrixBase::triangularView inherits from SparseMatrixBase already, (last time I checked).
I admit that things like DenseMatrix::triangularView<Upper>().triangularView<Lower>() are indeed debatable, but API-wise I would keep them 'Dense' until explicitly using .sparseView().

Motivation, why someone would like to use MatrixBase here instead of EigenBase is that MatrixBase::eval()/MatrixBase::nested<2>() will (should) always result in Direct-Access-types, which is not the case for Sparse matrices.
Comment 5 Christoph Hertzberg 2015-05-01 20:09:51 UTC
Another problem is that the output stream operator does not work for triangularView (although it used to work with 3.2)
Comment 6 Gael Guennebaud 2015-05-04 12:39:28 UTC
To comment #5:

The following does not work for me, neither with 3.1, nor 3.2:

int main() {
  Eigen::MatrixXd A;
  std::cout << A.triangularView<Eigen::Upper>() << "\n";
}



To comment #4:

A counter argument is that you cannot access to a "dense" TriangularView as a general dense matrix: accesses to half of the coefficients are forbidden.

I had for a long time the idea to introduce the concept of matrices defined by "columns (resp. rows) for which a unique dense sub-segment is non-zero". Here, "dense" really means that the respective sub-segment inherits MatrixBase with efficient random coeff/packet access. This way we could unify dense matrices, triangular-view, band matrices, etc.
Comment 7 Christoph Hertzberg 2015-05-04 13:33:12 UTC
(In reply to Gael Guennebaud from comment #6)
> The following does not work for me, neither with 3.1, nor 3.2: [...]
>   std::cout << A.triangularView<Eigen::Upper>() << "\n";

Yes, sorry I did not actually try this. The equivalent used to work with Eigen2.0 though (I noticed that while porting deprecated doxygen code snippets).
And apparently, this is hardly ever required (I can't remember any discussion about that on bz/ML/IRC).

> To comment #4:
> 
> A counter argument is that you cannot access to a "dense" TriangularView as
> a general dense matrix: accesses to half of the coefficients are forbidden.

Yes, that would still require to specializes stream operators for TriangularView. I guess we already had the discussion about letting coeff return Scalar(0) for the other half? (Unfortunately, this is not easily possible with coeffRef, though -- and introduces overhead ...)

> I had for a long time the idea to introduce the concept of matrices defined
> by "columns (resp. rows) for which a unique dense sub-segment is non-zero".
> Here, "dense" really means that the respective sub-segment inherits
> MatrixBase with efficient random coeff/packet access. This way we could
> unify dense matrices, triangular-view, band matrices, etc.

Sounds interesting -- this would also allow, e.g., accessing sub-blocks of triangular matrices.
Comment 8 Daniel Vollmer 2015-06-23 10:38:46 UTC
I'd also just realised that triangularView doesn't have the "usual" methods (in my case I wanted isZero() for checking the type of a Butcher-tableau). I think it would be nice and consistent (IMO) to have them.
Comment 9 Gael Guennebaud 2015-09-09 09:24:01 UTC
Alright, so I changed the summary of this entry to reflect that what is most needed is to extend the feature set supported by triangular and selfadjoint views. This include:

- is*()
- operator<<
- set*/fill
- sum(), prod(), max/minCoeff()
- *norm()
- row/column-wise operations

Of course all the reductions and checks will work on the triangular part only. Otherwise, prod() would always return 0, isOnes() would be always false, etc. (except for a 1x1 matrices...).

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