Merge pull request #20 from lcmd-epfl/regression-new-features

Add new features to regression tools

- sparse regression
- random state, read kernel, return prediction arguments

Author: briling

qstack/math/fps.py

@@ -0,0 +1,18 @@
+import numpy as np
+
+def do_fps(x, d=0):
+    # Code from Giulio Imbalzano
+    n = len(x)
+    if d==0:
+        d = n
+    iy = np.zeros(d,int)
+    measure = np.zeros(d-1,float)
+    iy[0] = 0
+    # Faster evaluation of Euclidean distance
+    n2 = np.sum(x*x, axis=1)
+    dl = n2 + n2[iy[0]] - 2.0*np.dot(x,x[iy[0]])
+    for i in range(1,d):
+        iy[i], measure[i-1] = np.argmax(dl), np.amax(dl)
+        nd = n2 + n2[iy[i]] - 2.0*np.dot(x,x[iy[i]])
+        dl = np.minimum(dl,nd)
+    return iy, measure

qstack/regression/final_error.py

@@ -5,22 +5,34 @@
 from qstack.regression.kernel_utils import get_kernel, defaults
 from qstack.tools import correct_num_threads
 
-def final_error(X, y, sigma=defaults.sigma, eta=defaults.eta, akernel=defaults.kernel, test_size=defaults.test_size, save_alpha=None):
+def final_error(X, y, read_kernel=False, sigma=defaults.sigma, eta=defaults.eta, akernel=defaults.kernel,
+                test_size=defaults.test_size,
+                random_state=defaults.random_state,
+                return_pred=False, save_alpha=None):
     """
 
     .. todo::
         Write the docstring
     """
-    kernel = get_kernel(akernel)
-    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_size, random_state=0)
-    K_all  = kernel(X_train, X_train, 1.0/sigma)
-    Ks_all = kernel(X_test,  X_train, 1.0/sigma)
+    if read_kernel is False:
+        kernel = get_kernel(akernel)
+        X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_size, random_state=random_state)
+        K_all  = kernel(X_train, X_train, 1.0/sigma)
+        Ks_all = kernel(X_test,  X_train, 1.0/sigma)
+    else:
+        idx_train, idx_test, y_train, y_test = train_test_split(np.arange(len(y)), y, test_size=test_size, random_state=random_state)
+        K_all  = X[np.ix_(idx_train,idx_train)]
+        Ks_all = X[np.ix_(idx_test, idx_train)]
     K_all[np.diag_indices_from(K_all)] += eta
     alpha = scipy.linalg.solve(K_all, y_train, assume_a='pos')
     y_kf_predict = np.dot(Ks_all, alpha)
     aes = np.abs(y_test-y_kf_predict)
-    if save_alpha: np.save(save_alpha, alpha)
-    return aes
+    if save_alpha:
+        np.save(save_alpha, alpha)
+    if return_pred:
+        return aes, y_kf_predict
+    else:
+        return aes
 
 def main():
     import sys
@@ -34,12 +46,13 @@ def main():
     parser.add_argument('--kernel',     type=str,   dest='kernel',      default=defaults.kernel,    help='kernel type (G for Gaussian, L for Laplacian, myL for Laplacian for open-shell systems) (default '+defaults.kernel+')')
     parser.add_argument('--save-alpha', type=str,   dest='save_alpha',  default=None,               help='file to write the regression coefficients to (default None)')
     parser.add_argument('--ll',     action='store_true', dest='ll',     default=False,              help='if correct for the numper of threads')
+    parser.add_argument('--random_state',  type=int, dest='random_state', default=defaults.random_state,  help='random state for test / train splitting')
     args = parser.parse_args()
     print(vars(args))
     if(args.ll): correct_num_threads()
     X = np.load(args.repr)
     y = np.loadtxt(args.prop)
-    aes = final_error(X, y, sigma=args.sigma, eta=args.eta, akernel=args.kernel, test_size=args.test_size, save_alpha=args.save_alpha)
+    aes = final_error(X, y, sigma=args.sigma, eta=args.eta, akernel=args.kernel, test_size=args.test_size, save_alpha=args.save_alpha, random_state=random_state)
     np.savetxt(sys.stdout, aes, fmt='%e')
 
 if __name__ == "__main__":

qstack/regression/hyperparameters.py

@@ -6,26 +6,38 @@
 from sklearn.model_selection import train_test_split, KFold
 from qstack.regression.kernel_utils import get_kernel, defaults, ParseKwargs
 from qstack.tools import correct_num_threads
