Hello,

I am interested in using eigen with non-default scalar types (interval arithmetic, multi-precision integers, etc). I haven't started experimenting yet, but from reading the doc I have a few questions.

* what operations are available if my scalar type is an integer type (no division, or at least not with the semantics eigen expects)? I guess multiplying matrices should be fine, but most reductions seem to contain divisions in their code, so I don't even know if I can get the determinant (except possibly for small sizes where explicit formulas are used). Similarly, in this context, computing the inverse of a matrix M means computing some other matrix N (something like the cofactor matrix) and a number d such that MN=dI, preferably with d not too large. And solving Ax=b also means finding (x,d) such that Ax=db. Using a fraction scalar type is different.

* my scalar type may occasionally throw an exception when doing comparisons like x>y or x==y, is it safe with eigen?

* NumTraits wants an estimate of the time it takes to perform additions and multiplications. For multiprecision arithmetic, I guess I should just put random things like 40 and 800?

* I am usually only interested in the sign of a determinant, not its actual value. For nice types like double, I understand computing the value is the fastest way to compute the sign. For other types, sometimes it would be nice to just multiply the signs of the diagonal values of the matrix, I guess I'll have to do that myself if I can get a reduction.

* if my scalar type also uses expression templates, is there something I should do to help it interact well with eigen?

* no optimization has been done for scalars in Z/2Z or in a boolean algebra?

Thank you in advance for any insight.

Hi,

marc1842fr wrote:Hello,

I am interested in using eigen with non-default scalar types (interval arithmetic, multi-precision integers, etc). I haven't started experimenting yet, but from reading the doc I have a few questions.

* what operations are available if my scalar type is an integer type (no division, or at least not with the semantics eigen expects)? I guess multiplying matrices should be fine, but most reductions seem to contain divisions in their code, so I don't even know if I can get the determinant (except possibly for small sizes where explicit formulas are used). Similarly, in this context, computing the inverse of a matrix M means computing some other matrix N (something like the cofactor matrix) and a number d such that MN=dI, preferably with d not too large. And solving Ax=b also means finding (x,d) such that Ax=db. Using a fraction scalar type is different.

Eigen does not forbid you to call determinant(), inverse() etc. on a matrix with an integral type, but of course the result is unlikely to be meaningful. I'm also not sure the divisions are performed at the right place to reduce underflow, i.e., at the very last moment, but even such a strategy would not be ideal since you are then likely to overflow... I don't know what are your objectives, but it sounds odd to me to want to do linear algebra with integral types.

* my scalar type may occasionally throw an exception when doing comparisons like x>y or x==y, is it safe with eigen?

Eigen won't catch your exceptions and it is up to you to deal with them.

* NumTraits wants an estimate of the time it takes to perform additions and multiplications. For multiprecision arithmetic, I guess I should just put random things like 40 and 800?

These values are mainly used to control unrolling, so putting an arbitrarily large value will simply disable unrolling for this scalar type that is probably a good idea for multiprec arithmetic.

* I am usually only interested in the sign of a determinant, not its actual value. For nice types like double, I understand computing the value is the fastest way to compute the sign. For other types, sometimes it would be nice to just multiply the signs of the diagonal values of the matrix, I guess I'll have to do that myself if I can get a reduction.

I see we are lacking a coefficient wise sign() function. Once it is added, you could simply do:

mat.diagonal().array().sign().prod();

Without a sign() function, currently you can also do:

(mat.diagonal().array()<0).count() % 2

* if my scalar type also uses expression templates, is there something I should do to help it interact well with eigen?

This is still an open issue that we definitely want to address in the future.

* no optimization has been done for scalars in Z/2Z or in a boolean algebra?

indeed, nothing special yet.

Thank you for you answer.

ggael wrote:I don't know what are your objectives, but it sounds odd to me to want to do linear algebra with integral types.

I want to do exact computations (yes, I'll use some filters so I don't do the slow exact computation all the time). That means double is out of the question. Even multi-precision floating-point numbers won't work with division. The obvious solution is to use fractions (pairs of multi-precision integers).

When you work with a point, instead of storing d fractions, it may be more interesting to store d+1 integers that represent homogeneous coordinates (in some sense, the last integer is a common denominator). With this representation, if I want to know the orientation of a (d+1)-tuple of points, I need to compute the sign of the determinant of a matrix with integer coefficients. I also need some inverse/solve for other operations. Now I indeed have no idea if there are more efficient techniques than providing eigen with a matrix where all the coefficients have been converted to fractions, but it feels worth asking.

* my scalar type may occasionally throw an exception when doing comparisons like x>y or x==y, is it safe with eigen?

Eigen won't catch your exceptions and it is up to you to deal with them.

Ok, I just wanted to know if there was a risk eigen would leak memory or do some other bad thing (which can easily happen when exceptions are thrown at unexpected places). Obviously I want to catch my exceptions myself

I see we are lacking a coefficient wise sign() function. Once it is added, you could simply do:

mat.diagonal().array().sign().prod();

Without a sign() function, currently you can also do:

(mat.diagonal().array()<0).count() % 2

Nice. Well I'd have to do the second version twice or some additional test on rank (sign is {-1,0,1} to me) but it looks ok.

* if my scalar type also uses expression templates, is there something I should do to help it interact well with eigen?

This is still an open issue that we definitely want to address in the future.

So it will work, just not as fast as possible (I don't mind that much, for now).

marc1842fr wrote:Thank you for you answer.
When you work with a point, instead of storing d fractions, it may be more interesting to store d+1 integers that represent homogeneous coordinates (in some sense, the last integer is a common denominator). With this representation, if I want to know the orientation of a (d+1)-tuple of points, I need to compute the sign of the determinant of a matrix with integer coefficients. I also need some inverse/solve for other operations. Now I indeed have no idea if there are more efficient techniques than providing eigen with a matrix where all the coefficients have been converted to fractions, but it feels worth asking.

Ok, I see. Indeed, computing the sign of the determinant is no big deal. For solving, it is possible to perform a division free LU factorization of the form:
PA = L D^-1 U
where P is the classical permutation and D^-1 is a diagonal matrix (all matrices are integer matrices).
Then I guess you could use such a decomposition to solve for Ax=b where x and b are in the homogenous form you described, but all of this requires special versions of the respective algorithms. It's definitely much simpler to use a rational scalar type.

Ok, I just wanted to know if there was a risk eigen would leak memory or do some other bad thing (which can easily happen when exceptions are thrown at unexpected places). Obviously I want to catch my exceptions myself

ah, ok sorry, I did not quite understood your concern, then yes memory leaks could happen if for instance an exception is triggered by a multiplication or an addition during a matrix product, or more generally anytime a temporary is created by an expression.

I see we are lacking a coefficient wise sign() function. Once it is added, you could simply do:

mat.diagonal().array().sign().prod();

Without a sign() function, currently you can also do:

(mat.diagonal().array()<0).count() % 2

Nice. Well I'd have to do the second version twice or some additional test on rank (sign is {-1,0,1} to me) but it looks ok.

adding such a sign function to the API is straightforward and can even be done outside Eigen source code.

News: Eigen should now be safe regarding custom scalar types and exceptions.