Add files via upload

Author: MWillettM

Demo.py

@@ -0,0 +1,29 @@
+import numpy as np
+import Precon_Conjugate_Gradient as PCG
+import Other_functions as wt
+import random
+
+
+#this demo version illustrates the increasing effectiveness of our preconditioner as the dimension of the underlying data grows. 
+#for our demo, we use the Quadratic Inverse RBF and draw from the multidimensional normal distribution, scaled appropriately. 
+
+n=1000
+
+for d in[2,3,6,10,25,100,158,1000,3981,10000,25118,63095]:
+    x_i = wt.points_normal(n,d)/(np.sqrt(2*d-1))
+    
+    f_i = [random.random() for _ in range(n)]
+
+    #Generating the interpolation matrix, phi.
+    phi = wt.interpolation_matrix_generator('Quadratic_inv',x_i)
+
+    #generating our preconditioner
+    precon_inv, const_1, const_2 = wt.precon_generator(n,1,'Quadratic_inv','')
+
+    #preconditioned solution
+    solution, count = PCG.untransformed_Pc_GC(phi,f_i,precon_inv,1e-10,const_1,const_2)
+
+    #unpreconditioned solution
+    solution2,count2 = PCG.CG_no_precon(phi,f_i,1e-10)
+
+    print('dimension:',d,'Preconditioned Iteration count:',count, 'unpreconditioned iteration count:',count2)

Full_implementation.py

@@ -0,0 +1,36 @@
+import Precon_Conjugate_Gradient as PCG
+import Other_functions as wt
+
+
+#this takes data centers, values at those points, the kind of RBF you wish to use and a value, x, at which you wish your interpolant to be evaluated.
+#this returns - iteration count for convergence of the CG algorithm, iterpolant coefficients and the interpolant evaluated at x.
+#it is suggested that you scale your data so that |x_i|<1 to prevent the RBFs from being too flat.
+#where appropriate, we use a shape parameter of 1 in the RBFs.
+
+def interpolant_generator_evaluator(data_centers, data_values, RBF,x):
+    x_i = [i for i in data_centers]
+    
+    n = len(x_i)
+
+    f_i = [i for i in data_values]
+
+    if n != len(f_i):
+        print('data_centers and data_values different length')
+        return
+    
+    else:
+
+        #Generating the interpolation matrix, phi.
+        phi = wt.interpolation_matrix_generator(RBF,x_i)
+
+        #generating our preconditioner
+        precon_inv, const_1, const_2 = wt.precon_generator(n,1,'Q','')
+
+        #generating our interpolation coefficients
+        solution, count = PCG.untransformed_Pc_GC(phi,f_i,precon_inv,1e-10,const_1,const_2)
+
+        #evaluating our interpolant at x
+        interpolant_eval = wt.interp(solution,x_i,x,RBF)
+
+        return count, solution, interpolant_eval
+    

Other_functions.py

@@ -0,0 +1,107 @@
+import numpy as np
+
+#Preconditioner generator - returns our preconditioner and its coefficients (the preconditioner is of the form aI + b e^Te)
+#n = number of data points
+#c = shape parameter for the RBFs - we are safe to use 1 here, but feel free to adjust to your liking.
+#distribution - if your data is uniformly distributed in the unit hypercube/unit sphere, select accordingly. Otherwise we take k=1.
+
+def precon_generator(n,c,RBF,distribution):
+    if distribution == 'cube':
+        k=1
+    if distribution =='ball':
+        k=np.sqrt(2)
+    else:
+        k=1
+
+    if RBF == 'Gaussian':
+        phi_rt2 = Gaussian(k,c)
+    if RBF == 'MQ':
+        phi_rt2 = MQ(k,c)
+    if RBF == 'MQ_inv':
+        phi_rt2 = MQ_inv(k,c)
+    if RBF == 'Quadratic_inv':
+        phi_rt2 = Quadratic_inv(k,c)
+    const_1 = phi_rt2 / ((phi_rt2-c)*(c-phi_rt2+n*phi_rt2))
+    const_2 = 1 / (phi_rt2-c)
+
+    preconditioner = const_1 * np.ones((n,n)) - const_2 * np.eye(n)
+    
+    return preconditioner, const_1, const_2
+
+#A function that evaluates our interpolant using 
+#coefficients - generated from the PCG algorithm output
+#centers - the same centers we used to generate the coefficients
+#x_value - the point we want to evaluate our interpolant at
+#RBF_type - make sure to use the same RBF type as the one used to generate the interpolant.
+
+#returns the value of the interpolant.
+
+def interp(coefficients, centers, x_value, RBF_type):
+    if RBF_type == 'Gaussian':
+        center_vector = [Gaussian(x_value - center,1) for center in centers]
+    if RBF_type == 'MQ':
+        center_vector = [MQ(x_value - center,1) for center in centers]
+    if RBF_type == 'MQ_inv':
+        center_vector = [MQ_inv(x_value - center,1) for center in centers]
+    if RBF_type == 'Quadratic_inv':
+        center_vector = [Quadratic_inv(x_value - center,1) for center in centers]
+    
+    return coefficients.T @ center_vector
+
+
+#RBF functions - note that our algorithm is not certain to converge for the multiquadric (MQ) as the interpolation matrix generated is not positive definite. 
+def MQ(point,c):
+    return np.sqrt(c**2+np.linalg.norm(point))
+
+def Gaussian(point, c):
+    return np.exp(-c * np.linalg.norm(point))
+
+def MQ_inv(point,c):
+    return 1/np.sqrt(c**2+np.linalg.norm(point))
+
+def Quadratic_inv(point,c):
+    return 1/(c**2+np.linalg.norm(point))
+
+#interpolation matrix generator
+def interpolation_matrix_generator(RBF_type, data):
+    phi = np.ones((len(data),len(data)))
+    for i, point in enumerate(data):
+        for j, center in enumerate(data):
+            phi[i,j] = np.linalg.norm(point-center)
+
+    if RBF_type == 'Gaussian':
+        vectorized_function = np.vectorize(Gaussian)
+    if RBF_type == 'MQ':
+        vectorized_function = np.vectorize(MQ)
+    if RBF_type == 'MQ_inv':
+        vectorized_function = np.vectorize(MQ_inv)
+    if RBF_type == 'Quadratic_inv':
+        vectorized_function = np.vectorize(Quadratic_inv)
+
+    return vectorized_function(phi,1)
+
+
+#data generator for a range of underlying distributions.
+def points_in_unit_ball(n,d):
+    cube = np.random.standard_normal(size=(n, d))
+    norms = np.linalg.norm(cube,axis=1)
+    surface_sphere = cube/norms[:,np.newaxis]
+    scales = np.random.uniform(0,1, size= n)
+    points = surface_sphere * (scales[:, np.newaxis])**(1/d)
+    return points
+
+def points_integer_grid(n,d, s_min, s_max):
+    points = np.random.randint(s_min,s_max, size = (n,d))
+    return points
+
+def points_in_cube(n,d,normalised):
+    points = np.random.uniform(0,1,size=(n,d))
+    if normalised == 'no':
+        return points
+    if normalised == 'yes':
+        return points/np.sqrt(d/6-17/120)
+
+def points_normal(n,d):
+    points = np.random.normal(0,1,size=(n,d))
+    return points
+

