Re: [eigen] triangular solve

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

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