Hi guys,

just started to look at Eigen2 to do some fast linear algebra for us... what I'm missing though is the ability to address super- and subdiagonals of a matrix. This could be added in the manner of octave (matlab) by having an integer argument to diagonal() and asDiagonal() which indicates which superdiagonal (subdiagonal if negative) to return, which defaults to zero, indicating the main diagonal.

Cheers

Another one: the ability to compute the log-determinant as the sum of logs of the absolute values of the diagonal entries of an LU (or other) decomposition so as to avoid over/underflow, as warned of in determinant().

As with the super/subdiagonals, I can work around this, but it would be nice to be able to say lu.logDeterminant() instead of

lu.matrixLU().diagonal().cwise().abs().cwise().log().sum()

I am impressed though that this expression worked out of the box... but is the double invocation of cwise() efficient?

Cheers

And a quick documentation bug report: in the sample code at
http://eigen.tuxfamily.org/api/Tutorial ... orialAdvLU
the line
Eigen::LUDecomposition lu(A);
should be
Eigen::LU lu(A);

Cheers

For kind of subDiagonal/superDiagonal methods, why not, but the index will have to be a template parameter, eg:

mat.diagonal()

Is that really useful to have a "asSubDiagonalMatrix" function ?? eg:
Vector2f(1,2).asSubDiagonalMatrix() would return:

0 0 0
1 0 0
0 2 0

do you have any use case ? because from my opinion it is better/faster to use combine block expressions with asDiagonal()....

For lu.logDeterminant(), this sounds to be a good idea.

ggael wrote:do you have any use case? because from my opinion it is better/faster to use combine block expressions with asDiagonal()....

So to get the superdiagonal of m for instance I can use "m.corner(TopRight, m.rows()-1, m.cols()-1).diagonal()"? Neat - I agree that a dedicated function is only syntactic sugar then. Hadn't actually occurred to me... thanks!

My use case is actually more involved though: I need to access (like a vector) the matrix elements m(2*i, 2*i+1) for integer i. (This kind of thing happens a lot when doing linear algebra with skew-symmetric matrices.) With offsets and strides this is easy and efficient to do. Eigen's blocks essentially implement offsets and ranges, but I can't find strides, i.e., make a vector that is every k'th element of some other vector, for given k? In Matlab-ese, Eigen seems to implement index ranges of the form i:j, but not i:k:j?

linoleo wrote:

ggael wrote:do you have any use case? because from my opinion it is better/faster to use combine block expressions with asDiagonal()....

So to get the superdiagonal of m for instance I can use "m.corner(TopRight, m.rows()-1, m.cols()-1).diagonal()"? Neat - I agree that a dedicated function is only syntactic sugar then. Hadn't actually occurred to me... thanks!

no, I agree that being able to easily extract as a vector a super/sub diagonal of a matrix would be nice (and I'm willing to add it), but for the other way round, i.e., converting a vector as a matrix with only one subdiagonal (and all other coeff==0)... really I don't see the purpose.

My use case is actually more involved though: I need to access (like a vector) the matrix elements m(2*i, 2*i+1) for integer i. (This kind of thing happens a lot when doing linear algebra with skew-symmetric matrices.)

we also plane to add a good support for skew matrix with a dedicated storage, but this won't happens very soon unless someone step in !

With offsets and strides this is easy and efficient to do. Eigen's blocks essentially implement offsets and ranges, but I can't find strides, i.e., make a vector that is every k'th element of some other vector, for given k? In Matlab-ese, Eigen seems to implement index ranges of the form i:j, but not i:k:j?

ah, I see, but if we extend diagonal to support sub/super diagonal this should be enough right ? do you think generic stride objects would still be useful ?

ggael wrote:no, I agree that being able to easily extract as a vector a super/sub diagonal of a matrix would be nice (and I'm willing to add it), but for the other way round, i.e., converting a vector as a matrix with only one subdiagonal (and all other coeff==0)... really I don't see the purpose.

If I need to (say) add a vector v of values to the superdiagonal of m, as long as diagonal() returns an lvalue I don't need asDiagonal() - I can just say m.diagonal() += v. I guess the purpose of asDiagonal() is to return a special type of matrix (namely, diagonal), and agree that having special types for super/subdiagonal matrices would be overkill.

ggael wrote:we also plane to add a good support for skew matrix with a dedicated storage, but this won't happens very soon unless someone step in !

Nice - there are very few packages out there that support skew symmetry.

ggael wrote:ah, I see, but if we extend diagonal to support sub/super diagonal this should be enough right ? do you think generic stride objects would still be useful ?

Yes - I for one would still need strides because my vector only contains every second element of the subdiagonal of the matrix. A workaround for vector strides is to reshape the vector into a 2 by n/2 matrix, then pick one of the rows. I don't think you have a matlab-style "reshape" command, but I can just define a reshaped matrix with the storage of the original one.

So to get a vector v of the matrix elements m(2*i, 2*i+1) I have to say:

VectorXd x = m.corner(TopRight, m.rows()-1, m.cols()-1).diagonal();
VectorXd v = Map >(x.data()).row(1);

This is cumbersome and constructs an intermediate with twice as many elements as needed. Moreover (unless I'm missing something) the intermediate means that I can't address the elements I want of m as lvalues - all I get is copies of them. What I'd really like to be able to say instead is things like

VectorXd a;
m.diagonal().every() += a;

where every() implements the stride (here: every second element). I think you already have all the internal machinery (namely: a stride value) to implement this for vectors.

For matrices, every() ideally would work like a reduction operator, with rowwise() and colwise() to select every k'th row or column; otherwise it could just pick every k'th element of the matrix data irrespective of its layout. I think internally you'd have to split your stride value into separate strides for rows and for columns to implement this. Not sure what the performance penalty would be for that.

Best

Yes, diagonal() returns a lvalue, so no need for any "asSubDiagonal()".

Thanks for the use case. So API wise having:

vec.every(n);
mat.colwise().every(n)
mat.rowwise().every(n)

would handle all use cases, right ?

About reshape you can currently mimic it using a Map object:

Code: Select all

MatrixXd mat(100,100);
Map rmat(mat.data(),50,200);
Map matAsVec(mat.data(),mat.cols()*mat.rows());

and you can use rmat just like any matrix.

Actually, you can also do:

Code: Select all

mat.resize(50,200);
// use the reshaped mat
// restore its initial size:
mat.resize(100,100);

this is because resize does nothing on the data if the amount of required memory does not change (this works for dynamic size matrix only).

ggael wrote:Thanks for the use case. So API wise having:

vec.every(n);
mat.colwise().every(n)
mat.rowwise().every(n)

would handle all use cases, right ?

Right, that would cover it.

Thanks for the resize() trick, very dirty

linoleo wrote:Thanks for the resize() trick, very dirty

yes the right way to go is to use Map, the resize trick is not something officially supported, it is more like an interesting side effect.