Re: [eigen] triangular solve
• To: eigen <eigen@xxxxxxxxxxxxxxxxxxx>
• Subject: Re: [eigen] triangular solve
• From: Gael Guennebaud <gael.guennebaud@xxxxxxxxx>
• Date: Fri, 21 Dec 2012 23:31:30 +0100
• Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=20120113; h=mime-version:in-reply-to:references:from:date:message-id:subject:to :content-type; bh=vRHtAF9bm4qkaGBFZXAi/eMotaJq+ZQyng/OMcyKqrU=; b=FiA8pF/ZXhqLc2SViX/Zpl1LekMGC7f0uNom4aD7k9ocgZj7qWmNGCvXkfziI67GXK NYoVo962WmPehE9vVMLh+4Oymiu4s/g4EwdXyYYnAQGJTGghJZPB09Rl7PiYmhd4Utdw hZC5GNG4FPjZyO2jvwdPB6xPmXW3UpBDs9QVInOoUEnuW8dxQX+ug7BYfWrEon8rKDHJ iM4UthH6pCbf+gM2YsfxrHYO4K1y467q2yOhSzhKIVuIdkz9j7o3/l7PUN3Y0A7HiPB9 tnpl3SOfqCLkC4F4r46HFpzMon8Xo117sW3MJn5RGiLlzJ8fESd6Cax8AxwYgCVRo8mv JMvg==

Since your matrices are quite large, the cost of this temporary is completely negligible, you can still avoid one by doing:

MatrixXf tmp(this->matrixL());
tmp = total_cov.matrixL().solve(tmp); // in-place solving
float denom = tmp.squaredNorm();

Higher performance could be achieved by taking into account that tmp is triangular and so total_cov.matrixL().solve(tmp) is triangular too. However, Eigen is not able to do so yet, and I'm not aware of any library being able to do so.

cheers,
Gael

On Fri, Dec 21, 2012 at 4:20 PM, M A wrote:
Thanks. If my understanding is correct putting in the
MatrixXf(this->matrixL()) results in the allocation of another matrix.
Is it possible to avoid doing this? After all the lower triangle
matrix already exists in the LLT object (albeit in some representation
that I don't fully understand). I'm going to be using this function
literally zillions of times and the matrices are on the order of 1000
x 1000 or larger, so efficiency is highly critical.

Thanks,
Mark A

On Thu, Dec 20, 2012 at 4:39 PM, Gael Guennebaud
<gael.guennebaud@xxxxxxxxx> wrote:
> Hi,
>
> you cannot take advantage of the triangular nature of the matrix for the
> right hand side, so you have to copy it into a full matrix:
>
> float denom =
> total_cov.matrixL().solve(MatrixXf(this->matrixL())).squaredNorm();
>
> cheers,
> Gael.
>
>
> On Thu, Dec 20, 2012 at 10:53 PM, M A <markoilcan@xxxxxxxxx> wrote:
>>
>> I'm having trouble compiling eigen code to do a solve with two lower
>> triangular matrices. That is I want L1^{-1} * L2 where L1 and L2 both
>> come from LLT objects by way of matrixL(). Here is my class and
>> function,
>>
>> class Covariance_comp : public LLT<MatrixXf> {
>> public:
>>     float variance;
>>
>>     Covariance_comp( const MatrixXf &pdmat, float var) : LLT<MatrixXf>
>> (pdmat), variance( var ) {}
>>
>>     float update_fixedpt( const Covariance_comp &total_cov, VectorXf
>> resids ) {
>>         float numer = (this->matrixU() *
>> total_cov.solve(resids)).array().square().sum();
>>         cout << "numer = " << numer << endl;
>>         // This next line is the problem:
>>         float denom =
>> total_cov.matrixL().solve(this->matrixL()).array().square().sum();
>>         cout << "denom = " << denom << endl;
>>         variance *= numer / denom;
>>         return variance;
>>     }
>> };
>>
>> The problem comes in the line float denom = .... with g++ complaining
>>
>> error: no matching function for call to ‘Eigen::TriangularView<const
>> Eigen::Matrix<float, -0x00000000000000001, -0x00000000000000001, 0,
>> -0x00000000000000001, -0x00000000000000001>, 1u>::solve(const
>> Eigen::TriangularView<const Eigen::Matrix<float, -0x00000000000000001,
>> -0x00000000000000001, 0, -0x00000000000000001, -0x00000000000000001>,
>> 1u>) const’
>>
>> I don't understand the problem. Both total_cov and this are
>> essentially LLT objects. Is the problem using this->matrixL() as the
>> argument to solve? What is the best way to do this that takes
>> advantage of the triangular nature of the matrices? This is using
>> eigen 3.1.2.
>>
>> Thanks!
>>
>> Mark A
>>
>>
>

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