functions.py

import numpy as np
from scipy.sparse import csc_matrix

"""# Objective Functions

"""
#Basis generator.
#Generates a set of n-orthonormal vectors.
def rvs(dim=3):
     random_state = np.random
     H = np.eye(dim)
     D = np.ones((dim,))
     for n in range(1, dim):
         x = random_state.normal(size=(dim-n+1,))
         D[n-1] = np.sign(x[0])
         x[0] -= D[n-1]*np.sqrt((x*x).sum())
         # Householder transformation
         Hx = (np.eye(dim-n+1) - 2.*np.outer(x, x)/(x*x).sum())
         mat = np.eye(dim)
         mat[n-1:, n-1:] = Hx
         H = np.dot(H, mat)
         # Fix the last sign such that the determinant is 1
     D[-1] = (-1)**(1-(dim % 2))*D.prod()
     # Equivalent to np.dot(np.diag(D), H) but faster, apparently
     H = (D*H.T).T
     return H

#Generate a random PSD quadratic with eigenvalues between certain numbers.
def randomPSDGenerator(dim, Mu, L):
    eigenval = np.zeros(dim)
    eigenval[0] = Mu
    eigenval[-1] = L
    eigenval[1:-1] = np.random.uniform(Mu, L, dim - 2)
    M = np.zeros((dim, dim))
    A = rvs(dim)
    for i in range(dim):
        M += eigenval[i]*np.outer(A[i], A[i])
    return M

#Random PSD matrix with a given sparsity.
def randomPSDGeneratorSparse(dim, sparsity):
    mask = np.random.rand(dim,dim)> (1- sparsity)
    mat = np.random.normal(size = (dim,dim))
    Aux = np.multiply(mat, mask)
    return np.dot(Aux.T, Aux) + np.identity(dim)

def calculateEigenvalues(M):
    from scipy.linalg import eigvalsh
    dim = len(M)
    L = eigvalsh(M, eigvals = (dim - 1,dim - 1))[0]
    Mu = eigvalsh(M, eigvals = (0,0))[0]
    return L, Mu
    
#Takes a random PSD matrix generated by the functions above and uses them as a function.
class funcQuadratic:
    import numpy as np
    def __init__(self, size, matrix, vector, Mu, L):
        self.len = size
        self.M = matrix.copy()
        self.b = vector.copy()
        self.L = L
        self.Mu = Mu
        return       
           
    def dim(self):
        return self.len
        
    #Evaluate function.
    def fEval(self, x):
        return 0.5*np.dot(x, self.M.dot(x)) + np.dot(self.b, x) 
    
    #Evaluate gradient.
    def fEvalGrad(self, x):
        return self.M.dot(x) + self.b
    
    #Line Search.
    def lineSearch(self, grad, d):
        return -np.dot(grad, d)/np.dot(d, self.M.dot(d))
      
    #Return largest eigenvalue.
    def largestEig(self):
        return self.L

    #Return smallest eigenvalue.
    def smallestEig(self):
        return self.Mu
    
    #Return largest eigenvalue.
    def returnM(self):
        return self.M

    #Return smallest eigenvalue.
    def returnb(self):
        return self.b

class funcQuadraticDiag:
    import numpy as np
    def __init__(self, size, xOpt, Mu = 1.0, L = 2.0):
        self.len = size
        self.matdim = int(np.sqrt(size))
        self.eigenval = np.zeros(size)
        self.eigenval[0] = Mu
        self.eigenval[-1] = L
        self.eigenval[1:-1] = np.random.uniform(Mu, L, size - 2)
        self.L = L
        self.Mu = Mu
        self.xOpt = xOpt
        self.b = - np.multiply(self.xOpt, self.eigenval)
        return       
           
    def dim(self):
        return self.len
        
    #Evaluate function.
    def fEval(self, x):
        return 0.5*np.dot(x, np.multiply(self.eigenval,x)) + np.dot(self.b, x) 
    
    #Evaluate gradient.
    def fEvalGrad(self, x):
        return np.multiply(x, self.eigenval) + self.b
      
    #Return largest eigenvalue.
    def largestEig(self):
        return self.L

    #Return smallest eigenvalue.
    def smallestEig(self):
        return self.Mu
    
    #Line Search.
    def lineSearch(self, grad, d):
        return -np.dot(grad, d)/np.dot(d, np.multiply(self.eigenval, d))
    
        #Return largest eigenvalue.
    def returnM(self):
        return self.eigenval

    #Return smallest eigenvalue.
    def returnb(self):
        return self.b

