* RE: (HH^t)^{-1}
@ 2001-12-19 13:20 Mikael Adlers
2001-12-19 13:20 ` (HH^t)^{-1} Faheem Mitha
0 siblings, 1 reply; 3+ messages in thread
From: Mikael Adlers @ 2001-12-19 13:20 UTC (permalink / raw)
To: 'Faheem Mitha', gsl-discuss
> -----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
^ permalink raw reply [flat|nested] 3+ messages in thread
* RE: (HH^t)^{-1}
2001-12-19 13:20 (HH^t)^{-1} Mikael Adlers
@ 2001-12-19 13:20 ` Faheem Mitha
0 siblings, 0 replies; 3+ messages in thread
From: Faheem Mitha @ 2001-12-19 13:20 UTC (permalink / raw)
To: Mikael Adlers; +Cc: gsl-discuss
On Thu, 11 Oct 2001, Mikael Adlers wrote:
> The problem of computing \mu is a least squares problem. What you
> have written is the normal equations to
>
> min || H\mu - h||_2
>
> 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.
Dear Mikael Adlers,
Thanks, this is a really good suggestion. I was handed these expressions
by someone else (I think they come from Bayesian computations) but I
haven't checked them myself. Like a ninny, I didn't realise /mu was a
solution to a least squares problem. This way of handling it is certainly
much better than the clumsy things I was doing. Thanks a lot for pointing
it out.
Sincerely, Faheem Mitha.
^ permalink raw reply [flat|nested] 3+ messages in thread
* Re: (HH^t)^{-1}
[not found] <20011010212301.23259.qmail@web20705.mail.yahoo.com>
@ 2001-12-19 13:20 ` Faheem Mitha
0 siblings, 0 replies; 3+ messages in thread
From: Faheem Mitha @ 2001-12-19 13:20 UTC (permalink / raw)
To: Michael Meyer; +Cc: gsl-discuss
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.
^ permalink raw reply [flat|nested] 3+ messages in thread
end of thread, other threads:[~2001-12-19 13:20 UTC | newest]
Thread overview: 3+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2001-12-19 13:20 (HH^t)^{-1} Mikael Adlers
2001-12-19 13:20 ` (HH^t)^{-1} Faheem Mitha
[not found] <20011010212301.23259.qmail@web20705.mail.yahoo.com>
2001-12-19 13:20 ` (HH^t)^{-1} Faheem Mitha
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).