MATH662 Homework 2, Issued 1/24/14, Due 1/30/14
Theme: αγεωµετρητ oς µηδεις εισιτω, Ageometretos medeis eisito, Let no one untrained in geometry
enter . Or, geometric insight yields elegant simple solutions of linear algebra problems.
1. Use Householder reflectors to show that det (I + xyT ) = 1 + xTy, x, y ∈ Rm.
Solution. Recall that transformations by unitary matrices Q ∈ Rm×m preserve lengths, and that
the determinant of a matrix A = [a1, , am] ∈ Rm×m gives the volume of parallelepiped whose
edges are the matrix column vectors ai. In particular, for unitary matrices, det Q = ±1, and
det (QA) = det Q det A = det A (assume columns of Q are right-ordered such that det Q = 1).
Let H be the Householder reflector that transforms e1 kxk into x, x = kxkHe1,
H =I −2
vv T
, v = x − kxke1 ,
vT v
and let w = H Ty, y = Hw. The determinant in this problem can be rewritten as
det (I + xy T ) = det (HH T + kxkHe1wTH T ) = det H det (I + kxke1w T ) det H T = det (I + kxke1w T ).
Note that
I + kxke1w T = ( (1 + kxkw1)e1 e2 + kxkw2e1
em + kxkwme1 ),
and the resulting determinant can be expanded along column 1 to give
det (I + kxke1wT ) = 1 + kxkw1 = 1 + kxkw Te1 = 1 + y THH T x = 1 + xT y.
2. Determine the Householder reflector H that maps x ∈ Cm into a multiple of an arbitrary y ∈ Cm.
Solution. The Householder reflector preserves lengths, and therefore transform x into αy with
α = kxk/ky k. The reflector is
H =I −2
vv ∗
v ∗v
with v = αy − x.
3. Prove the block multiplication rule
A B
C D
S T
U V
AS + BU AT + BV
=
,
CS + DU CT + DV
for arbitrary, dimension compatible, submatrices.
Solution. Let
X=
A B
C D
∈C
m×n
,Y =
S T
U V
n× p
∈C
,Z =
AS + BU AT + BV
CS + DU CT + DV
∈ Cm×p ,
and let xi ∈ Cm, yj ∈ Cn, zk ∈ Cm denote the column vectors of X , Y , Z. The XY multiplication
gives
XY = ( Xy1
1
Xy p ),
i.e. a juxtaposition of p linear combinations of the columns of X with coefficients in vectors yj .
The k th column vector is
zk = Xyk =
n
X
yjkx j
j =1
4. Construct Q unitary such that Q∗A − DQ∗ = R, with A, D, R ∈ Cm×m, D diagonal, R upper
triangular.
Solution. The standard QR = A factorization is built from
( a1
am ) = ( q 1

qm )
r11
0
r1 m
rm m

