Re: [eigen] Why not Matrix<T, const int rows, const int cols> ??

[Benoît Jacob <jacob@xxxxxxxxxxxxxxx>]

Le vendredi 05 janvier 2007 13:59, Andre Krause a écrit :
> Why is EIGEN restricted to square matrices?

Ah, very good question! At the beginning, Eigen had rectangular matrices, with 
template arguments Width and Height as you suggest. Then it created 
complications with the multiplication operators. Eigen does not only provide 
fixed-size matrices, it also provides dynamic-size matrices (class MatrixX) 
and both inherit the same MatrixBase using a 'curiously recurring template 
pattern'. The problem is that multiplication of rectangular matrices works 
differently for fixed-size matrices than it does for dynamic-size matrices. 
Indeed, for fixed-size matrices this operator takes an additional template 
parameter, int other_width. The type of the 'other' matrix is not the same as 
the type of *this. On the other hand, with dynamic-size matrices, there is no 
need for an additional template parameter and all matrices are of the same 
type, MatrixX<T>.

The only solution is to take the matrix-matrix product out of MatrixBase and 
implement it separately in Matrix and in MatrixX. Good, that's what I did. 
But there are other operator* and multiply() methods in MatrixBase, for 
matrix-vector product. They were now hidden by the operator* and multiply() 
implemented in the child class. So I had to remove them also from MatrixBase 
and implement them also in Matrix and MatrixX separately. 

At this point I started wondering whether all that was really useful. It was 
really starting to make Eigen less elegant, fatter, and partly defeated the 
point of having this curiously recurring template pattern. The goal of Eigen 
being to cover the needs of applications, I started to ask people around me 
what matrix width and height they planned to actually use in their apps. All 
answers mentionned only square matrices. Same for me. Only one person 
mentionned rectangular matrices, that was for linear regression, but since 
then I have implemented linear regression in Eigen using only square matrices 
(by unfolding the rectangular matrix product one needs to compute at some 
point. By the way that allowed for a +50% speed optimization).

So I removed all the rectangular stuff, and was amazed at how much it 
simplified the source code, the API, the documentation, and the tests (the 
test app already weighs 2 MB, think how much it'd be with all the 
width/height combinations).

Do you have a real need for rectangular matrices?

> by the way - is there any code that shows how to implement
> pseudo-inverse calculation using eigen?

Hmm, the best way to do that is by singular values decomposition, which Eigen 
doesn't do. However, Eigen does LU decomposition with complete pivoting. I 
don't know, maybe that could is sounds precisely like something that could 
require rectangular matrices. Or better, Singular Values Decomposition, which 
Eigen doesn't have either help here. It decomposes you matrix as a product

PLUQ

where P and Q are permutation matrices, hence unitary, hence it is enough to 
compute the pseudoinverse of LU. As L and U are triangular, I guess their 
pseudoinverses are easily computed, and it remains to deduce from that the 
pseudoinverse of LU. I don't know if that's feasible.

Hope that helps,
Benoit

Attachment: pgpJl2cYp6y1p.pgp
Description: PGP signature