+from qstack.math.fps import do_fps
 
 def hyperparameters(X, y,
            sigma=defaults.sigmaarr, eta=defaults.etaarr, gkernel=defaults.gkernel, gdict=defaults.gdict,
            akernel=defaults.kernel, test_size=defaults.test_size, splits=defaults.splits,
-           printlevel=0, adaptive=False, read_kernel=False):
+           printlevel=0, adaptive=False, read_kernel=False, sparse=None):
     """
 
     .. todo::
         Write the docstring
     """
 
-    def k_fold_opt(K_all):
+    def k_fold_opt(K_all, eta):
         kfold = KFold(n_splits=splits, shuffle=False)
         all_maes = []
         for train_idx, test_idx in kfold.split(X_train):
             y_kf_train, y_kf_test = y_train[train_idx], y_train[test_idx]
-            K  = K_all [np.ix_(train_idx,train_idx)]
-            Ks = K_all [np.ix_(test_idx,train_idx)]
+
+            if not sparse:
+                K_solve = np.copy(K_all [np.ix_(train_idx,train_idx)])
+                K_solve[np.diag_indices_from(K_solve)] += eta
+                y_solve = y_kf_train
+                Ks = K_all [np.ix_(test_idx,train_idx)]
+            else:
+                K_NM    = K_all [np.ix_(train_idx,sparse_idx)]
+                K_solve = K_NM.T @ K_NM
+                K_solve[np.diag_indices_from(K_solve)] += eta
+                y_solve = K_NM.T @ y_kf_train
+                Ks = K_all [np.ix_(test_idx,sparse_idx)]
+
             try:
-                alpha = scipy.linalg.solve(K, y_kf_train, assume_a='pos', overwrite_a=True)
+                alpha = scipy.linalg.solve(K_solve, y_solve, assume_a='pos', overwrite_a=True)
             except scipy.linalg.LinAlgError:
                 print('singular matrix')
                 all_maes.append(np.nan)
@@ -43,9 +55,7 @@ def hyper_loop(sigma, eta):
                 K_all = X_train
 
             for e in eta:
-                K_all[np.diag_indices_from(K_all)] += e
-                mean, std = k_fold_opt(K_all)
-                K_all[np.diag_indices_from(K_all)] -= e
+                mean, std = k_fold_opt(K_all, e)
                 if printlevel>0 :
                     sys.stderr.flush()
                     print(s, e, mean, std, flush=True)
@@ -63,6 +73,11 @@ def hyper_loop(sigma, eta):
         X_train = X[np.ix_(idx_train,idx_train)]
         sigma = [np.nan]
 
+    if sparse:
+        if read_kernel:
+            raise RuntimeError('Cannot do FPS with kernels')
+        sparse_idx = do_fps(X_train)[0][:sparse]
+
     work_sigma = sigma
     errors = []
     direction = None
@@ -111,14 +126,16 @@ def main():
     parser.add_argument('--ll',   action='store_true', dest='ll',       default=False,  help='if correct for the numper of threads')
     parser.add_argument('--ada',  action='store_true', dest='adaptive', default=False,  help='if adapt sigma')
     parser.add_argument('--readkernel', action='store_true', dest='readk', default=False,  help='if X is kernel')
+    parser.add_argument('--sparse',     type=int, dest='sparse', default=None,  help='regression basis size for sparse learning')
     args = parser.parse_args()
     if(args.readk): args.sigma = [np.nan]
     print(vars(args))
     if(args.ll): correct_num_threads()
 
     X = np.load(args.repr)
     y = np.loadtxt(args.prop)
-    errors = hyperparameters(X, y, read_kernel=args.readk, sigma=args.sigma, eta=args.eta, akernel=args.akernel, test_size=args.test_size, splits=args.splits, printlevel=args.printlevel, adaptive=args.adaptive)
+    errors = hyperparameters(X, y, read_kernel=args.readk, sigma=args.sigma, eta=args.eta, akernel=args.akernel, sparse=args.sparse,
+                             test_size=args.test_size, splits=args.splits, printlevel=args.printlevel, adaptive=args.adaptive)
 
     print()
     print('error        stdev          eta          sigma')

qstack/regression/kernel_utils.py

@@ -30,7 +30,8 @@ def __call__(self, parser, namespace, values, option_string=None):
   train_size=[0.125, 0.25, 0.5, 0.75, 1.0],
   etaarr=list(numpy.logspace(-10, 0, 5)),
   sigmaarr=list(numpy.logspace(0,6, 13)),
