Re: [eigen] patch for tridiagonal matrices in eigen

["Benoit Jacob" <jacob.benoit.1@xxxxxxxxx>]

Hi Mauro,

Thanks for the update.
I'll let Gael keep commenting on the technical points as he does that
much better than me, but here's my take on the release-management
part:

1) The Upper ---> UpperTriangular move should happen as soon as
possible in order to let people some time to adjust to this API
change. So if you can send me a patch with just that, i'd like to
include it in beta3 this weekend. Yes it affects a lot of places
indeed.

2) About whether to add this to Core or add a new "experimental" area,
i'd suggest instead to just wait until we branch 2.0, which should
happen in januay, in less than 1 month, at the same time that we
release eigen-2.0-rc1. Then trunk will be wide open for invasive
commits like this. It's true that currently we're in heavy pre-2.0
bugfixing mode and it's not the best time for such a change until we
branch.

Until now a lot of people use trunk to build trunk/KDE but that's
going to change, see the "kdesupport rearrangement" thread on
release-team@kde. So i have good hope that after the branching we can
freely break things in trunk.

Cheers,
Benoit

2008/12/19 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
>>>
>>> ---
>>>
>>
>>
>
> ---
>
>

---