This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
RE: (HH^t)^{-1}
- To: 'Faheem Mitha' <faheem at email dot unc dot edu>, gsl-discuss at sources dot redhat dot com
- Subject: RE: (HH^t)^{-1}
- From: Mikael Adlers <mikael at mathcore dot com>
- Date: Thu, 11 Oct 2001 10:08:53 +0200
> -----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