-  sigmaarr_mult=list(numpy.logspace(0,2, 5))
+  sigmaarr_mult=list(numpy.logspace(0,2, 5)),
+  random_state=0,
   )
 
 

qstack/regression/regression.py

@@ -5,25 +5,36 @@
 from sklearn.model_selection import train_test_split
 from qstack.regression.kernel_utils import get_kernel, defaults, ParseKwargs
 from qstack.tools import correct_num_threads
+from qstack.math.fps import do_fps
 
-def regression(X, y, read_kernel=False, sigma=defaults.sigma, eta=defaults.eta, akernel=defaults.kernel, gkernel=defaults.gkernel, gdict=defaults.gdict, test_size=defaults.test_size, train_size=defaults.train_size, n_rep=defaults.n_rep, debug=False):
+
+def regression(X, y, read_kernel=False, sigma=defaults.sigma, eta=defaults.eta,
+               akernel=defaults.kernel, gkernel=defaults.gkernel, gdict=defaults.gdict,
+               test_size=defaults.test_size, train_size=defaults.train_size, n_rep=defaults.n_rep,
+               random_state=defaults.random_state,
+               sparse=None, debug=False):
     """
-    
+
     .. todo::
         Write the docstring
     """
     if read_kernel is False:
         kernel = get_kernel(akernel, [gkernel, gdict])
-        X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_size, random_state=0)
+        X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_size, random_state=random_state)
         K_all  = kernel(X_train, X_train, 1.0/sigma)
         Ks_all = kernel(X_test,  X_train, 1.0/sigma)
     else:
-        idx_train, idx_test, y_train, y_test = train_test_split(np.arange(len(y)), y, test_size=test_size, random_state=0)
+        idx_train, idx_test, y_train, y_test = train_test_split(np.arange(len(y)), y, test_size=test_size, random_state=random_state)
         K_all  = X[np.ix_(idx_train,idx_train)]
         Ks_all = X[np.ix_(idx_test, idx_train)]
 
-    K_all[np.diag_indices_from(K_all)] += eta
     all_indices_train = np.arange(len(y_train))
+    if not sparse:
+        K_all[np.diag_indices_from(K_all)] += eta
+    else:
+        if read_kernel:
+            raise RuntimeError('Cannot do FPS with kernels')
+        sparse_idx = do_fps(X_train)[0][:sparse]
 
     if debug:
         np.random.seed(666)