Precon_Conjugate_Gradient.py

@@ -0,0 +1,99 @@
+import numpy as np
+
+#Setup for 'Untransformed Preconditioned Conjugate Gradient Method' to solve
+#Ax=b with a preconditioner M, which gives us
+#inv(M)Ax = inv(M)b
+
+
+#A is the interpolation matrix
+#x is the solution
+#M is the preconditioner
+#precon_inv can be generated from Other_functions.precon_generator, or feel free to use your own - you will need to adjust the code for the action of your preconditioner on a vector.
+
+#See README for full description of the algorithm.
+
+#Returns the coefficients of the interpolant and the iteration count.
+
+def untransformed_Pc_GC(A,b,precon_inv,tolerance,c1,c2):
+    
+    #Setup
+    n=np.size(b)
+    initial_guess = np.zeros(n)
+    residual = b - A@initial_guess
+
+    #efficient calculation of our preconditioner's action on a vector
+    def precon_action_on_vector(const_1,const_2,vector):
+        vector_of_ones = np.ones(len(vector))
+        return const_1*np.sum(vector) * vector_of_ones - const_2 * vector
+
+    search_direction = precon_action_on_vector(c1,c2,residual)
+
+    #x will be our solution
+    x = initial_guess
+
+    dummy_variable = residual.T @ search_direction
+
+    #iteration count
+    k=0
+    #Implementation
+    while np.linalg.norm(residual) > tolerance:
+        k+=1
+        #dummy variable is used twice to avoid recalculation
+        
+        step_size = dummy_variable / (search_direction.T @ A @ search_direction)
+        x += step_size * search_direction
+
+        residual -= step_size * A@search_direction
+
+        new_dummy = residual.T @ precon_action_on_vector(c1,c2,residual)
+
+        scale_factor = new_dummy / dummy_variable
+
+        dummy_variable = new_dummy
+
+        search_direction = precon_action_on_vector(c1,c2,residual)+ scale_factor*search_direction
+
+        if k> 1000:
+            break
+
+    return x,k
+
+#Unpreconditioned algorithm for comparison
+
+def CG_no_precon(A,b,tolerance):
+
+    #Setup
+    n=np.size(b)
+    initial_guess = np.zeros(n)
+    residual = b - A@initial_guess
+    search_direction = residual
+
+    #x will be our solution
+    x = initial_guess
+
+    dummy_variable = residual.T @ residual
+
+    #iteration coun
+    k=0
+    #Implementation
+    while np.linalg.norm(residual) > tolerance:
+        k+=1
+        #dummy variable is used twice to avoid recalculation
+
+        step_size = dummy_variable / (search_direction.T @ A @ search_direction)
+        x += step_size * search_direction
+
+        residual -= step_size * A@search_direction
+
+        new_dummy = residual.T @ residual
+
+        scale_factor = new_dummy / dummy_variable
+
+        dummy_variable = new_dummy
+
+        search_direction = residual + scale_factor*search_direction
+
+        if k>1000:
+            break
+
+    return x,k