RE: (HH^t)^{-1}

[Mikael Adlers <mikael at mathcore dot com>]

> -----Original Message-----
> From: Faheem Mitha [mailto:faheem@email.unc.edu] 
> Sent: den 11 oktober 2001 03:43
> To: Michael Meyer
> Cc: gsl-discuss@sources.redhat.com
> Subject: Re: (HH^t)^{-1}
> 
> 
> 
> 
> On Wed, 10 Oct 2001, Michael Meyer wrote:
> 
> > I saw your post regarding matrix inverses on
> > gsl_discuss.
> >
> > Sampling from a multinormal distribution (given the
> > mean and coavariance  matrix) usually does not involve
> > computation of inverses.
> >
> > As long as you can draw from a standard normal
> > distribution all one needs is a Cholesky factorization
> > of the covariance matrix. This is a very easy
> > algorithm one can easily write oneself (10 lines).
> >
> > I do not fully understand how (HH^t)^{-1} is related
> > to the mean of the distribution. This mean should be a
> > vector not a matrix.
> > Please provide details of your problem.
> 
> Well, to be precise...
> 
> Given a vector h=(h_1, ... h_T) of length T, then we define a 
> matrix H of
> dimension T X 2 whose kth row is (1,h_{k-1}), where h_0 = 0. 
> Then I want
> to simulate from a bivariate normal distribution with
> 
> mean \mu = (H^t H)^{-1} H^t h and variance (H^t H)_{-1}.
> 
> This can be written as
> 
> X = H^{t}Z + \mu where Z is standard normal and X has the required
> bivariate normal distribution. I don't see any way of 
> avoiding calculation
> of the inverse here.
> 
> I need to write a function which takes the value Z (the 
> standard normal)
> and returns X.
> 
> If you've any suggestions, let me know.
> 
>                                   Sincerely, Faheem Mitha.
> 

The problem of computing \mu is a least squares problem. What you
have written is the normal equations to 

 min || H\mu - h||_2

There are several way to solve this problem. One counld form the normal
equations
as you have and compute the matrix A=H^TH and solve the system A \my = H^Th.
You may look a lot of precision in this case when you form the cross-product
H^TH
(the condition number is squared). However if may work just fine.

Another way is to solve the problem by QR factorization, H = Q*R, 
||H\mu -h||_2 = ||R \mu - Q^Th||_2. To solve this you solve R \mu = Q^Tb
This method has much better error bounds and should be preferred.

The following two functions should do the trique:

Compute the QR factorization:

int gsl_linalg_QR_decomp (gsl_matrix * A, gsl_vector * tau) 
This function factorizes the M-by-N matrix A into the QR decomposition A = Q
R. On output the diagonal and upper triangular part of the input matrix
contain the matrix R. The vector tau and the columns of the lower triangular
part of the matrix A contain the Householder coefficients and Householder
vectors which encode the orthogonal matrix Q. The vector tau must be of
length k=\min(M,N). The matrix Q is related to these components by, Q = Q_k
... Q_2 Q_1 where Q_i = I - \tau_i v_i v_i^T and v_i is the Householder
vector v_i = (0,...,1,A(i+1,i),A(i+2,i),...,A(m,i)). This is the same
storage scheme as used by LAPACK. 

Solve the system:

Function: int gsl_linalg_QR_lssolve (const gsl_matrix * QR, const gsl_vector
* tau, 
          const gsl_vector * b, gsl_vector * x, gsl_vector * residual) 
This function finds the least squares solution to the overdetermined system
A x = b where the matrix A has more rows than columns. The least squares
solution minimizes the Euclidean norm of the residual, ||Ax - b||.The
routine uses the QR decomposition of A into (QR, tau) given by
gsl_linalg_QR_decomp. The solution is returned in x. The residual is
computed as a by-product and stored in residual. 


However, using the QR factorization you will not compute the (H^TH)^(-1) so
you can't 
look at the variance matrix explicitly.

Regards,
/Mikael Adlers