Re: [eigen] inverse() method for TriangularView class?

[Douglas Bates <bates@xxxxxxxxxxxxx>]

On Mon, Oct 31, 2011 at 1:46 PM, Helmut Jarausch
<jarausch@xxxxxxxxxxxxxxxxxxx> wrote:
> On 10/31/2011 06:10:16 PM, Douglas Bates wrote:
>> 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
>
> Isn't it more stable to use the SVD of X = U S V'  (S diagonal)
> Then (X'X)^(-1)= V S^(-2) V'

You can use the SVD but, depending on the size of the model matrix, X,
the SVD can be much, much slower than using a LDLT decomposition of
X'X.  The latter method does involve calculating X'X whose condition
number is the square of the condition number of X, which is a
disadvantage.  But even if you don't want to form X'X you could use a
column-pivoted QR decomposition to get close to the speed of the
Choleksy decomposition.  And both are much faster than the SVD on
large examples.