Re: [eigen] patch for tridiagonal matrices in eigen

["Mauro Iazzi" <mauro.iazzi@xxxxxxxxx>]

Ok, I'm sorry for the long time it took and thanks for the comments. I
address them below.

Following what Benoit said I changed Upper to UpperTriangular and now
UpperTriangular==UpperTriangularBit
Unfortunately this affects a lot of places so I think it deserves a
patch on its own. I will check again before submitting. In any case I
agree with Gael about being cautious over the patch because a lot of
things are still unclear to me and I may have missed the proper way of
doing something.

I do not understand properly the meaning of these flags. For example
part needs assumes some to be exclusive but not others (i.e.
UpperTriangular is not UpperHessenberg, but StrictlyUpper is
UpperTridiagonal...). If part() checked using == instead of & I think
it would be safe to use the Flags to mark the type of the matrix (it
already does in some points).

This would automatically allow to skip a QR step if the matrix is
already StrictlyUpper just by checking for UpperTriangular... Or
tridiagonalizing a diagonal would be O(1) just checking for a
Tridiagonal (which was why I assumed it worked like that). What do you
think?

2008/12/15 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
>
> Hi,
>
> wow that sounds really great. I have not much to say.
>
> About the shortcuts, yes Upper must remain unchanged (Upper ==
> UpperTriangularBit), and I'm also OK to rename it to UpperTriangular.
>
> About the vectorization of DiagonalProduct, actually some cases was already
> vectorized for a long time. Also note that it does not have the
> EvaluateBeforeAssign bit, so it can already works in-place. What is missing
> is the vectorization of the other cases with all the mess we discussed last
> week... Now we have a TridiagonalProduct as well, I also think that's better
> to go for a C-style function with in-place evaluation.
>
> Marco: in ei_traits<Product<TridiagonalProduct> >::Flags you should remove
> the VectorizableBit flag, we will vectorize it later using an auxiliary
> C-style function.

you are right: done

> Also, I think we need two different implementations of TridiagonalProduct:
> one when we take the tridiagonal part of a complete matrix (or when we mark
> a complete matrix as tridiagonal), and another one when the tridiagonal
> matrix relies on a sparse storage (set of 3 vectors). Currently it relies on
> dynamic branching in TridiagonalProduct::coeff(). That's ok for fixed-size
> because our meta-unrollers will allow the compiler to remove them, but
> that's pretty bad for large sizes. Actually, an easy and generic way to
> implement it would be to use a sum of 3 diagonal products, but that would
> not be cache friendly, so...
>
> I also think the TridiagonalMatrix expression should overload the diagonal()
> function (to return a vector) and provide the respective sub-diagonal access
> functions.

yes, I will do that.

> We probably also need to add a specialization of assign to optimize simple
> operations on tridiagonal matrices (sum, etc.), more like we did for
> symmetric matrices.

It would be very good. I will look into it. Making a new Matrix type
and specializing this operations should be easy to do.

> Then there is the story of tridiagonal symmetric matrices....
>
> We should also come up with a mechanism to define TriagonalMatrix with their
> own storage. For instance with DiagonalMatrix, one can currently do:
>
> typedef DiagonalMatrix<NestByValue<VectorXf> > DiagonalMatrixXf;
>
> to almost get a "real" diagonal matrix of float with dynamic size. I say
> "almost", because we cannot set its size, resize it, etc... Well one can add
> take a vector, and add a ".asDiagonal()" all the time, but, 1) that's not
> super convenient, and 2) this would be even less convenient for
> TridiagonalMatrix, so we really need to come up with something here.

I do not follow. What's wrong with DiagonalMatrix<VectorXf>?

> To conclude, I'm obviously 100% OK that you commit your stuff whenever you
> can, but I'm not in favor to add such new beta features in Core. So I would
> suggest to put the support for tridiagonal matrix in the Sparse module or
> into a new module that we could call Experimental to give us the time we
> need to refine it a bit. What do you think ?

given how trivial it was to create this patch I think it is not worth
to commit it yet. Anyone in a strong need for tridiagonals can do it
by himself or take this from the ml. in the tree it means just more
code to mantain while at present I'm really not that sure I will have
time.

> About the QR module, I think that MatrixType is already square, isn't it ?
> so I don't really follow you... About the "if", don't bother too much since
> this need to be redesigned anyway (see the other thread), but your use case
> (the matrix is already tridiagonal) will help to make the right choices...

Yes, but is particular cases it is not assignable:
SelfAdjointEigenSolver< TridiagonalMatrix<VectorXd, VectorXd, VectorXd> >
fails because it defines m_eivec as a triangular matrix. But if you
already have a tridiagonal matrix you do not want to cast it into a
square matrix before.

not counting that then the solver cannot understand that it can skip
tridiagonalization unless is checks (n-1)(n-2) coefficients, while it
could get that info from the compiler...

The SquareMatrixType is the corresponding square type with assignable
coefficients. so I think it is the correct type for m_eivec. Same idea
could apply elsewhere.

I hope that makes sense.

cheers,

mauro

>
> cheers,
> Gael.
>
>
>
> On Mon, Dec 15, 2008 at 2:59 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx>
> wrote:
>>
>> >>> -const unsigned int Upper = UpperTriangularBit;
>> >>> -const unsigned int StrictlyUpper = UpperTriangularBit | ZeroDiagBit;
>> >>> -const unsigned int Lower = LowerTriangularBit;
>> >>> -const unsigned int StrictlyLower = LowerTriangularBit | ZeroDiagBit;
>> >>> +const unsigned int Upper = UpperTriangularBit | UpperHessenbergBit;
>> >>> +const unsigned int StrictlyUpper = Upper | ZeroDiagBit;
>> >>> +const unsigned int Lower = LowerTriangularBit | LowerHessenbergBit;
>> >>> +const unsigned int StrictlyLower = Lower | ZeroDiagBit;
>> >>>
>> >>> here I don't really follow. These are parameters for part(). So unless
>> >>> you want to add Hessenberg support in part (which would be a good
>> >>> idea) no need to change that. If you want to add Hessenberg support,
>> >>> you need to make the corresponding changes in Part.h.
>> >>
>> >> I saw that the product helper checked for diagonals using the flags,
>> >> so I did the same for tridiagonals.
>> >
>> > OK right, sorry. These are always convenient shortcuts. What needs to
>> > be updated then is the comment saying that these are the parameters
>> > for part(), actually they're not used only for part().
>>
>> Wait, no, I answered too fast. We really can't do such a thing as
>>
>> +const unsigned int Upper = UpperTriangularBit | UpperHessenbergBit;
>>
>> for example how would we then in part() distinguish between upper
>> triangular and upper hessenberg...
>>
>> the real problem is that in the first place we allowed ourselves this
>> silly shortcut Upper for upper triangular. Now that there is upper
>> hessenberg too, this is getting ambiguous. I think we just rename
>> Upper to UpperTriangular and then we add UpperHessenberg alongside.
>>
>> Yes this is an incompatible API change, there's a reason why we
>> haven't released yet. KDE / KOffice is not affected, as far as I can
>> remember.
>>
>> Cheers,
>> Benoit
>>
>> ---
>>
>
>

---