Re: [eigen] notes on Householder and Givens |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] notes on Householder and Givens*From*: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>*Date*: Wed, 4 Feb 2009 19:16:19 +0100*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:in-reply-to:references :date:message-id:subject:from:to:content-type :content-transfer-encoding; bh=nS61zHlPMeqrr6weBXL+c8Mwdk8Ve94+AMhvYaSJ2Pw=; b=idT6NXdVaxQX8CAiwyxObkUdltjM5THnR0aiZ3J9DiV4DtfnmVBpcaN+kCuT5pWSmM H4dDauLKxl0ZDFjIfzNhhSj8pgPmDqb6qS7urwVo5wt+AxKzTWE/JaOGtT75mgvoJxPS 1vQzfRTl9UQodmXAy+DKkiQ6R2g+hG1e+uW98=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type:content-transfer-encoding; b=KPRT2z63UMi5GgTqwW2fyuXY08PlmBXg7hthGcInvB04qYKt2CvN8yVomsqqorBMTt +q7TKQQsXU948hnv0pXSmwTY2vHZcCBHVVJUs4zzzJrNlPTNpyCnToZPgHvBLGgr4efF i1fK3wHtep9HxhhY93/qo75WdiqjFRqKZN5h4=

hm, ok, now i can see that Golub & van Loan's code does handle correctly the case x0=0. What I still disagree with is their " if(sigma==0) ". Why does this come into play ? Seems redundant to me. The code I propose has only 1 branching. Cheers, Benoit 2009/2/4 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>: > Something just occured to me. > > If x_0 is "zero", that is muchSmallerThan(norm(x)), then there is > nothing wrong with this formula: > > v = x + norm(x) e_0 > > Indeed, no cancellation can occur since x0 is much smaller than > norm(x), and after that v0 is guaranteed to be approximately equal to > norm(x). > > So, i'm sorry to not have thought of that earlier, but it seems that > no pivoting is needed: i think that the following > ei_compute_householder will work for all nonzero vectors x even if x0 > = 0: > > void ei_compute_householder(const Vector& x of size n, Householder *h) > { > // compute natural householder vector v > Vector v = x; > if(ei_isMuchSmallerThan(x0, norm(x)) > { > v0 += norm(x); > } > else > { > v0 += ( x_0 / abs(x_0) ) norm(x); > } > > Vector v_starting_with_1 = v/v0; > > h->essential = v_starting_with_1.end(n-1); > h->beta = 2 / v_starting_with_1.squaredNorm(); > }; > > Notice that this is still very different from Golub's code since our > code takes good care of the case x0=0 while they don't. > > Cheers, > Benoit >

**Follow-Ups**:**Re: [eigen] notes on Householder and Givens***From:*Benoit Jacob

**References**:**[eigen] notes on Householder and Givens***From:*Benoit Jacob

**Re: [eigen] notes on Householder and Givens***From:*Keir Mierle

**Re: [eigen] notes on Householder and Givens***From:*Benoit Jacob

**Re: [eigen] notes on Householder and Givens***From:*Benoit Jacob

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] notes on Householder and Givens** - Next by Date:
**Re: [eigen] [RFC] Full pivoting LDLT** - Previous by thread:
**Re: [eigen] notes on Householder and Givens** - Next by thread:
**Re: [eigen] notes on Householder and Givens**

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