#! /usr/bin/env python2 # -*- coding: utf-8; mode: python -*- """ A Python 2 script to test kernels, for the Kernel Methods homework. Source : http://lear.inrialpes.fr/people/mairal/teaching/2015-2016/MVA/ PDF: http://lear.inrialpes.fr/people/mairal/teaching/2015-2016/MVA/fichiers/homework_mva_2016.pdf - *Date:* Saturday 06 February 2016. - *Author:* Lilian Besson, for the MVA Master, (C) 2015-16. - *Licence:* MIT Licence (http://lbesson.mit-license.org). """ from __future__ import print_function, division # Python 2 compatibility if needed import numpy as np try: from sympy import lcm, gcd except: # if True: print("Error: sympy is not available, switching back to fractions implementation of lcm/gcd") from fractions import gcd def lcm(x, y): """ lcm(x, y) : lowest common multiple, = x * y / gcd(x, y). """ if gcd(x, y) == 0: return 0 else: return (x * y) / gcd(x, y) #: Number of random tests to do N = 1000 #: Size of the sum_ij T = 5 T = 50 #: Lower bound for random number in IR LOWER_BOUND = -300.0 #: Upper bound for random number in IR UPPER_BOUND = 300.0 #: Upper bound for random number in IN UPPER_BOUND_INT = 1000 #: Small upper bound for random number in IN SMALL_UPPER_BOUND_INT = 30 #: Tolerance for real values tests to zero REAL_TOLERANCE = 1e-12 #: Tolerance for probabilities check to 100% PROBA_TOLERANCE = 0.01 #: Standard deviation of the weights WEIGHT_STD = 100 #: Distribution for weights # WEIGHT_DISTRIBUTION = "NORMAL" # WEIGHT_DISTRIBUTION = "UNIFORM_0_1" # Useless! WEIGHT_DISTRIBUTION = "UNIFORM_-1_1" #: Map a domain to a random generator on this domain random_generators = { "REAL": lambda: np.random.uniform(LOWER_BOUND, UPPER_BOUND), "REALPLUS": lambda: np.random.uniform(0.0, UPPER_BOUND), "INTEGER": lambda: np.random.randint(0, UPPER_BOUND_INT), "SMALLINT": lambda: np.random.randint(0, SMALL_UPPER_BOUND_INT), "MONEONE": lambda: np.random.uniform(-1.0, 1.0) } def test_kernel(K, N=N, T=T, verb=False): """ Test if a function is a kernel on lots of random values (500 by default): - test if K(x, y) = K(y, x) (symmetry), - test if K(x, x) = 0 ==> x = 0 (definiteness), - test if sum_ij a_i a_j K(x_i, x_j) >= 0 for 'all' t, and a_i in IR, x_i in DOMAIN (positiveness) (for all i,j in range(t), and t drawn in range(T) randomly). Return numerical values, average of success rate on N = 500 tests. """ # Get the domain of K and an appropriate random generator. domain = K.domain newvalue = random_generators[domain] # Random generator for weights def newweights(t): if WEIGHT_DISTRIBUTION == "UNIFORM_0_1": return np.random.rand(t) * WEIGHT_STD if WEIGHT_DISTRIBUTION == "UNIFORM_-1_1": return (2*np.random.rand(t) - 1) * (WEIGHT_STD / 2) if WEIGHT_DISTRIBUTION == "NORMAL": return WEIGHT_STD * np.random.randn(t) # Test for symmetry prop_sym = 0 # Test for definiteness prop_def = 0 # Test for positiveness prop_pos = 0 for i in range(N): x = newvalue() y = newvalue() if verb: print(" - Random values x = {}, y = {}.".format(x, y)) if K(x, y) == K(y, x): prop_sym += 1 # Test : (A ==> B) <=> (B or not(A)) if (not (abs(K(x, x)) < REAL_TOLERANCE)) or (abs(x) > REAL_TOLERANCE): prop_def += 1 if (not (abs(K(y, y)) < REAL_TOLERANCE)) or (abs(y) > REAL_TOLERANCE): prop_def += 1 # Generate lots of samples for checking positiveness t = np.random.randint(1, T) X = [newvalue() for i in range(t)] a = newweights(t) if verb: print(" - Random values t = {}, a = {}, X = {}.".format(t, a, X)) if sum(a[i] * a[j] * K(X[i], X[j]) for i in range(t) for j in range(t)) >= 0: prop_pos += 1 if verb: print(" - After {} random tests, sym = {:>3}, def = {:>3}, pos = {:>3}.".format(N, prop_sym, prop_def, prop_pos)) # Return the results return { "symmetry": (prop_sym / N), "definiteness": (prop_def / (2*N)), "positiveness": (prop_pos / N), } # %% Tweak to define a new kernel #: List of all kernels to test, dynamically created by @kernel decorator all_kernels = [] def kernel(domain="REAL"): """ Decorator for a new kernel (just an elegant way to add it to all_kernels). The parameter 'domain' can be: - "REAL" for (-oo,+oo) = IR real number, - "REALPLUS" for [0,+oo) = IR+ positive real number, - "INTEGER" for [|0,+oo) = IN positive integer, - "SMALLINT" for [|0,+oo) = IN (small) positive integer, - "MONEONE" for (-1,+1) = UR the unit open disk of real values. """ # global all_kernels def kernel_decorator(func): """ Decorator for a new kernel, with a certain domain. """ global all_kernels # Add a new kernel ! all_kernels += [func] func.domain = domain return func return kernel_decorator # %% List of Kernels # Classical example of kernel: linear kernel! @kernel("REAL") def abs_k(x, y): """ K(x, y) = < x, y > = x * y """ return x * y @kernel("MONEONE") def one_by_k(x, y): """ K(x, y) = 1 / (1 - x * y) """ return 1 / (1 - x * y) # Here we could have an overflow error, so using small integers # From the slides @kernel("SMALLINT") def power_two_dot_k(x, y): """ K(x, y) = 2 ** (x + y) """ return 2 ** (x * y) # Here we could have an overflow error, so using small integers @kernel("SMALLINT") def power_two_sum_k(x, y): """ K(x, y) = 2 ** (x * y) """ return 2 ** (x * y) @kernel("REALPLUS") def log_k(x, y): """ K(x, y) = log(1 + x * y) """ return np.log(1 + x * y) @kernel("REAL") def exp_k(x, y): """ K(x, y) = exp(-(x-y)^2) """ return np.exp(-(x-y)**2) @kernel("REAL") def cos_plus_k(x, y): """ K(x, y) = cos(x + y) """ return np.cos(x + y) @kernel("REAL") def cos_minus_k(x, y): """ K(x, y) = cos(x - y) """ return np.cos(x - y) @kernel("REALPLUS") def min_k(x, y): """ K(x, y) = min(x, y) """ return min(x, y) @kernel("REALPLUS") def max_k(x, y): """ K(x, y) = max(x, y) """ return max(x, y) @kernel("REALPLUS") def min_by_max_k(x, y): """ K(x, y) = min(x, y) / max(x, y) """ if abs(max(x, y)) < REAL_TOLERANCE: return 0.0 else: return min(x, y) / max(x, y) @kernel("INTEGER") def gcd_k(x, y): """ K(x, y) = gcd(x, y) """ return gcd(x, y) @kernel("INTEGER") def lcm_k(x, y): """ K(x, y) = lcm(x, y) """ return lcm(x, y) @kernel("INTEGER") def gcd_by_lcm_k(x, y): """ K(x, y) = gcd(x, y) / lcm(x, y) """ if abs(lcm(x, y)) < REAL_TOLERANCE: return 0.0 else: return gcd(x, y) / lcm(x, y) # %% Testing all the kernels def main(kernels_to_test=all_kernels, T=T, verb=False): """ Test a lot of different kernels. """ nb = len(all_kernels) print("\nStarting to test {} kernels...".format(nb)) # Store results results = [1] * nb not_kernels = [] # Start for i, K in enumerate(all_kernels): print("\n- Testing the kernel {} {}:".format(K.__name__, K.__doc__)) print(" It has domain = {}.".format(K.domain)) result = test_kernel(K, N=N, T=T, verb=verb) print(" Results of the N = {} tests:".format(N)) # bad_test = None # min_rate = None for test, rate in result.items(): print(" - Tested for {:^12} : success rate {:.2%}.".format(test, rate)) if rate < 1.0 - PROBA_TOLERANCE: print(" ===> I guess this kernel {} {} is NOT a positive definite kernel !".format(K.__name__, K.__doc__)) not_kernels.append({'K': K, 'test': test, 'rate': rate}) # bad_test, min_rate = test, rate print(" Finished for the kernel {}...".format(K.__name__)) results[i] = result print("\n\nDone testing some kernels...") return results, not_kernels if __name__ == '__main__': print("\nScript to test if functions are kernels, by Lilian Besson (C) 2016") print("MIT Licensed (http://lbesson.mit-license.org)") print("\n\nUsing N = {}, T = {}.".format(N, T)) print("For test domains:") print(" - IR ~= [{}, {}],".format(LOWER_BOUND, UPPER_BOUND)) print(" - IR+ ~= [ 0, {}],".format(UPPER_BOUND)) print(" - IN ~= [ 0, {}],".format(UPPER_BOUND_INT)) print(" - IN ~= [ 0, {}] (small to reduce overflow),".format(SMALL_UPPER_BOUND_INT)) print("Tolerance for small real values = {}.".format(REAL_TOLERANCE)) print("Tolerance for success rates = {}.".format(PROBA_TOLERANCE)) if WEIGHT_DISTRIBUTION == "UNIFORM_0_1": print("And sampling the a_i from uniform in [0, {}].".format(WEIGHT_STD)) if WEIGHT_DISTRIBUTION == "UNIFORM_-1_1": print("And sampling the a_i from uniform in [-{}, {}].".format(WEIGHT_STD/2, WEIGHT_STD/2)) if WEIGHT_DISTRIBUTION == "NORMAL": print("And sampling the a_i from a Gaussian of mean 0 and std {}.".format(WEIGHT_STD)) # all_kernels = [absk] # results, not_kernels = main(kernels_to_test=all_kernels, verb=True) results, not_kernels = main() for item in not_kernels: K = item['K'] test = item['test'] rate = item['rate'] print("\n- Apparently, the kernel {} {} is NOT a positive definite kernel : because {} was only observed with probability {:.2%}.".format(K.__name__, K.__doc__, test, rate)) print("\n- All the other kernel appeared to be positive definite. Now prove it.") print("\nFinished for this script 'test_kernel.py'.") # End of test_kernel.py