#Function used in NAGD
class funcSimplexLambdaNormalizedEigen:
    #Assemble the matrix from the active set.
    def __init__(self, activeSet, z, A, L, Mu):
        from scipy.sparse.linalg import eigsh
        self.len = len(activeSet)
        Mat = np.zeros((len(activeSet[0]), self.len))
        self.c = Mu*A + L - Mu
        self.b = np.zeros(len(activeSet))
        for i in range(0, self.len):
            Mat[:, i] = activeSet[i]
            self.b[i] = -np.dot(z, activeSet[i])
        self.b /= self.c
        self.M = np.dot(np.transpose(Mat),Mat)
        #Create a sparse matrix from the data.
        self.M = csc_matrix(self.M)
        if(self.M.shape == (1,1)):
            self.L = 1.0
            self.Mu = 1.0
        else:
            self.L = eigsh(self.M, 1, which='LM', return_eigenvectors = False)[0]
            self.Mu = eigsh(self.M, 1, sigma=1.0e-10, which='LM', return_eigenvectors = False)[0]
        return
        
    def fEval(self, x):
        return 0.5*np.dot(x.T, self.M.dot(x)) + np.dot(self.b, x)
    
    def fEvalGrad(self, x):
        return self.M.dot(x) + self.b
    
    #Line Search.
    def lineSearch(self, grad, d):
        return -np.dot(grad, d)/np.dot(d, np.dot(self.M, d))
    
    def returnM(self):
        return self.M
    
    def largestEig(self):
        return self.L
    
    def smallestEig(self):
        return self.Mu
    
    def FWGap(self, x):
        grad = self.fEvalGrad(x)
        v = np.zeros(len(x))
        minVert = np.argmin(grad)
        v[minVert] = 1.0
        return np.dot(grad, x - v)
    
    def returnlen(self):
        return self.len
    
    def update(self, activeSet, z, A, L, Mu):
        self.c = Mu*A + L - Mu
        self.b = np.zeros(len(activeSet))
        for i in range(0, len(activeSet)):
            self.b[i] = -np.dot(z, activeSet[i])
        self.b /= self.c
        return
        
#Function that will be used in the optimization problem.
#Creates a random quadratic with a positive semidefinite matrix.
class funcAccelScheme:
    import numpy as np
    #Input the initial point, the dimension of the set and the Q matrix.
    #Check that the initial matrix is positive semidefinite.
    def __init__(self, size, matrix, vector, L, Mu):
        self.len = size
        self.matdim = int(np.sqrt(size))
        self.M = matrix.copy()
        self.y = np.zeros(size, dtype = float)
        self.kappa = 0.0
        self.b = vector.copy()
        self.L = L
        self.Mu = Mu
        
    def sety(self, yVal):
        self.y = yVal
        return
    
    def setKappa(self, kappaVal):
        self.kappa = kappaVal
        return
        
    def fEval(self, x):
        return 0.5*np.dot(x, self.M.dot(x)) + np.dot(self.b, x) + self.kappa/2 * np.dot(x - self.y, x - self.y)

    def fEvalBaseProblem(self, x):
        return 0.5*np.dot(x, self.M.dot(x)) + np.dot(self.b, x)
        
    def fEvalGrad(self, x):
        return self.M.dot(x) + self.b + self.kappa * (x - self.y)
    
    def fEvalGradBaseProblem(self, x):
        return self.M.dot(x) + self.b
    
    #Line Search.
    def lineSearch(self, grad, d):
        return -np.dot(grad, d)/(np.dot(d, self.M.dot(d)) + self.kappa*np.dot(d,d))
    
    def FWGap(self, x, feasReg):
        grad = self.fEvalGrad(x)
        v = feasReg.LPOracle(grad)
        return np.dot(grad, (x - v))
    
    def FWGapBaseProblem(self, x, feasReg):
        grad = self.fEvalGradBaseProblem(x)
        v = feasReg.LPOracle(grad)
        return np.dot(grad, (x - v))
    
    def largestEig(self):
        return self.L

    def smallestEig(self):
        return self.Mu

#Function that will be used in the optimization problem.
#Creates a random quadratic with a positive semidefinite matrix.
class funcAccelSchemeDiag:
    import numpy as np
    #Input the initial point, the dimension of the set and the Q matrix.
    #Check that the initial matrix is positive semidefinite.
    def __init__(self, size, matrix, vector, L, Mu):
        self.len = size
        self.matdim = int(np.sqrt(size))
        self.M = matrix
        self.y = np.zeros(size, dtype = float)
        self.kappa = 0.0
        self.b = vector
        self.L = L
        self.Mu = Mu
        
    def sety(self, yVal):
        self.y = yVal
        return
    
    def setKappa(self, kappaVal):
        self.kappa = kappaVal
        return
        
    def fEval(self, x):
        return 0.5*np.dot(x,np.multiply(self.M,x)) + np.dot(self.b, x) + self.kappa/2 * np.dot(x - self.y, x - self.y)

    def fEvalBaseProblem(self, x):
        return 0.5*np.dot(x,np.multiply(self.M,x)) + np.dot(self.b, x)
        
    def fEvalGrad(self, x):
        return np.multiply(x, self.M) + self.b + self.kappa * (x - self.y)
    
    def fEvalGradBaseProblem(self, x):
        return np.multiply(x, self.M) + self.b

    def FWGap(self, x, feasReg):
        grad = self.fEvalGrad(x)
        v = feasReg.LPOracle(grad)
        return np.dot(grad, (x - v))
    
    def FWGapBaseProblem(self, x, feasReg):
        grad = self.fEvalGradBaseProblem(x)
        v = feasReg.LPOracle(grad)
        return np.dot(grad, (x - v))
    
    #Line Search.
    def lineSearch(self, grad, d):
        return -np.dot(grad, d)/np.dot(d, np.multiply(self.M, d))
    
    def largestEig(self):
        return self.L

    def smallestEig(self):
        return self.Mu