Étant donné une matrice triangulaire supérieure $R$ et un vecteur $b$, on peut facilement résoudre le système $R x = b$. La matrice triangulaire est inversible si et seulement si les coefficients de sa diagonale sont non-nuls. L'algorithme backsubs implémente la susbstitution inversée pour calculer $x = R^{-1}b$. Ici, on ne calcule pas explicitement $R^{-1}$, seulement son produit par $b$.
La factorisation QR est une factorisation d'une matrice $A$ de taille mxn de la forme $A = QR$, où $Q$ est une matrice unitaire et $R$ est une matrice triangulaire supérieure. On cherche donc à exprimer $A$ dans une base des vecteurs colonnes de $Q$ de telle sorte que $Q^{-1} A$ (la matrice $A$ dans la base $Q$), soit triangulaire.
from numpy import *
%matplotlib inline
import matplotlib.pyplot as plt
import pylab as pl
def backsubs(R,b):
"""Solution d'un système linéaire - méthode substitution inverse
Usage:
x = backsubs(R,b)
Arguments:
R: matrice carrée mxm triangulaire supérieure inversible
b: vecteur mx1
Retour:
x: solution de Rx = b
"""
m, n = R.shape
x = zeros((m,1))
for j in range(m-1,-1,-1):
xr = dot(R[j,(j+1):m],x[(j+1):m])
x[j] = (b[j]-xr)/R[j,j]
return x
Le premier algorithme, clgs
, est la méthode de Gram-Schmidt classique, qui est instable numériquement.
def clgs(A):
"""factorisation QR Gram-Schmidt classique - instable
Usage:
Q, R = clgs(A)
Argument:
A: matrice mxn avec m>=n
Retour:
Q: matrice Q de la factorisation QR
R: matrice R de la factorisation QR
Note:
Cette factorisation n'est pas stable numériquement
Tire de: LN Trefethen & D Bau III, Numerical Linear Algebra, 1997 SIAM Philadelphia
"""
m,n = A.shape
v = zeros((m,n))
R = zeros((n,n))
Q = zeros((m,n))
for j in range(0,n):
v[:,j] = A[:,j]
for i in range(0,j):
R[i,j] = dot( (Q[:,i].conjugate()).T, A[:,j] )
v[:,j] = v[:,j] - R[i,j]*Q[:,i]
R[j,j] = linalg.norm(v[:,j],2)
Q[:,j] = v[:,j]/R[j,j]
return (Q,R)
Cette méthode, mgs
, est stable numériquement.
def mgs(A):
"""factorisation QR Gram-Schmidt modifiée - stable
Usage:
Q, R = mgs(A)
Argument:
A: matrice mxn avec m>=n
Retour:
Q: matrice Q de la factorisation QR
R: matrice R de la factorisation QR
Note:
Cette factorisation est stable numériquement
Tire de: LN Trefethen & D Bau III, Numerical Linear Algebra, 1997 SIAM Philadelphia
"""
m,n = A.shape
v = zeros((m,n))
R = zeros((n,n))
Q = zeros((m,n))
for i in range(0,n):
v[:,i] = A[:,i]
for i in range(0,n):
R[i,i] = linalg.norm(v[:,i],2)
Q[:,i] = v[:,i]/R[i,i]
for j in range(i,n):
R[i,j] = dot( (Q[:,i].conjugate()).T, v[:,j])
v[:,j] = v[:,j] - R[i,j]*Q[:,i]
return (Q,R)
Cette méthode se divise en 2 étapes. Le premier algorithme, householder
, calcule une matrice $W$ et une matrice $R$. La matrice $R$ est la matrice de la factorisation QR. L'algorithme householder
ne calcule pas directement la matrice unitaire $Q$. La matrice $W$ est une matrice de vecteurs de reflexions de Householder.
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
L'algorithme formQ
peut être utilisée pour calculer $Q$, mais dans plusieurs applications, il n'est pas nécessaire de calculer explicitement $Q$, car on peut calculer $Q^*b$ à partir de $W$.
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
Le calcul du produit $Q^*b$ intervient dans la résolution d'un système linéaire $Ax = b$. En factorisant $A = QR$, le système linéaire devient $QRx =b$, et $Rx = Q^{-1}b$. Comme $Q$ est unitaire $Q^{-1} = Q^*$ et le problème revient à résoudre un système triangulaire $Rx = Q^*b$.
def Qadjb(W,b):
"""Produit matrice-vecteur Q adjoint fois b
Usage:
qb = Qadjb(W,b)
Arguments:
W: matrice mxn de reflexions obtenue de la décomposition Householder
b: vecteur mx1
Retour:
Produit Q adjoint par b, où Q est la matrice de la factorisation QR
Note:
La matrice W se calcule avec W, R = householder(A)
Tire de: LN Trefethen & D Bau III, Numerical Linear Algebra, 1997 SIAM Philadelphia
"""
m,n = W.shape
qb = copy(b)
for k in range(0,n):
qb[k:m]=qb[k:m]-2*dot(outer(W[k:m,k],W[k:m,k]),qb[k:m])
return qb
Pour tester les algorithmes, on génère une matrice (carrée) $A$ 5x5, un vecteur colonne $b$ de taille 5 et une matrice $B$ de taille 5x3.
A = random.randn(5,5)
print A
b = random.randn(5,1)
print b
B = random.randn(5,3)
print B
[[ 1.69341731 0.69711409 1.04559261 -1.19192191 -0.20859709] [-0.40558211 -0.3911715 0.20974031 0.59123117 0.9950144 ] [-0.7304288 0.4336937 0.40317186 -0.13578827 -2.19259613] [-1.76333446 -0.69704146 0.85621643 0.57051842 -1.68546467] [ 1.55064814 0.40427411 -1.40579341 1.97408806 -1.12813358]] [[-1.93928333] [ 0.88296536] [ 0.61127156] [-0.0710438 ] [-0.01844664]] [[-1.50835853 1.32661551 0.40371587] [ 1.9044091 1.00369094 -0.65311313] [ 0.28752332 -0.19939505 2.07919572] [ 1.11977711 0.2450445 0.32092334] [ 0.40043858 1.17794625 1.10765667]]
Les étapes pour résoudre le système $Ax = b$ à l'aide de la factorisation QR et la méthode de Householder sont les suivantes:
householder
.Qadjb
.backsubs
.La reconstruction de la matrice $Q$ à l'aide de formQ
n'est pas nécessaire pour résoudre le système linéaire.
W, R = householder(A)
qb = Qadjb(W,b)
x = backsubs(R,qb)
print dot(A,x)
print b
[[-1.93928333] [ 0.88296536] [ 0.61127156] [-0.0710438 ] [-0.01844664]] [[-1.93928333] [ 0.88296536] [ 0.61127156] [-0.0710438 ] [-0.01844664]]
La factorisation se fait en deux étapes:
householder
fromQ
W, R = householder(B)
Q = formQ(W)
print B
print dot(Q,R)
[[-1.50835853 1.32661551 0.40371587] [ 1.9044091 1.00369094 -0.65311313] [ 0.28752332 -0.19939505 2.07919572] [ 1.11977711 0.2450445 0.32092334] [ 0.40043858 1.17794625 1.10765667]] [[-1.50835853 1.32661551 0.40371587] [ 1.9044091 1.00369094 -0.65311313] [ 0.28752332 -0.19939505 2.07919572] [ 1.11977711 0.2450445 0.32092334] [ 0.40043858 1.17794625 1.10765667]]
On génère un jeux de données, deux vecteurs $x$ et $y$ de taille $N = 500$. Afin de bien voir la régression, on construit $y$ comme un fonction linéaire plus un bruit gaussien: $y = ax + b + \epsilon$, où $\epsilon$ est un vecteur de taille $N$ de variables normales indépendantes, identiquement distribuées de moyennes nulles et de variance $\sigma^2 = 0.1$, $a = -0.2$ et $b = 0.5$.
N = 500
sigma2 = 0.1
a = -0.2
b = 0.5
epsilon = sigma2*random.randn(N)
x = random.randn(N)
y = a*x + b + epsilon
pl.plot(x,y,'b.');
pl.xlabel('x');
pl.ylabel('y');
Le problème de régression linéaire est $Ac = y$, où $c = (a,b)$ et $A = [x, 1]$, une matrice Nx2 avec colonnes $x$ et $1$.
# construction de la matrice A pour la régression linéaire
A = ones((N,2))
A[:,0] = x
La solution $x = A^+y$ se calcule à partir des équations normales $A^*Ax = A^* y$. On a $A^*Ax = (\hat Q \hat R)^* (\hat Q \hat R)x = R^* Rx = R^* Q^* y$, d'où $Rx = Q^* y$.
La méthode des moindres carrés par factorisation QR se fait de la façon suivante:
# étape 1.
W, R = householder(A)
R_hat = R[0:2][:,0:2]
# étape 2.
qy = Qadjb(W,y)
qy_hat = qy[0:2]
# étape 3.
c = backsubs(R_hat, qy_hat)
La solution $c$ devrait être très similaire aux valeurs de $(a,b)$:
print c
print a, b
[[-0.19919044] [ 0.49872184]] -0.2 0.5
pl.plot(x,y,'b.');
xx = linspace(min(x),max(x),50)
pl.plot(xx,c[0]*xx + c[1],'k')
pl.xlabel('x');
pl.ylabel('y');