| Re: [eigen] Matrix product between some MatrixType and triangularView<UnitLower> |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
Thanks !
When looking at the new code, I have some questions though :
- Is .col(i).head(i) better than .block(0,i,i-1,1) for performance, or is
it just for readability purpose ? I was using block() because I was sure it
was doing computation inplace (and because I'm not used to this syntax too)
- V.const_cast_derived() is the same as const_cast<...>(V) ?
- Is there any difference between .transpose() and .adjoint for double
real matrix ? I didn't try to get the routine work for complex values, so I
can't ensure that BlockHouseholderTriangularFactor is valid for complex
value yet.
- What coding standard should I change if I want this file to be included
in the experimental/ directory of Eigen ?
Thanks, Vincent
On Tue, 15 Jun 2010 22:03:01 +0200, Gael Guennebaud
<gael.guennebaud@xxxxxxxxx> wrote:
> ok, the TriangularView class has been extended to trapezoidal matrices
> but not the underlying matrix products whence your issue. Problem
> solved. I've also attached a modified version of your file which works
> here.
>
> gael
>
> On Tue, Jun 15, 2010 at 5:01 PM, <vincent.lejeune@xxxxxxxxxx> wrote:
>> &On Tue, 15 Jun 2010 16:29:12 +0200, Gael Guennebaud
>>
>> <gael.guennebaud@xxxxxxxxx> wrote:
>>
>>> Actually, the pb is line 50 where you store a reference to the
>>
>>> temporary matrix returned by BlockHouseholderTriangularFactor.
>>
>>>
>>
>>> I would rather write it like this:
>>
>>>
>>
>>> const TriangularView<Derived,UnitLower>& Vtri(V);
>>
>>> TriangularFactorType T = BlockHouseholderTriangularFactor(Vtmp,
>>
>> coeffs,V);
>>
>>> A -= Vtri * T.transpose() * Vtri.transpose() * A;
>>
>>
>>
>> When using "const TriangularView<Derived,UnitLower>& Vtri(V)" the
result
>>
>> is changed :/
>>
>>
>>
>> I've joined the result of a computation on some 8x8 matrix (with a
block
>>
>> factor of 2) and,
>>
>> as you can the, the 4x4 right bottom corner is wrong...
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>>
>>
>>> But here there are still many many temporaries: one for T and one for
>>
>>> each matrix product! Also here the matrix T is triangular right ? So
>>
>>> I'd suggest the following:
>>
>>>
>>
>>> Derived::PlainObject tmp(Vtri.transpose() * A);
>>
>>> tmp.noalias() = T.transpose().triangularView<Upper>() * tmp;
>>
>>> A.noalias() -= Vtri * tmp;
>>
>>>
>>
>>> There are still 2 temporaries, one for T and for tmp, but that's still
>>
>>> better than 4. I'm also not 100% sure about the second line, because
>>
>>> in theory a triangular matrix times a matrix can be performed in-place
>>
>>> but I have to check the implementation. If you have troubles just
>>
>>> remove the .noalias().
>>
>>
>>
>> What is the purpose of noalias() ?
>>
>>
>>
>>
>>
>>>
>>
>>> But this still looks relatively expensive in term of temporaries,
>>
>>> there must be something smarter.
>>
>>>
>>
>>> gael
>>
>>>
>>
>>> On Tue, Jun 15, 2010 at 4:08 PM, <vincent.lejeune@xxxxxxxxxx> wrote:
>>
>>>> Here is my source code, the part of interest is at the top
>>
>>>>
>>
>>>> (HHUtilities<..>).
>>
>>>>
>>
>>>>
>>
>>>>
>>
>>>> The incrimined lines is the 47/48 : with line 47, everything works,
but
>>
>>>> it
>>
>>>>
>>
>>>> does a copy I would like to avoid.
>>
>>>>
>>
>>>> With line 48, I do no copy, but it crashs.
>>
>>>>
>>
>>>>
>>
>>>>
>>
>>>> "Derived" is Matrix<float,Dynamic,Dynamic> and OtherDerived is
>>
>>>>
>>
>>>> Block<MatrixType> where MatrixType is the previous
>>
>>>>
>>
>>>> Matrix<float,Dynamic,Dynamic> one.
>>
>>>>
>>
>>>>
>>
>>>>
>>
>>>> On Tue, 15 Jun 2010 16:02:04 +0200, Gael Guennebaud
>>
>>>>
>>
>>>> <gael.guennebaud@xxxxxxxxx> wrote:
>>
>>>>
>>
>>>>> Hi,
>>
>>>>
>>
>>>>>
>>
>>>>
>>
>>>>> could you resend your piece of code in a text file please, because
>>
>>>>
>>
>>>>> inlined in the email, it is really hard to read !
>>
>>>>
>>
>>>>>
>>
>>>>
>>
>>>>> Regarding the previous issues, what was it ?
>>
>>>>
>>
>>>>>
>>
>>>>
>>
>>>>>
>>
>>>>
>>
>>>>> cheers,
>>
>>>>
>>
>>>>> gael
>>
>>>>
>>
>>>>>
>>
>>>>
>>
>>>>>
>>
>>>>
>>
>>>>> On Tue, Jun 15, 2010 at 3:42 PM, <vincent.lejeune@xxxxxxxxxx>
wrote:
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> Hi,
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> I fixed the issue I had that prevent me from getting a working
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> implementation of QR Householder decomposition, and apply
Householder
>>
>>>>
>>
>>>> on
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> the left/right function.
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> However, Im' currently facing another issue : my implementation is
>>
>>>>
>>
>>>> slower
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> than the unblocked version of theses functions, because I'm relying
>>
>> on
>>
>>>>
>>
>>>>>> the
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> "triangularView" template.
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> As I understand, matrix.triangularView() does not copy the matrix
>>
>>>>
>>
>>>> object,
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> and provides some comfort. However, the following
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> applyBlockHouseholderOnTheLeft function does trigger a segfault :
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> template<typename Derived, typename OtherDerived = Matrix<
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> typename Derived::Scalar, Derived::ColsAtCompileTime, Dynamic>
>>>>>> >
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> class HHutilities
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> {
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> public:
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> typedef Matrix<typename Derived::Scalar,
>>
>>>>
>>
>>>>>> Derived::ColsAtCompileTime,
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> Derived::ColsAtCompileTime> TriangularFactorType;
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> typedef typename ei_plain_diag_type<Derived>::type CoeffVect;
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> TriangularFactorType
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> static
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> BlockHouseholderTriangularFactor(const Derived& V,
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> const CoeffVect& coeffs,const Derived& VV)
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> {
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> const int k = V.cols();
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> ei_assert(V.cols()==k && V.rows()>=k && "Vector storing
>>>>>> Matrix
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> must have same columns than coeffs size and more rows than
columns");
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> TriangularFactorType ret = TriangularFactorType::Zero(k, k);
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> const int j = V.rows();
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> for (int i = 0; i < k; i++)
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> {
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> ret.block(0, i, i, 1) = - coeffs(i)
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> * V.block(i, 0, j - i, i).transpose() * V.block(i,
i,
>>
>> j
>>
>>>>
>>
>>>> -
>>
>>>>
>>
>>>>>>
>>
>> <script type="text/javascript"
>>
>>
src="https://webmail.scilab.org/roundcube/program/js/tiny_mce/themes/advanced/langs/fr.js?s=1240817786";></script>gt;>
>>
>>>>>> i, 1);
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> ret.block(0, i, i, 1) = ret.block(0, 0, i, i) *
>>
>>>>>> ret.block(0,
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> i, i,
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> 1);
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> ret(i, i) = coeffs(i);
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> }
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> return ret;
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> }
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> static
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> void
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> applyBlockHouseholderOnTheLeft(OtherDerived& A, const Derived&
>>
>> V,
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> const CoeffVect& coeffs)
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> {
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> const TriangularView<Derived,UnitLower>& Vtmp = V.template
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> triangularView<UnitLower> ();
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> const TriangularFactorType& T =
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> BlockHouseholderTriangularFactor(Vtmp,
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> coeffs,V);
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> A -= Vtmp * T.transpose() * Vtmp.transpose() * A;
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> }
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> }
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> With OtherDerived=Block<MatrixType>, and
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> MatrixType=Matrix<float,Dynamic,Dynamic>.
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> I'm forced to make a copy of the TriangularView :
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> const Derived& Vtmp(V.template triangularView<UnitLower> ());
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> (instead of
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> const TriangularView<Derived,UnitLower>& Vtmp = V.template
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> triangularView<UnitLower> ()
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> )
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> to get my code to work. However, this is a bad workaround, as it
>>
>> makes
>>
>>>>
>>
>>>> an
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> additional copy.
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> The overhead is not that big (0.1 s on an old P4 3Ghz for a 256x256
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> matrix), but it should not be there, as theoricaly, blocked
>>
>> computation
>>
>>>>
>>
>>>>>> are
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> better on big matrix...So I hope to improve things if I can get rid
>>
>> of
>>
>>>>
>>
>>>>>> this
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> avoidable matrix copy.
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> Thanks,
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>> Vincent
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>
>>
>>>>
>>
>>>>>>