Householder vector

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Householder vector

Alexey Cherkaev
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.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Householder vector

werner kassens-2
On Tue, May 3, 2016 at 11:12 AM, Alexey Cherkaev <[hidden email]> 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.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Householder vector

Alexey Cherkaev

On 3 May 2016 at 16:09, werner kassens <[hidden email]> 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?
werner

Let'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.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Householder vector

werner kassens-2
>(PS. Is there an easier way to get identity matrix?)
DhbSymmetricMatrix>>identity: aDimension
werner

On Tue, May 3, 2016 at 5:12 PM, Alexey Cherkaev <[hidden email]> wrote:

On 3 May 2016 at 16:09, werner kassens <[hidden email]> 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?
werner

Let'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.
Loading...