Re: [eigen] inverse() method for TriangularView class?
• To: eigen@xxxxxxxxxxxxxxxxxxx
• Subject: Re: [eigen] inverse() method for TriangularView class?
• From: Douglas Bates <bates@xxxxxxxxxxxxx>
• Date: Mon, 31 Oct 2011 12:10:16 -0500
• Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=mime-version:sender:in-reply-to:references:date :x-google-sender-auth:message-id:subject:from:to:content-type :content-transfer-encoding; bh=CqGBan1tbsg0DBRWW0O3nFwuJGZMkgaJm/ReF+IZhUI=; b=XP/0GsMgGlTXeSKtd+plBA157N1vbNSWiA+9tdohq1jgPjHgifJXA4zQQgMHhHze7c Csmbm/QeS8NdpJzvSb75xDoZ+QWWXpqpkW6Ll2WZToyKmsZfe6kLLZZ6VqCBFnY+bytU S169Oe/BvFlm9G8GvIqwlZypoglBFKRQ5h41k=

On Mon, Oct 31, 2011 at 11:51 AM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:
> I don't think there's a big reason not to do it, it's just that
> probably few people have felt a need for it so far. For most tasks,
> solve() is a better tool than inverse().
>
> How big is the speed difference between triangular.solve(rhs) and
> triangular*rhs ?

In my case it is not that I want to use the inverse to form solutions
of systems.  I actually want the inverse and at present am applying
solve to the identity matrix to obtain it.  My proposal for an inverse
method is simply to make my code more concise and transparent.

A common operation in statistical computing is solving a least squares
problem, $\hat{\beta}=\arg\min_{\beta}\|y - X\beta\|$ and calculating
the standard errors of the estimated coefficients.  That calculation
requires the square roots of the diagonal elements of X'X.  If you use
a QR, pivoted QR, LLT Cholesky or LDLT Cholesky decomposition to
evaluate the solution, you can do this by evaluating the inverse of R
(or L) and taking the norm of the rows (columns for L).

> If you would like this functionality to be implemented, you could file
> a bug, explain your use case, and give a simple benchmark result
> showing a performance advantage for  triangular*rhs over
> triangular.solve(rhs).
>
> Cheers,
> Benoit
>
> 2011/10/31 Douglas Bates <bates@xxxxxxxxxxxxx>:
>> The TriangularView class has solve methods but not an explicit inverse
>> method, which could be as simple as
>> solve(DenseMatrixType::Identity(cols(), cols())) (at least I think
>> that is how one would write the general form, but I'm not sure).  If
>> one wanted to get fancier there may be a slight advantage in taking
>> account of the triangularity of the returned value and solving reduced
>> linear systems for each column, as is done in the Lapack *trtri
>> subroutines.
>>
>> Has such a method been considered, or even implemented and I have
>> somehow managed to overlook it?
>>
>>
>>
>
>
>



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