From d79a92accb05776b5b3a4b586aac921d62ff13f3 Mon Sep 17 00:00:00 2001 From: SamuelDiai Date: Thu, 1 Apr 2021 11:56:48 +0200 Subject: [PATCH] First push --- main.py | 140 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ utils.py | 46 ++++++++++++++++++ 2 files changed, 186 insertions(+) create mode 100644 main.py create mode 100644 utils.py diff --git a/main.py b/main.py new file mode 100644 index 0000000..d67f2e7 --- /dev/null +++ b/main.py @@ -0,0 +1,140 @@ +import math +import time +import cvxpy as cp +from utils import Laplacian_dual, Laplacian_inv, Laplacian + +class LearnGraphTopolgy: + def __init__(self, S, alpha=0, beta=1e4, n_iter=10000, c1=0., c2=1e10, tol = 1e-6): + self.tol = tol + self.S = S + self.p = S.shape[0] + self.n = np.int(self.p * (self.p - 1)/2) + self.alpha = alpha + self.beta = beta + self.n_iter = n_iter + self.c1 = c1 + self.c2 = c2 + + def w_init(self, w0_init, Sinv): + """ + Initialization of w + """ + if w0_init == 'naive': + w0 = Laplacian_inv(Sinv, self.p, self.n) + w0[w0<0] = 0 + else: + raise ValueError('Method not implemented') + return w0 + + def update_w(self, w, Lw, U, lambda_, K): + """ + Compute w update as equation 12 + """ + c = Laplacian_dual(U @ np.diag(lambda_) @ U.T - K / self.beta, self.p, self.n) + grad_f = Laplacian_dual(Lw, self.p, self.n) - c + M_grad_f = Laplacian_dual(Laplacian(grad_f, self.p), self.p, self.n) + wT_M_grad_f = sum(w * M_grad_f) + dwT_M_dw = sum(grad_f * M_grad_f) + ## exact line search + t = (wT_M_grad_f - sum(c * grad_f)) / dwT_M_dw + ## no line search : + # p = int(0.5*(1 + np.sqrt(1 + 8*w.shape[0]))) + # t = 0.5*p + w_new = w - t * grad_f + w_new[w_new < 0] = 0 + return w_new + + def adjacency(self, w): + ''' + Compute the Adjacency matrix from w + ''' + Aw = np.zeros((self.p,self.p)) + k=0 + for i in range(0, self.p): + for j in range(i+1,self.p): + Aw[i][j] = w[k] + k = k + 1 + Aw = Aw + Aw.T + return Aw + + def update_lambda(self, U, Lw, k): + """ + Compute lambda update as proposed in the supplementary + """ + q = Lw.shape[1] - k + d = np.diag(np.dot(U.T, np.dot(Lw, U))) + assert(d.shape[0] == q) + lambda_ = (d + np.sqrt(d**2 + 4/self.beta))/2 + + cond = (lambda_[q-1] - self.c2 <= 1e-9) and (lambda_[0] - self.c1 >= -1e-9) and np.all(lambda_[1:q] - lambda_[0:q-1] >= -1e-9) + + if cond: + return lambda_ + else: + lambda_[lambda_ < self.c1] = self.c1 + lambda_[lambda_ > self.c2] = self.c2 + cond = (lambda_[q-1] - self.c2 <= 1e-9) and (lambda_[0] - self.c1 >= -1e-9) and np.all(lambda_[1:q] - lambda_[0:q-1] >= -1e-9) + + if cond: + return lambda_ + else: + raise ValueError("Consider increasing value of beta") + + def update_U(self, Lw, k): + """ + Compute U update as equation 14 + """ + _, eigvec = np.linalg.eigh(Lw) + assert(eigvec.shape[1] == self.p) + return eigvec[:, k:] + + def objective(self, Lw, lambda_, K, U): + """ + Compute objective function - equation 8 + """ + term1 = np.sum(-np.log(lambda_)) + term2 = np.trace(np.dot(K, Lw)) + term3 = 0.5 * self.beta * np.linalg.norm(Lw - np.dot(U, np.dot(np.diag(lambda_), U.T)), ord='fro')**2 + return term1 + term2 + term3 + + def learn_graph(self, k=1, w0_init='naive', eps = 1e-4): + + # find an appropriate inital guess + Sinv = np.linalg.pinv(self.S) + # if w0 is either "naive" or "qp", compute it, else return w0 + w0 = self.w_init(w0_init, Sinv) + # compute quantities on the initial guess + Lw0 = Laplacian(w0, self.p) + # l1-norm penalty factor + H = self.alpha * (np.eye(self.p) - np.ones((self.p, self.p))) + K = self.S + H + U0 = self.update_U(Lw = Lw0, k = k) + lambda_0 = self.update_lambda(U = U0, Lw = Lw0, k = k) + + objective_seq = [] + + for _ in range(self.n_iter): + w = self.update_w(w = w0, Lw = Lw0, U = U0, lambda_ = lambda_0, K = K) + Lw = Laplacian(w, self.p) + U = self.update_U(Lw = Lw, k = k) + lambda_ = self.update_lambda(U = U, Lw = Lw, k = k) + + # check for convergence + convergence = np.linalg.norm(w0 - w, ord=2) < self.tol + objective_seq.append(self.objective(Lw, lambda_, K, U)) + + if convergence: + break + + # update estimates + w0 = w + U0 = U + lambda_0 = lambda_ + Lw0 = Lw + K = self.S + H / (-Lw + eps) + + # compute the adjancency matrix + Aw = self.adjacency(w) + results = {'laplacian' : Lw, 'adjacency' : Aw, 'w' : w, 'lambda' : lambda_, 'U' : U, + 'convergence' : convergence, 'objective_seq' : objective_seq} + return results diff --git a/utils.py b/utils.py new file mode 100644 index 0000000..f019947 --- /dev/null +++ b/utils.py @@ -0,0 +1,46 @@ +import numpy as np + +def Laplacian_inv(M, p, n): + ''' + Computes the inverse operation of L + M : shape p*p + w : shape p*(p - 1)/2 + ''' + ind = 0 + w = np.zeros(n) + for i in range(0, p): + for j in range(i+1, p): + w[ind] = -M[i][j] + ind += 1 + return w + +def Laplacian(w, p): + """ + Compute laplacian from a vector + + w : shape p*(p-1)/2 + LW : shape p*p + """ + Lw = np.zeros((p, p)) + ind = 0 + for i in range(p-1): + for j in range(i+1, p): + Lw[i, j] = -w[ind] + ind += 1 + Lw = Lw + Lw.T + np.fill_diagonal(Lw, -(Lw.sum(axis = 0) - np.trace(Lw))) + return Lw + +def Laplacian_dual(M, p, n): + """ + Compute the dual of L operator + + M : shape p*p + (L*)M : shape p*(p-1)/2 + """ + L_Y = np.zeros(n) + for j in range(p): + for i in range(j+1, p): + k = np.int(i - j + j*(2*p - j - 1)/2) - 1 + L_Y[k] = M[i, i] - M[i, j] - M[j, i] + M[j, j] + return L_Y