@@ -35,40 +46,55 @@ def regression(X, y, read_kernel=False, sigma=defaults.sigma, eta=defaults.eta,
         for rep in range(n_rep):
             train_idx = np.random.choice(all_indices_train, size = size_train, replace=False)
             y_kf_train = y_train[train_idx]
-            K  = K_all [np.ix_(train_idx,train_idx)]
-            Ks = Ks_all[:,train_idx]
-            alpha = scipy.linalg.solve(K, y_kf_train, assume_a='pos')
+
+            if not sparse:
+                K_solve = K_all [np.ix_(train_idx,train_idx)]
+                y_solve = y_kf_train
+                Ks = Ks_all[:,train_idx]
+            else:
+                K_NM    = K_all [np.ix_(train_idx,sparse_idx)]
+                K_solve = K_NM.T @ K_NM
+                K_solve[np.diag_indices_from(K_solve)] += eta
+                y_solve = K_NM.T @ y_kf_train
+                Ks = Ks_all[:,sparse_idx]
+
+            alpha = scipy.linalg.solve(K_solve, y_solve, assume_a='pos')
             y_kf_predict = np.dot(Ks, alpha)
             maes.append(np.mean(np.abs(y_test-y_kf_predict)))
+
         maes_all.append((size_train, np.mean(maes), np.std(maes)))
     return maes_all
 
+
 def main():
     import argparse
     parser = argparse.ArgumentParser(description='This program computes the learning curve.')
-    parser.add_argument('--x',          type=str,   dest='repr',       required=True, help='path to the representations file')
-    parser.add_argument('--y',          type=str,   dest='prop',       required=True, help='path to the properties file')
-    parser.add_argument('--test',       type=float, dest='test_size',  default=defaults.test_size, help='test set fraction (default='+str(defaults.test_size)+')')
-    parser.add_argument('--eta',        type=float, dest='eta',        default=defaults.eta,       help='eta hyperparameter (default='+str(defaults.eta)+')')
-    parser.add_argument('--sigma',      type=float, dest='sigma',      default=defaults.sigma,     help='sigma hyperparameter (default='+str(defaults.sigma)+')')
-    parser.add_argument('--akernel',     type=str,   dest='akernel',     default=defaults.kernel,    help='local kernel type (G for Gaussian, L for Laplacian, myL for Laplacian for open-shell systems) (default '+defaults.kernel+')')
-    parser.add_argument('--gkernel',     type=str,   dest='gkernel',     default=defaults.gkernel,    help='global kernel type (avg for average kernel, rem for REMatch kernel) (default '+str(defaults.gkernel)+')')
-    parser.add_argument('--gdict',     nargs='*',   action=ParseKwargs, dest='gdict',     default=defaults.gdict,    help='dictionary like input string to initialize global kernel parameters')
-    parser.add_argument('--splits',     type=int,   dest='splits',     default=defaults.n_rep,     help='number of splits (default='+str(defaults.n_rep)+')')
-    parser.add_argument('--train',      type=float, dest='train_size', default=defaults.train_size, nargs='+', help='training set fractions')
-    parser.add_argument('--debug',      action='store_true', dest='debug', default=False,  help='enable debug')
-    parser.add_argument('--ll',         action='store_true', dest='ll',    default=False,  help='if correct for the numper of threads')
-    parser.add_argument('--readkernel', action='store_true', dest='readk', default=False,  help='if X is kernel')
+    parser.add_argument('--x',             type=str,   dest='repr',       required=True, help='path to the representations file')
+    parser.add_argument('--y',             type=str,   dest='prop',       required=True, help='path to the properties file')
+    parser.add_argument('--test',          type=float, dest='test_size',  default=defaults.test_size, help='test set fraction (default='+str(defaults.test_size)+')')
+    parser.add_argument('--eta',           type=float, dest='eta',        default=defaults.eta,       help='eta hyperparameter (default='+str(defaults.eta)+')')
+    parser.add_argument('--sigma',         type=float, dest='sigma',      default=defaults.sigma,     help='sigma hyperparameter (default='+str(defaults.sigma)+')')
+    parser.add_argument('--akernel',       type=str,   dest='akernel',     default=defaults.kernel,    help='local kernel type (G for Gaussian, L for Laplacian, myL for Laplacian for open-shell systems) (default '+defaults.kernel+')')
+    parser.add_argument('--gkernel',       type=str,   dest='gkernel',     default=defaults.gkernel,    help='global kernel type (avg for average kernel, rem for REMatch kernel) (default '+str(defaults.gkernel)+')')
+    parser.add_argument('--gdict',         nargs='*',   action=ParseKwargs, dest='gdict',     default=defaults.gdict,    help='dictionary like input string to initialize global kernel parameters')
+    parser.add_argument('--splits',        type=int,   dest='splits',     default=defaults.n_rep,     help='number of splits (default='+str(defaults.n_rep)+')')
+    parser.add_argument('--train',         type=float, dest='train_size', default=defaults.train_size, nargs='+', help='training set fractions')
+    parser.add_argument('--debug',         action='store_true', dest='debug', default=False,  help='enable debug')
+    parser.add_argument('--ll',            action='store_true', dest='ll',    default=False,  help='if correct for the numper of threads')
+    parser.add_argument('--readkernel',    action='store_true', dest='readk', default=False,  help='if X is kernel')
+    parser.add_argument('--sparse',        type=int, dest='sparse', default=None,  help='regression basis size for sparse learning')
+    parser.add_argument('--random_state',  type=int, dest='random_state', default=defaults.random_state,  help='random state for test / train splitting')
     args = parser.parse_args()
     print(vars(args))
     if(args.ll): correct_num_threads()
     X = np.load(args.repr)
     y = np.loadtxt(args.prop)
     maes_all = regression(X, y, read_kernel=args.readk, sigma=args.sigma, eta=args.eta, akernel=args.akernel,
-                          test_size=args.test_size, train_size=args.train_size, n_rep=args.splits, debug=args.debug)
+                          test_size=args.test_size, train_size=args.train_size, n_rep=args.splits, sparse=args.sparse,
+                          debug=args.debug)
     for size_train, meanerr, stderr in maes_all:
         print("%d\t%e\t%e" % (size_train, meanerr, stderr))
 
+
 if __name__ == "__main__":
     main()
-