# Householder vector Classic List Threaded 4 messages Open this post in threaded view
|

## Householder vector

 Hi all!I am planning to implement Broyden's method (F(x)=0 for vector problems without supplied Jacobian). In Numerical Recipes it is pointed out, that Jacobian approximation update in Broyden's method can be done cheaply by updating QR-decomposition (due to the update having a special form). Currently, PolyMath offers QR-decomposition but not the update on existing decomposition. So, I've been looking into its implementation and by looks of things, it uses Householder algorithm. But then, there is a line there which calls `DhbVector>>householder`. I was trying to figure out what it does (in comments it says it returns a pair of "beta" and "householder vector"), but not with much success.PS. Householder reflection of the vector around hyperplane with normal V is defined as P = I - 2VV' (for real numbers). Then each vector X is transformed via X* = PX. In QR for Householder one chooses column X of the matrix and gets U = X - a E1, where | a | = || X ||, and sign(a) = -sing(pivot(X)) for numerical stability (but generally is not important for transformation itself) and E1 = (1, 0, 0, ..., 0)'. Then, V = U / || U || becomes the normal vector for Householder reflection. So, I'm trying to figure out what out of all this is meant by "householder vector".Thanks, Alexey -- You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout.
Open this post in threaded view
|

## Re: Householder vector

 On Tue, May 3, 2016 at 11:12 AM, Alexey Cherkaev wrote:...I've been looking into its implementation and by looks of things, it uses Householder algorithm. But then, there is a line there which calls `DhbVector>>householder`. I was trying to figure out what it does (in comments it says it returns a pair of "beta" and "householder vector"), but not with much success.becomes the normal vector for Householder reflection. So, I'm trying to figure out what out of all this is meant by "householder vector". Hi Alexey,i guess you had a look at the docu https://github.com/SergeStinckwich/SciSmalltalk/wiki/Math-DHB-wk .do you have access to the mentioned book (Golub)? unfortunately i dont remember exactly what i have done there, but DhbVector>>householder returns {b.V} with orthogonal P = I - bVV' and ||X||E1 = PX and (V at:1)=1. does that halfway answer your question?werner -- You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout.
 On 3 May 2016 at 16:09, werner kassens wrote:Hi Alexey,i guess you had a look at the docu https://github.com/SergeStinckwich/SciSmalltalk/wiki/Math-DHB-wk .do you have access to the mentioned book (Golub)? I didn't have. I've got it now from MIT website. unfortunately i dont remember exactly what i have done there, but DhbVector>>householder returns {b.V} with orthogonal P = I - bVV' and ||X||E1 = PX and (V at:1)=1. does that halfway answer your question?wernerLet's see if I understood :)Let’s say `x = #(3 4 0) asDhbVector`. The result of `householder` is: ``````#(3 4 0) asDhbVector householder. "an Array((2/5) a DhbVector(1 -2 0))" `````` So, `b = 2/5` and `v = #(1 -2 0) asDhbVector`. Then, `P = I - b v v'` must be orthogonal (PS. Is there an easier way to get identity matrix?): ``````p := (DhbMatrix rows: #(#(1 0 0) #(0 1 0) #(0 0 1))) - (2/5 * (#(1 -2 0) asDhbVector tensorProduct: #(1 -2 0) asDhbVector)). "a DhbVector((3/5) (4/5) 0) a DhbVector((4/5) (-3/5) 0) a DhbVector(0 0 1)" `````` which looks correct. And `p * x` should return `#(5 0 0) asDhbVector`: ``````p * #(3 4 0) asDhbVector. "a DhbVector(5 0 0)" `````` Worked! Thanks! ​ -- You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout.
 >(PS. Is there an easier way to get identity matrix?)DhbSymmetricMatrix>>identity: aDimensionwernerOn Tue, May 3, 2016 at 5:12 PM, Alexey Cherkaev wrote:On 3 May 2016 at 16:09, werner kassens wrote:Hi Alexey,i guess you had a look at the docu https://github.com/SergeStinckwich/SciSmalltalk/wiki/Math-DHB-wk .do you have access to the mentioned book (Golub)? I didn't have. I've got it now from MIT website. unfortunately i dont remember exactly what i have done there, but DhbVector>>householder returns {b.V} with orthogonal P = I - bVV' and ||X||E1 = PX and (V at:1)=1. does that halfway answer your question?wernerLet's see if I understood :)Let’s say `x = #(3 4 0) asDhbVector`. The result of `householder` is: ``````#(3 4 0) asDhbVector householder. "an Array((2/5) a DhbVector(1 -2 0))" `````` So, `b = 2/5` and `v = #(1 -2 0) asDhbVector`. Then, `P = I - b v v'` must be orthogonal (PS. Is there an easier way to get identity matrix?): ``````p := (DhbMatrix rows: #(#(1 0 0) #(0 1 0) #(0 0 1))) - (2/5 * (#(1 -2 0) asDhbVector tensorProduct: #(1 -2 0) asDhbVector)). "a DhbVector((3/5) (4/5) 0) a DhbVector((4/5) (-3/5) 0) a DhbVector(0 0 1)" `````` which looks correct. And `p * x` should return `#(5 0 0) asDhbVector`: ``````p * #(3 4 0) asDhbVector. "a DhbVector(5 0 0)" `````` Worked! Thanks! ​ -- You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout. -- You received this message because you are subscribed to the Google Groups "SciSmalltalk" group. To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email]. For more options, visit https://groups.google.com/d/optout.