Re: [eigen] triangular solve

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

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 <markoilcan@xxxxxxxxx> 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
>> about the solve:
>>
>> 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
>>
>>
>