From afb64297845db5c0492bf921456f2c1923907b43 Mon Sep 17 00:00:00 2001
From: Navjot Singh <navjot@vpnpool-10-250-19-203.near.illinois.edu>
Date: Thu, 1 Oct 2020 20:34:45 -0500
Subject: [PATCH]  adding explicit als, gn and their test

---
 CP_GN.py             | 229 ++++++++++++++++++++++++++++++++++++
 arg_defs.py          |   6 +
 backend/ctf_ext.py   | 189 ++++++++++++++++++++++++++++++
 backend/numpy_ext.py | 273 +++++++++++++++++++++++++++++++++++++++++++
 explicit_als.py      |  74 ++++++++++++
 script_GN.py         | 213 +++++++++++++++++++++++++++++++++
 6 files changed, 984 insertions(+)
 create mode 100644 CP_GN.py
 create mode 100644 backend/ctf_ext.py
 create mode 100644 backend/numpy_ext.py
 create mode 100644 explicit_als.py
 create mode 100644 script_GN.py

diff --git a/CP_GN.py b/CP_GN.py
new file mode 100644
index 0000000..6168b92
--- /dev/null
+++ b/CP_GN.py
@@ -0,0 +1,229 @@
+ #!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Sat Apr 18 14:30:39 2020
+
+@author: navjot
+"""
+import numpy as np
+import numpy.linalg as la
+import time
+
+
+class CP_GN_Completer():
+    
+    def __init__(self,tenpy, T, Omega, A ):
+        self.tenpy = tenpy
+        self.T = T
+        self.Omega = Omega
+        self.A = A
+        self.atol = 0
+        self.cg_tol = 1e-03
+        self.maxiter = 500
+        self.atol = 0
+        self.total_cg = 0
+
+    def fast_block_diag_precondition(self,x,regu):
+        V = []
+        for i in range(len(self.A)):
+            lst_mat = []
+            for j in range(len(self.A)):
+                if i != j :
+                    lst_mat.append(self.A[j])
+                else:
+                    lst_mat.append(self.tenpy.zeros(self.A[i].shape)) 
+            if self.tenpy.name() == "numpy":
+                V.append(self.tenpy.Solve_Factor(self.Omega,lst_mat,x[i],i,regu))
+            else:
+                self.tenpy.Solve_Factor(self.Omega,lst_mat,x[i],i)
+                V.append(lst_mat[i])
+        return V
+
+        
+    def matvec(self,regu,delta):
+        N = len(self.A)
+        ret = []
+        lst_mat = self.A[:]
+        lst_mat[0] = delta[0].copy()
+        inter = self.tenpy.TTTP(self.Omega, lst_mat)
+        for n in range(1,N):
+            ret.append(self.tenpy.zeros(self.A[n].shape))
+            lst_mat= self.A[:]
+            lst_mat[n] = delta[n].copy()
+            inter += self.tenpy.TTTP(self.Omega, lst_mat)
+
+        lst_mat = self.A[:]
+        lst_mat[0] = self.tenpy.zeros(self.A[0].shape)
+        self.tenpy.MTTKRP(inter,lst_mat,0)
+        ret.append(self.tenpy.zeros(self.A[0].shape))
+        ret[0]+=lst_mat[0]
+        ret[0]+= regu*delta[0]
+        for n in range(1,N):
+            lst_mat = self.A[:]
+            lst_mat[n] = self.tenpy.zeros(self.A[n].shape)
+            self.tenpy.MTTKRP(inter,lst_mat,n)
+            ret[n]+=lst_mat[n]
+            ret[n]+= regu*delta[n]
+
+        return ret
+
+
+    def fast_precond_conjugate_gradient(self,g,Regu):
+        start = time.time()
+        
+        x = [self.tenpy.zeros(A.shape) for A in g]
+        
+        g_norm = self.tenpy.list_vecnorm(g)
+            
+
+        tol = np.max([self.atol,np.min([self.cg_tol,np.sqrt(g_norm)])])*g_norm
+        
+        if g_norm<tol:
+            return x
+
+        z = self.fast_block_diag_precondition(g,Regu)
+
+        p = z
+
+        counter = 0
+        while True:
+            mv = self.matvec(Regu,p)
+
+            mul = self.tenpy.mult_lists(g,z)
+
+            alpha = mul/self.tenpy.mult_lists(p,mv) 
+
+            x =self.tenpy.scl_list_add(alpha,x,p)
+
+            g = self.tenpy.scl_list_add(-1*alpha,g,mv)
+            
+            
+            if self.tenpy.list_vecnorm(g)<tol:
+                counter+=1
+                #end = time.time()
+                break
+
+            z = self.fast_block_diag_precondition(g,Regu)
+
+            beta = self.tenpy.mult_lists(g,z)/mul
+
+            p = self.tenpy.scl_list_add(beta,z,p)
+
+            counter += 1
+            
+            if counter == self.maxiter:
+                #end = time.time()
+                break
+                
+        end = time.time()
+        self.tenpy.printf("cg took:",end-start)
+        self.tenpy.printf("CG iterations is",counter)
+
+        return x,counter
+
+
+    def fast_conjugate_gradient(self,g,Regu):
+        start = time.time()
+
+        x = [self.tenpy.zeros(A.shape) for A in g]
+        
+        g_norm = self.tenpy.list_vecnorm(g)
+
+        tol = np.max([self.atol,np.min([self.cg_tol,np.sqrt(g_norm)])])*g_norm
+        
+        
+        r = g
+        
+        #self.tenpy.printf('starting res in cg is',self.tenpy.list_vecnorm(r))
+        if g_norm <tol:
+            return x
+        
+        p = r
+        counter = 0
+
+        while True:
+            mv = self.matvec(Regu,p)
+
+            prod = self.tenpy.mult_lists(p,mv)
+
+            alpha = self.tenpy.mult_lists(r,r)/prod
+
+            x = self.tenpy.scl_list_add(alpha,x,p)
+
+            r_new = self.tenpy.scl_list_add(-1*alpha,r,mv)
+                
+            #self.tenpy.printf('res in cg is',self.tenpy.list_vecnorm(r_new))
+
+            if self.tenpy.list_vecnorm(r_new)<tol:
+                counter+=1
+                end = time.time()
+                break
+            beta = self.tenpy.mult_lists(r_new,r_new)/self.tenpy.mult_lists(r,r)
+
+            p = self.tenpy.scl_list_add(beta,r_new,p)
+            r = r_new
+            counter += 1
+
+            if counter == self.maxiter:
+                end = time.time()
+                break
+                
+        self.tenpy.printf('cg took',end-start)
+        self.tenpy.printf('Number of cg iterations is :',counter)
+        
+
+        return x,counter
+
+    def Get_RHS(self):
+        grad = []
+        inter = self.tenpy.TTTP(self.Omega, self.A)
+        inter = self.T - inter
+        for i in range(len(self.A)):
+            lst_mat = self.A[:]
+            lst_mat[i] = self.tenpy.zeros(self.A[i].shape)
+            self.tenpy.MTTKRP(inter,lst_mat,i)
+            grad.append(lst_mat[i])
+        return grad
+    
+    def update_A(self,delta):
+        for i in range(len(delta)):
+            self.A[i] += delta[i]
+            
+    def step(self,Regu):
+        g = self.Get_RHS()
+        #P = self.Compute_preconditioner(Regu)
+        [delta,counter] = self.fast_precond_conjugate_gradient(g,Regu)
+        #[delta,counter] = self.fast_conjugate_gradient(g,Regu)
+        self.total_cg+= counter
+        print("TOTAL CG ITERATIONS :",self.total_cg)
+        self.update_A(delta)
+        
+        return self.A
+
+def getCPGN(tenpy, T_in, T, O, U, V, W, reg_GN,I,J,K,R, num_iter_GN):
+    opt = CP_GN_Completer(tenpy, T_in, O, [U,V,W])
+
+    #if T_in.sp == True:
+    #    nnz_tot = T_in.nnz_tot
+    #else:
+    #    nnz_tot = ctf.sum(omega)
+    nnz_tot = tenpy.sum(O)
+        
+    regu = reg_GN
+    print("--------------------------------GN WIth  CG-----------------------------")
+    start= time.time()
+    # T_in = backend.einsum('ijk,ijk->ijk',T,O)
+    it = 0
+
+    for i in range(num_iter_GN):
+        it += 1
+        [U,V,W] = opt.step(regu)
+        print("After " + str(it) + " iterations,")
+        print("RMSE",(tenpy.vecnorm(tenpy.TTTP(O,[U,V,W])-T_in))/(tenpy.sum(O))**0.5)
+        print("Full Tensor Objective",(tenpy.norm(tenpy.einsum('ir,jr,kr->ijk',U,V,W)-T)))
+
+    end= time.time()
+
+    print('GN time taken is ',end - start)
+    
+    return [U,V,W]
\ No newline at end of file
diff --git a/arg_defs.py b/arg_defs.py
index f695c45..d37691d 100644
--- a/arg_defs.py
+++ b/arg_defs.py
@@ -24,6 +24,12 @@ def add_general_arguments(parser):
         default=10,
         metavar='int',
         help='Input CP decomposition rank (default: 10)')
+    parser.add_argument(
+        '--num-iter-GN',
+        type=int,
+        default=10,
+        metavar='int',
+        help='Number of iterations (sweeps) to run GN (default: 10)')
     parser.add_argument(
         '--num-iter-ALS-implicit',
         type=int,
diff --git a/backend/ctf_ext.py b/backend/ctf_ext.py
new file mode 100644
index 0000000..40833c5
--- /dev/null
+++ b/backend/ctf_ext.py
@@ -0,0 +1,189 @@
+import ctf
+import numpy as np
+
+def _get_num_str(n):
+    allstr = "abcdefghijklmonpqrstuvwzyx0123456789,./;'][=-`"
+    return allstr[0:n]
+
+
+def name():
+    return 'ctf'
+
+def diag(v):
+    return ctf.diag(v)
+
+def save_tensor_to_file(T, filename):
+    np.save(filename, T.to_nparray())
+
+def load_tensor_from_file(filename):
+    try:
+        T = np.load(filename)
+        print('Loaded tensor from file ', filename)
+    except FileNotFoundError:
+        raise FileNotFoundError('No tensor exist on: ', filename)
+    return ctf.from_nparray(T)
+
+def from_nparray(arr):
+    return ctf.from_nparray(arr)
+
+def TTTP(T, A):
+    return ctf.TTTP(T, A)
+
+def MTTKRP(T,A,idx):
+    return ctf.MTTKRP(T, A, idx)
+
+def is_master_proc():
+    if ctf.comm().rank() == 0:
+        return True
+    else:
+        return False
+
+def printf(*string):
+    if ctf.comm().rank() == 0:
+        print(string)
+
+def tensor(shape, sp, *args):
+    return ctf.tensor(shape, sp, *args)
+
+def sparse_random(shape, begin, end, sp_frac):
+    tensor = ctf.tensor(shape, sp=True)
+    tensor.fill_sp_random(begin, end, sp_frac)
+    return tensor
+
+def vecnorm(T):
+    return ctf.vecnorm(T)
+
+def Solve_Factor(T,mats,rhs,num):
+    return ctf.Solve_Factor(T,mats,rhs,num)
+    
+def list_vecnormsq(list_A):
+    l = [i**2 for i in list_A]
+    s = 0
+    for i in range(len(l)):
+        s+=ctf.sum(l[i])
+    return s
+    
+def list_vecnorm(list_A):
+    l = [i**2 for i in list_A]
+    s = 0
+    for i in range(len(l)):
+        s+=ctf.sum(l[i])
+        
+    return s**0.5
+    
+    
+def list_add(list_A, list_B):
+    return [A + B for (A, B) in zip(list_A, list_B)]
+
+
+def scalar_mul(sclr, list_A):
+    return [sclr * A for A in list_A]
+
+
+def mult_lists(list_A,list_B):
+    s = 0
+    for i in range(len(list_A)):
+        s+= ctf.einsum('ij,ij->',list_A[i],list_B[i])
+    return s
+
+def scl_list_add(scl,list_A,list_B):
+    x= []
+    for i in range(len(list_A)):
+        x.append(list_A[i]+scl*list_B[i])
+            
+    return x
+    
+def norm(v):
+    return v.norm2() 
+
+def dot(A,B):
+    return ctf.dot(A,B)
+
+def svd(A,r=None):
+    return ctf.svd(A,r)
+
+def svd_rand(A,r):
+    return ctf.svd_rand(A,r)
+
+def cholesky(A):
+    return ctf.cholesky(A)
+
+def solve_tri(A, B, lower=True, from_left=False, transp_L=False):
+    return ctf.solve_tri(A, B, lower, from_left, transp_L)
+
+def einsum(string, *args):
+    if "..." in string:
+        left = string.split(",")
+        left[-1], right = left[-1].split("->")
+        symbols = "".join(list(set([chr(i) for i in range(48, 127)]) - set(string.replace(".", "").replace(",", "").replace("->", ""))))
+        symbol_idx = 0
+        for i, (s, tsr) in enumerate(zip(left, args)):
+            num_missing = tsr.ndim - len(s.replace("...", ""))
+            left[i] = s.replace("...", symbols[symbol_idx:symbol_idx+num_missing])
+            symbol_idx += num_missing
+        right = right.replace("...", symbols[:symbol_idx])
+        string = ",".join(left) + "->" + right
+
+    return ctf.einsum(string, *args)
+
+def ones(shape):
+    return ctf.ones(shape)
+
+def zeros(shape):
+    return ctf.zeros(shape)
+
+def dot_product(a,b):
+    mystr = _get_num_str(a.ndim)
+    einstr = mystr +"," + mystr + "->"
+    return einsum(einstr,a,b)
+
+def sum(A, axes=None):
+    return ctf.sum(A, axes)
+
+def random(shape):
+    return ctf.random.random(shape)
+
+def seed(seed):
+    return ctf.random.seed(seed)
+    
+def list_add(list_A,list_B):
+    return [A+B for (A,B) in zip(list_A,list_B)]
+    
+def scalar_mul(sclr,list_A):
+    return [sclr*A for A in list_A]
+
+def speye(*args):
+    return ctf.speye(*args)
+
+def eye(*args):
+    return ctf.eye(*args)
+
+def transpose(A):
+    return ctf.transpose(A)
+
+def argmax(A, axis=0):
+    return abs(A).to_nparray().argmax(axis=axis)
+
+def asarray(T):
+    return ctf.astensor(T)
+    
+
+def reshape(A,shape,order='F'):
+    return ctf.reshape(A,shape,order)
+
+def einsvd(einstr, A, r=None, transpose=True, compute_uv=True, full_matrices=True, mult_sv=False):
+    str_a, str_uv = einstr.split("->")
+    str_u, str_v = str_uv.split(",")
+    U, S, Vh = A.i(str_a).svd(str_u, str_v, rank=r)
+    if not compute_uv:
+        return S
+    if mult_sv:
+        char_i = list(set(str_v) - set(str_a))[0]
+        char_s = list(set(str_a) - set(str_v))[0]
+        Vh = ctf.einsum(char_s + char_i + "," + str_v + "->" + str_v.replace(char_i, char_s), ctf.diag(S), Vh)
+    if not transpose:
+        Vh = Vh.T
+    return U, S, Vh
+
+def squeeze(A):
+    return A.reshape([s for s in A.shape if s != 1])
diff --git a/backend/numpy_ext.py b/backend/numpy_ext.py
new file mode 100644
index 0000000..0dc55c4
--- /dev/null
+++ b/backend/numpy_ext.py
@@ -0,0 +1,273 @@
+import numpy as np
+import numpy.linalg as la
+import scipy.linalg as sla
+
+def name():
+    return 'numpy'
+
+def fill_diagonal(matrix, value):
+    return np.fill_diagonal(matrix, value)
+
+def diag(v):
+    return np.diag(v)
+
+def save_tensor_to_file(T, filename):
+    np.save(filename, T)
+
+def load_tensor_from_file(filename):
+    try:
+        T = np.load(filename)
+        print('Loaded tensor from file ', filename)
+    except FileNotFoundError:
+        raise FileNotFoundError('No tensor exist on: ', filename)
+    return T
+
+def TTTP(T, A):
+    T_inds = "".join([chr(ord('a')+i) for i in range(T.ndim)])
+    einstr = ""
+    A2 = []
+    for i in range(len(A)):
+        if A[i] is not None:
+            einstr += chr(ord('a')+i) + chr(ord('a')+T.ndim) + ','
+            A2.append(A[i])
+    einstr += T_inds + "->" + T_inds
+    A2.append(T)
+    return np.einsum(einstr, *A2)
+
+def MTTKRP(T,A,idx):
+    T_inds = "".join([chr(ord('a')+i) for i in range(T.ndim)])
+    einstr = ""
+    A2 = []
+    for i in range(len(A)):
+        if i != idx:
+            einstr += chr(ord('a')+i) + chr(ord('a')+T.ndim) + ','
+            A2.append(A[i])
+    einstr += T_inds + "->" + chr(ord('a')+idx) + chr(ord('a')+T.ndim)
+    A2.append(T)
+    A[idx][:] = np.einsum(einstr, *A2,optimize=True)
+    
+    
+def is_master_proc():
+    return True
+
+def Solve_Factor(Omega,A,RHS,num,reg):
+    N = len(A)
+    R = A[0].shape[1]
+    lst_mat = []
+    T_inds = "".join([chr(ord('a')+i) for i in range(Omega.ndim)])
+    einstr=""
+    for i in range(N):
+        if i != num:
+            einstr+=chr(ord('a')+i) + 'r' + ','
+            lst_mat.append(A[i])
+            einstr+=chr(ord('a')+i) + 'z' + ','
+            lst_mat.append(A[i])
+    einstr+= T_inds + "->"+chr(ord('a')+num)+'rz'
+    lst_mat.append(Omega)
+    P = np.einsum(einstr,*lst_mat,optimize=True)
+    o = np.zeros_like(RHS)
+    for j in range(A[num].shape[0]):
+        o[j,:] = la.solve(P[j]+reg*np.eye(R),RHS[j,:])
+    return o
+
+def printf(*string):
+    print(string)
+
+def tensor(shape,sp, *args2):
+    return np.ndarray(shape, *args2)
+
+
+def list_add(list_A,list_B):
+    return [A+B for (A,B) in zip(list_A,list_B)]
+
+def scalar_mul(sclr,list_A):
+    return [sclr*A for A in list_A]
+
+def mult_lists(list_A,list_B):
+    s = 0
+    for i in range(len(list_A)):
+        s+= np.einsum('ij,ij->',list_A[i],list_B[i],optimize=True)
+    return s
+
+def list_vecnormsq(list_A):
+    l = [i**2 for i in list_A]
+    s = 0
+    for i in range(len(l)):
+        s+=np.sum(l[i])
+    return s
+
+def list_vecnorm(list_A):
+    l = [i**2 for i in list_A]
+    s = 0
+    for i in range(len(l)):
+        s+=np.sum(l[i])
+        
+    return s**0.5
+
+def scl_list_add(scl,list_A,list_B):
+    x= []
+    for i in range(len(list_A)):
+        x.append(list_A[i]+scl*list_B[i])
+            
+    return x
+    
+
+def sparse_random(shape, begin, end, sp_frac):
+    tensor = np.random.random(shape)*(end-begin) + begin
+    mask = np.random.random(shape)<sp_frac
+    tensor = tensor * mask
+    return tensor
+
+def vecnorm(T):
+    return la.norm(np.ravel(T))
+
+def norm(v):
+    return la.norm(v)
+
+def dot(A,B):
+    return np.dot(A,B)
+
+def eigvalh(A):
+    return la.eigvalh(A)
+
+def eigvalsh(A):
+    return la.eigvalsh(A)
+
+def svd(A,r=None):
+    U,s,VT = la.svd(A,full_matrices=False)
+    if r is not None:
+        U = U[:,:r]
+        s = s[:r]
+        VT = VT[:r,:]
+    return U,s,VT
+
+def inv(A):
+    return la.inv(A)
+
+def svd_rand(A,r=None):
+    return svd(A,r)
+
+def cholesky(A):
+    return la.cholesky(A)
+
+def solve_tri(A, B, lower=True, from_left=True, transp_L=False):
+    if not from_left:
+        B = B.T
+        A = A.T
+        llower = not lower
+        X = sla.solve_triangular(A, B, trans=transp_L, lower=llower)
+        return X.T
+    else:
+        return sla.solve_triangular(A, B, trans=transp_L,lower=lower)
+
+def einsum(string, *args):
+    out = np.einsum(string, *args)
+    return out
+
+def ones(shape):
+    return np.ones(shape)
+
+def zeros(shape):
+    return np.zeros(shape)
+
+def dot_product(a,b):
+    return np.inner(a.reshape(-1),b.reshape(-1))
+
+def sum(A, axes=None):
+    return np.sum(A, axes)
+
+def random(shape):
+    return np.random.random(shape)
+
+def seed(seed):
+    return np.random.seed(seed)
+
+def asarray(T):
+    return np.array(T)
+
+def speye(*args):
+    return np.eye(*args)
+
+def eye(*args):
+    return np.eye(*args)
+
+def transpose(A):
+    return A.T
+
+def argmax(A, axis=0):
+    return abs(A).argmax(axis=axis)
+
+def qr(A):
+    return la.qr(A)
+
+def reshape(A,shape,order='F'):
+    return np.reshape(A,shape,order)
+
+def einsvd(operand, tns, r=None, transpose=True, compute_uv=True, full_matrices=True, mult_sv=False):
+    ''' compute SVD of tns
+
+    oprand (str): in form 'src -> tgta, tgtb'
+    tns (ndarray): tensor to be decomposed
+    transpose (bool): True iff VT(H) is required instead of V
+    compute_uv, full_matrices (bool): see numpy.linalg.svd
+
+    REQUIRE: only one contracted index
+    '''
+
+    src, _, tgt = operand.replace(' ', '').partition('->')
+    tgta, _, tgtb = tgt.partition(',')
+
+    # transpose and reshape to the matrix to be SVD'd
+    tgt_idx = set(tgta).union(set(tgtb))
+    contract_idx = str(list(tgt_idx.difference(set(src)))[0])
+    new_idx = (tgta + tgtb).replace(contract_idx, '')
+    trsped = np.einsum(src + '->' + new_idx, tns)
+
+    # do svd
+    shape = tns.shape
+    letter2size = {}
+    for i in range(len(src)):
+        letter2size[src[i]] = shape[i]
+    col_idx = tgtb.replace(contract_idx, '')
+    ncol = 1
+    for letter in col_idx:
+        ncol *= letter2size[letter]
+    mat = trsped.reshape((-1, ncol))
+    if not compute_uv:
+        return la.svd(mat, compute_uv=False)
+
+    # if u, v are needed
+    u, s, vh = la.svd(mat, full_matrices=full_matrices)
+
+    if r != None and r < len(s):
+        u = u[:,:r]
+        s = s[:r]
+        vh = vh[:r,:]
+    if mult_sv:
+        vh = np.dot(np.diag(s),vh)
+
+    # reshape u, v into shape (..., contract) and (contract, ...)
+    row_idx = tgta.replace(contract_idx, '')
+    shapeA = []
+    shapeB = [-1]
+    for letter in row_idx:
+        shapeA.append(letter2size[letter])
+    for letter in col_idx:
+        shapeB.append(letter2size[letter])
+    shapeA.append(-1)
+    u = u.reshape(shapeA)
+    vh = vh.reshape(shapeB)
+
+    # transpose u and vh into tgta and tgtb
+    preA = tgta.replace(contract_idx, '') + contract_idx
+    preB = contract_idx + tgtb.replace(contract_idx, '')
+    u = np.einsum(preA + '->' + tgta, u)
+    vh = np.einsum(preB + '->' + tgtb, vh)
+
+    # return
+    if not transpose:
+        vh = np.conj(vh.T)
+    return u, s, vh
+
+def squeeze(A):
+    return A.squeeze()
diff --git a/explicit_als.py b/explicit_als.py
new file mode 100644
index 0000000..db67c06
--- /dev/null
+++ b/explicit_als.py
@@ -0,0 +1,74 @@
+ #!/usr/bin/env python3
+
+import numpy as np
+import numpy.linalg as la
+import time
+
+
+class explicit_als_Completer():
+    
+    def __init__(self,tenpy, T, Omega, A ):
+        self.tenpy = tenpy
+        self.T = T
+        self.Omega = Omega
+        self.A = A
+
+        
+
+    def Get_RHS(self,num):
+        lst_mat = []
+        for j in range(len(self.A)):
+            if j != num :
+                lst_mat.append(self.A[j])
+            else:
+                lst_mat.append(self.tenpy.zeros(self.A[num].shape))
+
+        self.tenpy.MTTKRP(self.T,lst_mat,num)
+        grad = lst_mat[num]
+        
+        return grad
+            
+    def step(self,regu):
+        for i in range(len(self.A)):
+            lst_mat = []
+            for j in range(len(self.A)):
+                if i != j :
+                    lst_mat.append(self.A[j])
+                else:
+                    lst_mat.append(self.tenpy.zeros(self.A[i].shape))
+            g = self.Get_RHS(i)
+            if self.tenpy.name() == "numpy": 
+                self.A[i] = self.tenpy.Solve_Factor(self.Omega,lst_mat,g,i,regu)
+            else:
+                self.tenpy.Solve_Factor(self.Omega,lst_mat,g,i)
+                self.A[i] = lst_mat[i]
+    
+        return self.A
+
+def explicit_als(tenpy, T_in, T, O, U, V, W, reg_als,I,J,K,R, num_iter_als):
+    opt = explicit_als_Completer(tenpy, T_in, O, [U,V,W])
+
+    #if T_in.sp == True:
+    #    nnz_tot = T_in.nnz_tot
+    #else:
+    #    nnz_tot = ctf.sum(omega)
+    nnz_tot = tenpy.sum(O)
+        
+    regu = reg_als
+    print("--------------------------------explicit_als-----------------------------")
+    start= time.time()
+    # T_in = backend.einsum('ijk,ijk->ijk',T,O)
+    it = 0
+
+    for i in range(num_iter_als):
+        it += 1
+        [U,V,W] = opt.step(regu)
+        print("After " + str(it) + " iterations,")
+        print("RMSE",(tenpy.vecnorm(tenpy.TTTP(O,[U,V,W])-T_in))/(tenpy.sum(O))**0.5)
+        print("Full Tensor Objective",(tenpy.norm(tenpy.einsum('ir,jr,kr->ijk',U,V,W)-T)))
+
+    end= time.time()
+
+    print('Explicit als time taken is ',end - start)
+    
+    return [U,V,W]
\ No newline at end of file
diff --git a/script_GN.py b/script_GN.py
new file mode 100644
index 0000000..645a37f
--- /dev/null
+++ b/script_GN.py
@@ -0,0 +1,213 @@
+import ctf
+import time
+import random
+import sys
+import numpy as np
+import numpy.linalg as la
+from ctf import random as crandom
+
+import gzip
+import shutil
+import os
+import argparse
+import arg_defs as arg_defs
+
+import backend.ctf_ext as tenpy
+import backend.numpy_ext as tenpy_np
+
+from als import getALS_CG
+from sgd import sparse_SGD
+from CP_GN import getCPGN
+from explicit_als import explicit_als
+
+
+glob_comm = ctf.comm()
+#netflix_tensor_dims = (370169, 500, 2170)
+
+def getOmega(T):
+    [inds, data] = T.read_local_nnz()
+    data[:] = 1.
+    Omega = ctf.tensor(T.shape, sp=T.sp)
+    Omega.write(inds, data)
+    return Omega
+
+def create_lowr_tensor(I, J, K, r, sp_frac, use_sp_rep):
+    np.random.seed(1122)
+    ctf.random.seed(100)
+    #U = ctf.random.random((I, r))
+    U_np = np.random.randn(I,r)
+    U = ctf.astensor(U_np, dtype=np.float64)
+    #V = ctf.random.random((J, r))
+    V_np = np.random.randn(J,r)
+    V = ctf.astensor(V_np,dtype=np.float64)
+    
+    #W = ctf.random.random((K, r))
+    W_np = np.random.randn(K,r)
+    W = ctf.astensor(W_np, dtype=np.float64)
+    
+    T_in = ctf.tensor((I, J, K), sp=use_sp_rep)
+    T = ctf.ones((I, J, K))
+    T = tenpy.TTTP(T, [U, V, W])
+    
+    T_in.fill_sp_random(1, 1, sp_frac)
+    
+    T_in = ctf.TTTP(T_in,[U,V,W])
+    
+    #T_in+= ctf.einsum('ijk,ijk->ijk',noise,T_in)
+    
+    [inds, data] = T_in.read_local_nnz()
+    data[:] = 1.
+    Omega = ctf.tensor(T_in.shape,sp=T.sp)
+    Omega.write(inds,data)
+#    T = ctf.ones((I,J,K))
+    #T = ctf.tensor((I, J, K), sp=use_sp_rep)
+    #T.fill_sp_random(1, 1, sp_frac)
+#    T = ctf.TTTP(T, [U, V, W])
+    #np.random.seed(100)
+    #omega_np = np.ones(I*J*K)
+    #omega_np[:int(sp_frac*(I*J*K))] = 0
+    #np.random.shuffle(omega_np)
+   #omega = ctf.astensor(omega_np, dtype=np.float64)
+    #omega = omega_np.reshape((I, J, K))
+
+    #T_in = np.einsum('ijk,ijk->ijk',T,omega)
+    
+    return T,T_in,Omega
+
+
+def get_objective(T, U, V, W, omega, regParam):
+    t_obj = ctf.timer("ccd_get_objective")
+    t_obj.start()
+    L = ctf.tensor(T.shape, sp=T.sp)
+    t0 = time.time()
+    L.i("ijk") << T.i("ijk") - ctf.TTTP(omega, [U, V, W]).i("ijk")
+    t1 = time.time()
+    normL = ctf.vecnorm(L)
+    if T.sp:
+        RMSE = normL / (T.nnz_tot**.5)
+    else:
+        nnz_tot = ctf.sum(omega)
+        RMSE = normL / (nnz_tot**.5)
+    objective = normL + (ctf.vecnorm(U) + ctf.vecnorm(V) +
+                         ctf.vecnorm(W)) * regParam
+    t2 = time.time()
+    if glob_comm.rank() == 0:
+        print('generate L takes {}'.format(t1 - t0))
+        print('calc objective takes {}'.format(t2 - t1))
+    t_obj.stop()
+    return [objective, RMSE]
+
+def read_tensor_from_file(sp_frac):
+    #global netflix_tensor_dims
+    T = ctf.tensor((480189, 17770, 2182), sp=True)
+    T.read_from_file('/scratch/06720/tg860195/tensor.txt')
+
+    #T.fill_sp_random(1., 1., sp_frac)
+
+    #omega_np = np.ones(I*J*K)
+    #omega_np[:int(sp_frac*(I*J*K))] = 0
+    #np.random.shuffle(omega_np)
+    #omega = ctf.astensor(omega_np, dtype=np.float64)
+    #omega = omega.reshape(I, J, K)
+
+    #T_in = ctf.einsum('ijk,ijk->ijk',T,omega)
+
+    return T
+
+
+if __name__ == "__main__":
+    global netflix_tensor_dims
+
+    parser = argparse.ArgumentParser()
+    arg_defs.add_general_arguments(parser)
+    args, _ = parser.parse_known_args()
+
+    I = args.I
+    J = args.J
+    K = args.K
+    R = args.R
+
+    numiter_GN = args.num_iter_GN
+    numiter_ALS_imp = args.num_iter_ALS_implicit
+    time_limit = args.time_limit
+    err_thresh = args.err_thresh
+    sp_frac = args.sp_fraction
+    use_sp_rep = args.use_sparse_rep
+    block_size_ALS_imp = args.block_size_ALS_implicit
+    reg_ALS = args.regularization_ALS
+    use_func_tsr = args.function_tensor
+    
+    #T = read_tensor_from_file(sp_frac)
+    #omega = getOmega(T)
+    #I = 480189
+    #J = 17770
+    #K = 2182
+    T,T_in,omega= create_lowr_tensor(I, J, K, R, sp_frac, use_sp_rep)
+
+    
+
+    
+
+
+    ctf.random.seed(10)
+    U = ctf.random.random((I, R))
+    V = ctf.random.random((J, R))
+    W = ctf.random.random((K, R))
+
+
+    [_, RMSE] = get_objective(T, U, V, W, omega, 0)
+    if ctf.comm().rank() == 0:
+        print("Initial RMSE is ", RMSE)
+
+    if numiter_ALS_imp > 0:
+        if ctf.comm().rank() == 0:
+            print(
+                "Performing expicit ALS and regularization parameter is ",
+                reg_ALS)
+        U_copy = ctf.copy(U)
+        V_copy = ctf.copy(V)
+        W_copy = ctf.copy(W)
+
+        T_np = T.to_nparray()
+        T_in_np = T_in.to_nparray()
+        omega_np = omega.to_nparray()
+        U_np = U.to_nparray()
+        V_np = V.to_nparray()
+        W_np = W.to_nparray()
+        U_copy,V_copy,W_copy= explicit_als(tenpy_np,
+            T_in_np,
+            T_np,
+            omega_np,
+            U_np,
+            V_np,
+            W_np,
+            1e-04,
+            I,
+            J,
+            K,
+            R,
+            numiter_ALS_imp
+            )
+
+
+    if numiter_GN>0:
+        T_np = T.to_nparray()
+        T_in_np = T_in.to_nparray()
+        omega_np = omega.to_nparray()
+        U_np = U.to_nparray()
+        V_np = V.to_nparray()
+        W_np = W.to_nparray()
+        
+        U_np,V_np,W_np= getCPGN(tenpy_np,
+            T_in_np,
+            T_np,
+            omega_np,
+            U_np,
+            V_np,
+            W_np,
+            1e-04,
+            I,
+            J,
+            K,
+            R,
+            numiter_GN)