=
r11 q1 r12 q1 + r22 q2
Pm
j =1
r1jqm ,
by choosing q1 colinear with a1 and of unit 2-norm, q2 perpendicular to q1 and of unit 2-norm, and
so on. Generalize the approach to the requested factorization A = QDQ∗ + R or
( a1
∗
=
am ) − d11 q1 q1∗ − − dm mqmqm
r11 q1 r12 q1 + r22 q2
Pm
j =1
r1jqm .
(1)
Multiply (1) by q1∗ to obtain
( q1∗a1 q1∗a2
q1∗am ) − d11 q1∗ = ( r11 r12
r1m ).
(2)
As in standard QR factorization, choose q1 = a1/ka1k (all norms are 2-norms), so the first row of
R is
d11 q1∗
1
( r11 r12 r1m ) =
.
ka1k2 a∗1a2 a∗1am −
ka1k
ka1k
The procedure can be repeated for successive rows, e.g., multiply (1) by q2∗ to obtain
( q2∗a1 q2∗a2
q2∗am ) − d22 q2∗ = ( 0 r22
r2m ),
and define q2 as
q2 =
a2 − r12 q1
.
ka2 − q1∗a2k
Computational problem. Implement the Householder reduction to upper triangular form. Compare
the accuracy/performance with classical and modified Gram-Schmidt on the Hilbert matrices of
size m = 2p, p = 4, 5, 6, 7.
Algorithm - Householder to factor QR = A, A ∈ Cm×n , Q ∈ Cm×n , R ∈ Cn×n
R = A, Q = I
for k = 1 to n
x = Ak:m,k
v = sign(x1) kxk2e1 + x
v = v/kv k2
Rk:m,k:n = Ak:m,k:n − 2v(v ∗Ak:m,k:n)
Qk:m,1:n = Qk:m,1:n − 2v(v ∗ Qk:m,1:n)
Here’s a Python implementation
Python] from pylab import *
Python]
2
Python] def qrH(A):
m,n=shape(A); Q=eye(m); R=zeros([m,n])
R=A.copy(); v=zeros(m); x=zeros(m)
for k in range(n):
x[k:m]=R[k:m,k]; xnorm=norm(x[k:m],2)
v[k:m]=x[k:m]; v[k]=sign(x[k])*xnorm+x[k]; vnorm=norm(v[k:m],2)
v[k:m]=v[k:m]/vnorm
R[k:m,k:n]=R[k:m,k:n]-2*outer(v[k:m],dot(v[k:m].T,R[k:m,k:n]))
Q[k:m,:]=Q[k:m,:]-2*outer(v[k:m],dot(v[k:m].T,Q[k:m,:]))
return Q.T,R
Python]
Test the implementation on a small example
Python] from scipy.linalg import hilbert; H=hilbert(3); Qh,Rh=qrH(H)
Python] norm(H-dot(Qh,Rh),2)
0.000000
Python]
Implementation seems to be correct. Here are the classical, modified Gram-Schmidt from Lecture 5.
Python] def clgs(A):
m,n=shape(A); Q=zeros([m,n]); R=zeros([n,n]); v=zeros(m)
for j in range(n):
v=A[:,j]
for i in range(0,j-1):
R[i,j]=dot(Q[:,i],v)
v=v-R[i,j]*Q[:,i]
R[j,j]=sqrt(dot(v,v))
Q[:,j]=v/R[j,j]
return Q,R
Python] def mgs(A):
m,n=shape(A); Q=zeros([m,n]); R=zeros([n,n]); V=zeros([m,n])
V=A.copy()
for i in range(n):
R[i,i]=sqrt(dot(V[:,i],V[:,i]))
Q[:,i]=V[:,i]/R[i,i]
for j in range(i+1,m):
R[i,j]=dot(Q[:,i],V[:,j])
V[:,j]=V[:,j]-R[i,j]*Q[:,i]
return Q,R
Python]
Carry out the requested Hilbert matrix experiment
Python] for p in range(4,8):
H=hilbert(2**p)
Qc,Rc=clgs(H); Qm,Rm=mgs(H); Qh,Rh=qrH(H)
eQRc=norm(dot(Qc,Rc)-H,2); eQRm=norm(dot(Qm,Rm)-H,2)
eQRh=norm(dot(Qh,Rh)-H,2)
print ’p=’,p
print ’eQRc=’,eQRc,’ errm=’,eQRm,’ errh=’,eQRh
eQhc=norm(abs(Qh)-abs(Qc),2); eQhm=norm(abs(Qh)-abs(Qm),2)
print ’eQhc=’,eQhc,’ eQhm=’,eQhm
3
p= 4
eQRc=
eQhc=
p= 5
eQRc=
eQhc=
p= 6
eQRc=
eQhc=
p= 7
eQRc=
eQhc=
8.88476820789e-17 errm= 6.87850905907e-17
2.89454825087 eQhm= 1.34802460378
errh= 4.9230073005e-16
1.55161435469e-16 errm= 1.45296671094e-16
4.60920686037 eQhm= 1.6483199255
errh= 1.74880651206e-15
1.39368726555e-16 errm= 1.9599011519e-16
7.12469439585 eQhm= 1.74098708185
9.41341148897e-17 errm= 2.00374016764e-16
10.8092992822 eQhm= 1.79449074604
errh= 1.41699302751e-15
errh= 2.21445751118e-15
Python]
Result analysis:
•
All 3 algorithms exhibit the backward stability of Q R factorization, i.e. errors in the
reconstructed QR, kH − QRk = O(ǫmach)
•
Errors in Q are largest between classical and Householder QR factorizations
4

Scuola Categoria Tempo Tot Punteggio Totale Sede

Cross Reference

MIN-YOUNG LEE, PH.D. Curriculum Vitae EDUCATION

クイズ①