Smalltalk › Frameworks & Tools › SciSmalltalk

Householder vector

_‹ Previous Topic Next Topic _›

Classic

List

Threaded

4 messages Options

Alexey Cherkaev

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.

werner kassens-2

Re: Householder vector

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

Alexey Cherkaev

Re: Householder vector

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!

werner kassens-2

Re: Householder vector

>(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.