515 KiB
515 KiB
In [203]:
import numpy as np
import time
import scipy
from scipy.linalg import circulant
import scipy.io
import matplotlib.pyplot as plt
On veut résoudre: $$\min_x \frac{1}{2} ||Ax-b||^2_2$$ Le problème peut s'écrire comme: $$\min_x (a_{11}x_1+a_{12}x_2-b_1)^2 + ...$$ $$\frac{\delta f}{\delta x_2}=0 \Leftrightarrow 2a_{12}(a_{11}x_1+a_{12}x_2-b_1) + ...$$ $$(a_{12}^2+a_{22}^2+a_{32}^2)x_2 = -(A(:,2)^TA(:,1))x_1 + ...$$ $$(A(:,2)^TA(:,2))x_2 = -(A^TAx)_2 + (A^Tb)_2 + A(:,2)^TA(:,2)x_2$$ $$x_2=x_2-\frac{(A^TAx-A^Tb)_2}{(A^TA)_{2,2}}$$ $$x_i^+ = x_i-\frac{[\Delta f(x)]_i}{[A^TA]_{i,i}}$$
In [204]:
def CDLS(A,b,x0,maxiter=10**4,maxtime=2):
t0 = time.time()
t = np.zeros(maxiter)
fobj = np.zeros(maxiter)
x = x0.copy()
# O(mr^2)
AtA = A.T@A
# O(mr)
Atb = A.T@b
# O(m)
btb = b@b / 2
g = AtA@x-Atb
fobj[0] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2
k = 0
while k < maxiter-1 and t[k] < maxtime:
# On mets à jour le vecteur x entrée par entrée
for i in range(r):
xinew = x[i] - g[i] / AtA[i,i]
g = g - AtA[:,i]*x[i]+AtA[:,i]*xinew # O(r)
x[i] = xinew
fobj[k+1] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2
t[k+1] = time.time() - t0
k += 1
return t[:k+2], fobj[:k+2], x
In [205]:
def GRLS(A,b,x0,maxiter=10**4,maxtime=2):
t0 = time.time()
t = np.zeros(maxiter)
fobj = np.zeros(maxiter)
x = x0.copy()
# O(mr^2)
AtA = A.T@A
# O(mr)
Atb = A.T@b
# O(m)
btb = b@b / 2
g = AtA@x-Atb
L = max(np.linalg.eigvals(AtA))
alpha = 1/L
fobj[0] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2
k = 0
while k < maxiter-1 and t[k] < maxtime:
d = -g
x = x + alpha*d
g = AtA@x - Atb
fobj[k+1] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2
t[k+1] = time.time() - t0
k += 1
return t[:k+1], fobj[:k+1], x
In [206]:
def GALS(A,b,x0,maxiter=10**4,maxtime=2):
t0 = time.time()
t = np.zeros(maxiter)
fobj = np.zeros(maxiter)
x = x0.copy()
# O(mr^2)
AtA = A.T@A
# O(mr)
Atb = A.T@b
# O(m)
btb = b@b / 2
g = AtA@x-Atb
L = max(np.linalg.eigvals(AtA))
alpha = 1/L
fobj[0] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2
k = 0
xprec = x.copy()
while k < maxiter-1 and t[k] < maxtime:
y = x + ((k-1)/(k+2))*(x-xprec)
xprec = x
g = AtA@y - Atb
d = -g
x = y + alpha*d
fobj[k+1] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2
t[k+1] = time.time() - t0
k += 1
return t[:k+1], fobj[:k+1], x
In [207]:
#Pour la Question 3
m, r = 1000, 100
A = np.random.rand(m,r)
b = np.random.rand(m)
x0 = np.ones(r)
maxiter, maxtime = 10**4, 2
t_cd, f_cd, x_cd = CDLS(A,b,x0,maxiter=maxiter,maxtime=maxtime)
t_gr, f_gr, x_gr = GRLS(A,b,x0,maxiter=maxiter,maxtime=maxtime)
t_ga, f_ga, x_ga = GALS(A,b,x0,maxiter=maxiter,maxtime=maxtime)
In [208]:
#Pour la Question 3
xopt = np.linalg.solve(A.T@A,A.T@b)
fopt = 0.5*xopt@((A.T@A)@xopt)-((A.T@b)@xopt)+b@b/2
#Affichage première figure
plt.figure()
plt.semilogy(f_cd-fopt, label="Coord. Descent")
plt.semilogy(f_gr-fopt, label="Meth.Gradient")
plt.semilogy(f_ga-fopt, label="Meth.Gradient ACC")
plt.legend()
plt.plot()
#Affichage deuxième figure
plt.figure()
plt.semilogy(t_cd, f_cd-fopt, label="Coord. Descent")
plt.semilogy(t_gr, f_gr-fopt, label="Meth.Gradient")
plt.semilogy(t_ga, f_ga-fopt, label="Meth.Gradient ACC")
plt.legend()
plt.plot()
Out[208]:
In [209]:
def CDNNLS(A,b,x0,maxiter=10**4,maxtime=2):
t0 = time.time()
t = np.zeros(maxiter)
fobj = np.zeros(maxiter)
x = x0.copy()
# O(mr^2)
AtA = A.T@A
# O(mr)
Atb = A.T@b
# O(m)
btb = b@b / 2
g = AtA@x-Atb
fobj[0] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2
k = 0
while k < maxiter-1 and t[k] < maxtime:
# On mets à jour le vecteur x entrée par entrée
for i in range(r):
xinew = max(0, x[i] - g[i] / AtA[i,i])
g = g - AtA[:,i]*x[i]+AtA[:,i]*xinew # O(r)
x[i] = xinew
fobj[k+1] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2
t[k+1] = time.time() - t0
k += 1
return t[:k+2], fobj[:k+2], x
In [210]:
def GRNNLS(A,b,x0,maxiter=10**4,maxtime=2):
t0 = time.time()
t = np.zeros(maxiter)
fobj = np.zeros(maxiter)
x = x0.copy()
# O(mr^2)
AtA = A.T@A
# O(mr)
Atb = A.T@b
# O(m)
btb = b@b / 2
g = AtA@x-Atb
L = max(np.linalg.eigvals(AtA))
alpha = 1/L
fobj[0] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2
k = 0
while k < maxiter-1 and t[k] < maxtime:
d = -g
x = np.maximum(0, x + alpha*d)
g = AtA@x - Atb
fobj[k+1] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2
t[k+1] = time.time() - t0
k += 1
return t[:k+1], fobj[:k+1], x
In [211]:
def GANNLS(A,b,x0,maxiter=10**4,maxtime=2):
t0 = time.time()
t = np.zeros(maxiter)
fobj = np.zeros(maxiter)
x = x0.copy()
# O(mr^2)
AtA = A.T@A
# O(mr)
Atb = A.T@b
# O(m)
btb = b@b / 2
g = AtA@x-Atb
L = max(np.linalg.eigvals(AtA))
alpha = 1/L
fobj[0] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2
k = 0
xprec = x.copy()
while k < maxiter-1 and t[k] < maxtime:
y = x + ((k-1)/(k+2))*(x-xprec)
xprec = x
g = AtA@y - Atb
d = -g
x = np.maximum(0, y + alpha*d)
fobj[k+1] = 0.5*x@(AtA@x) - (Atb@x) + btb # 0.5*|A*x-b|_2^2
t[k+1] = time.time() - t0
k += 1
return t[:k+1], fobj[:k+1], x
In [212]:
#Pour la Question 4
m, r = 10000, 1000
A = np.random.rand(m,r)
b = np.random.rand(m)
x0 = np.ones(r)
maxiter, maxtime = 10**4, 2
t_cd, f_cd, x_cd = CDNNLS(A,b,x0,maxiter=maxiter,maxtime=maxtime)
t_gr, f_gr, x_gr = GRNNLS(A,b,x0,maxiter=maxiter,maxtime=maxtime)
t_ga, f_ga, x_ga = GANNLS(A,b,x0,maxiter=maxiter,maxtime=maxtime)
In [213]:
#Pour la Question 4
xopt, fopt = scipy.optimize.nnls(A,b)
fopt = (fopt**2)/2
#Affichage première figure
plt.figure()
plt.semilogy(f_cd-fopt, label="Coord. Descent")
plt.semilogy(f_gr-fopt, label="Meth.Gradient")
plt.semilogy(f_ga-fopt, label="Meth.Gradient ACC")
plt.legend()
plt.plot()
#Affichage deuxième figure
plt.figure()
plt.semilogy(t_cd, f_cd-fopt, label="Coord. Descent")
plt.semilogy(t_gr, f_gr-fopt, label="Meth.Gradient")
plt.semilogy(t_ga, f_ga-fopt, label="Meth.Gradient ACC")
plt.legend()
plt.plot()
Out[213]:
$$\min_{H\geq0, W\geq0} ||X-WH||_F^2$$
Si on considère $H$ fixe -> convexe
Si on considère $W$ fixe-> convexe
Pour l'algo on alterne entre les deux:
for i = maxiter
W <- Update(X,H)
H <- Update(X,H)
$||X-WH||_F^2=||X^T-H^TW^T||_F^2$
$H^* \leftarrow UpdateH(X,W)$
$W^* \leftarrow UpdateH(X^T,W^T)^T$
$||X-WH||_F^2 = ||x(:,1)-WH(:,1)||_2^2 + ||x(:,2)-WH(:,2)||_2^2 + \ldots$
Prenons un sous-problème:
$||x(:,2)-WH(:,2)||_2^2$
$\min_{x\geq0} ||Ax-b||_2^2$
=> NNLS
In [224]:
def initialisation(m,n,r):
W = np.random.random((m,r))
H = np.random.random((r,n))
return W, H
In [225]:
def NNLS(AtA,Atb,x,maxiter=10**2):
r = len(x)
g = AtA@x - Atb
for iter in range(maxiter):
for i in range(r):
xinew = np.maximum(0, x[i] - g[i] / AtA[i,i])
g = g - AtA[:,i]*x[i]+AtA[:,i]*xinew # O(r)
x[i] = xinew
return x
In [230]:
def optimizeH(X,W,H):
# Pour X et W donnés,
# ce code renvoie la matrice H minimisant ||X-WH||_F^T avec H>=0
# en résolvant un NNLS pour chaque colonne de H
n = X.shape[1]
AtA = W.T@W # O(mr^2)
AtB = W.T@X # O(mrn) <- le + cher
for j in range(n):
H[:,j] = NNLS(AtA, AtB[:,j], H[:,j])
In [236]:
def alternatealgo_l2nmf(X,r,maxiter):
m, n = X.shape
W, H = initialisation(m,n,r)
fobj = np.zeros(maxiter)
for iter in range(maxiter):
optimizeH(X,W,H) # optimisation de H étant donnés X et W
optimizeH(X.T, H.T, W.T) # optimisation de W étant donnés X^T et W^T
return W, H
In [237]:
#Pour la Question 5
X = np.array([[0 , 1 , 2 , 2 , 1 , 0],
[0 , 0 , 1 , 2 , 2 , 1],
[1 , 0 , 0 , 1 , 2 , 2],
[2 , 1 , 0 , 0 , 1 , 2],
[2 , 2 , 1 , 0 , 0 , 1],
[1 , 2 , 2 , 1 , 0 , 0]])
In [238]:
#Pour la Question 5
def load_matrix(name):
S = scipy.io.loadmat(name)
X = S["M"]
return X
def monaffichage(X,nbimages,nbrow,nbcol,m_image,n_image):
if nbrow*nbcol==nbimages and m_image*n_image==X.shape[0]:
Xtot = np.zeros((nbrow*m_image,nbcol*n_image))
i, j = 0, 0
for k in range(nbimages):
reshaped_image = np.reshape(X[:,k], (m_image,n_image)) / max(X[:,k])
Xtot[(i)*m_image:(i+1)*m_image,(j)*n_image:(j+1)*n_image] = reshaped_image.T
if j>=nbcol-1:
i+=1
j = 0
else: j+=1
plt.figure()
plt.imshow(1 - Xtot, cmap='gray')
nbvisages = 200;
X = load_matrix('exemplenmf2.mat')
Xtemp = X[:,:nbvisages]
monaffichage(Xtemp,200,10,20,19,19)
#Lancement de l'algorithme sur la matrice Xtemp avec r=5
r = 10
W, H = alternatealgo_l2nmf(Xtemp,r,maxiter=5)
#plt.figure()
#plt.semilogy(f)
#plt.xlabel("Itérations")
#plt.ylabel("Erreur ||X-WH||")
#plt.title("Evolution de l'erreur")
#plt.show()
In [239]:
#Affichage de W
monaffichage(W,10,2,5,19,19)
In [240]:
#Affichage de W
monaffichage(W@H,200,10,20,19,19)
