Re: [eigen] Implementation of M.transpose()*M & M*M.transpose() |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] Implementation of M.transpose()*M & M*M.transpose()*From*: Gael Guennebaud <gael.guennebaud@xxxxxxxxx>*Date*: Sat, 20 Feb 2010 18:25:49 +0100*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:in-reply-to:references :from:date:message-id:subject:to:content-type; bh=L2UOQFQSPr7TiB5y/gxAlFvMISQU3YvVF32JYc3AI7w=; b=JzAWfR0TNaRJM/aJjm5+FATTzS2OchYZ3xFshvwr5aGAWU52IppnixWqH6L3KTFw43 lLsHcSvPKxFyK7ghUdAEOAF7Em4M/+dq1Z88rT9k5vJYEdnmryz5pDUPJObk5RcB0a8N fOgvPBsvti0cKYXY1xI07X4zXK2dfCWJY/w24=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:from:date:message-id:subject:to :content-type; b=o2y/p19P3KcJeXXrtXWDtKda20zKm7WoMsqkQ65hk9serJ5VFJrMu1vebS+qIZc0F5 SiJjYVlTkR1djOxQkHrfFjqsabvguFltjcJl1FROAtl020XMIa0rc3lODX4VvCPJHqtb THTMBiExqitz+e4mjEKTdAN94ZzhnoNbSDjfU=

On Sat, Feb 20, 2010 at 5:38 PM, Jinglei Hu <jingleihu@xxxxxxxxx> wrote:

yes but that would be terribly slow because here you do not take care of the caches, vectorization, etc.

gael

Well, by trivial I mean one can write something as follows(not using Eigen):Matrix M_Mult_MT(const Matrix& Rhs) {

// [m][n] * [n][m] --> [m][m]

register int N(Rhs.num_rows());

double v[N][N];

double** v2d(Rhs);

for (int i = 0; i < N; ++i) {

double* vi = &(v[i][0]);

double* pi = v2d[i];

for (int j = i; j < N; ++j) {

double sum = 0.0;

double* pj = v2d[j];

for (int k = 0; k < Rhs.num_cols(); ++k)

sum += pi[k] * pj[k];

vi[j] = v[j][i] = sum;

}

}

return Matrix(N, N, &(v[0][0]));

}

yes but that would be terribly slow because here you do not take care of the caches, vectorization, etc.

gael

----- Original Message -----From:Gael GuennebaudSent:2010年2月20日 16:27Subject:Re: [eigen] Implementation of M.transpose()*M & M*M.transpose()

2010/2/20 Jinglei Hu <jingleihu@xxxxxxxxx>

Hi,I'm curious how M.transpose()*M and M*M.transpose() is implemented in Eigen by using _expression_ Template. Or it's just the same as the usual multiplication of two matrices?

Currently A = M.transpose()*M involves a full matrix multiplication. If you want to compute only an half, with the devel branch you can do:

A.setZero()

A.selfadjointView<Upper>().rankUpdate(M); // A += M.transpose()*M

and you can compute A * B as follow:

A.selfadjointView<Upper>() * B

This is the preferred way if you only need to have one half (in most of the case it is enough). If you really need to have the full result then check you really really need it ;) Ok so if you really really need it then do:

A.triangularView<Lower>() = A.transpose();

Actually, it is quite easy to detect that case in the matrix product, so in the future you can expect that A = M.transpose()*M will be automatically optimized as in Atlas for instance.

gaelThe fact that the resultant matrix is symmetric can be used to save computation cost especially for large matrix. For example, the Hessian matrix, often used in nonlinear optimization/minimization, is actually obtained as Jacobian.transpose() * Jacobian. It is trivial to write a efficient function to do this by considering the symmetry.

not that trivial when you consider cache performance as we do... have a look at the file in Eigen/src/Core/products/SelfadjointProduct.h.

But anyway it is less elegant in the sense that the code does not look nice, e.g. Matrix hessian = MT_Mult_M(Jacobian).Any idea to do this efficiently by sticking to _expression_ Template? A different but somewhat similar case is the multiplication of normal matrix by diagonal matrix.

this one is already optimized.

Cheers,Jinglei Hu

**Follow-Ups**:**Re: [eigen] Implementation of M.transpose()*M & M*M.transpose()***From:*Jinglei Hu

**References**:**[eigen] Implementation of M.transpose()*M & M*M.transpose()***From:*Jinglei Hu

**Re: [eigen] Implementation of M.transpose()*M & M*M.transpose()***From:*Gael Guennebaud

**Re: [eigen] Implementation of M.transpose()*M & M*M.transpose()***From:*Jinglei Hu

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] Implementation of M.transpose()*M & M*M.transpose()** - Next by Date:
**Re: [eigen] Implementation of M.transpose()*M & M*M.transpose()** - Previous by thread:
**Re: [eigen] Implementation of M.transpose()*M & M*M.transpose()** - Next by thread:
**Re: [eigen] Implementation of M.transpose()*M & M*M.transpose()**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |