def householder(A): # computes a QR factorization upper triangular matrix R
"""matrices W, R pour la factorisation QR'
Usage:
W, R = householder(A)
Argument:
A: matrice mxn avec m>=n
Retour:
W: matrice de vecteurs de reflexions de householder
R: matrice R de la factorisation QR
Tire de: LN Trefethen & D Bau III, Numerical Linear Algebra, 1997 SIAM Philadelphia
"""
m,n = A.shape
ident = eye(m)
W = zeros((m,n))
R = copy(A)
for k in range(0,n):
e1=ident[k:m,k]
x = R[k:m,k]
W[k:m,k]= ( int(x[0]>=0) - int(x[0]<0) )*linalg.norm(x,2)*e1+x
W[k:m,k]=W[k:m,k]/linalg.norm(W[k:m,k],2)
R[k:m][:,k:n]=R[k:m][:,k:n]-dot(outer(2*W[k:m,k],W[k:m,k]),R[k:m][:,k:n])
return W, R
def formQ(W):
"""Matrice Q de la factorisation QR
Usage:
Q = formQ(W)
Argument:
W: matrice mxn de reflexions obtenue de la décomposition Householder
Retour:
Q: Matrice Q de la factorisation QR
Note:
La matrice W se calcule avec W, R = householder(A)
"""
m,n = W.shape
Q=zeros((m,m))
# compute the product QI = Q, column by column
for i in range(0,m):
Q[i,i]=1
for k in range(n-1,-1,-1):
Q[k:m,i]=Q[k:m,i]-2*dot(outer(W[k:m,k],W[k:m,k]),Q[k:m,i])
return Q