Example for L-BFGS-B optimization of L1 regularization in kernel regression.

test_l1reg.py

# MIT License
#
# Copyright (c) 2017 Anders Steen Christensen
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

from __future__ import print_function

import os
import numpy as np
import qml
from qml.kernels import laplacian_kernel
from qml.math import cho_solve
from qml.math import bkf_solve

from qml.representations import get_slatm_mbtypes

from sklearn.linear_model import Lasso
from sklearn.linear_model import Lars
from sklearn.linear_model import Ridge
from  sklearn.linear_model import ElasticNet
from  sklearn.linear_model import LogisticRegression

from copy import deepcopy

import sys

__LAMBDA__ = float(sys.argv[1])

def get_energies(filename):
    """ Returns a dictionary with heats of formation for each xyz-file.
    """

    f = open(filename, "r")
    lines = f.readlines()
    f.close()

    energies = dict()

    for line in lines:
        tokens = line.split()

        xyz_name = tokens[0]
        hof = float(tokens[1])

        energies[xyz_name] = hof

    return energies

class L1(object):

    def __init__(self, K, Y):

        self.K = deepcopy(K)
        self.Y = deepcopy(Y)

        self.norm = -1.0
        self.iter = 0
        self.alpha = []

    def l1_norm(self, alpha):

        # self.norm = np.sum(np.abs(np.dot(self.K, alpha) - self.Y))

        Ka =  np.dot(self.K, alpha)

        T2 = __LAMBDA__/2.0 * np.matmul(alpha.T, Ka)

        self.norm = np.sum(np.abs(Ka - self.Y)) + T2
        # self.norm = np.linalg.norm(Ka - self.Y, ord=1) + T2

        return self.norm
    
    def l2_norm(self, alpha):

        self.norm = np.sqrt(np.sum(np.square(np.dot(self.K, alpha) - self.Y)))
        return self.norm

    def output(self, alpha):

        self.iter += 1
        self.alpha = deepcopy(alpha)
        print(" %7i  Norm = %20.10f   %20.10f" % (self.iter, self.norm, alpha[0]))

def test_krr_cmat():

    test_dir = os.path.dirname(os.path.realpath(__file__))

    # Parse file containing PBE0/def2-TZVP heats of formation and xyz filenames
    data = get_energies(test_dir + "/data/hof_qm7.txt")

    # Generate a list of qml.Compound() objects
    np.random.seed(666)

    xyz_files = sorted(data.keys())
    np.random.shuffle(xyz_files)
    
    mols = []

    for xyz_file in xyz_files[:1000]:

        # Initialize the qml.Compound() objects
        mol = qml.Compound(xyz=test_dir + "/qm7/" + xyz_file)

        # Associate a property (heat of formation) with the object
        mol.properties = data[xyz_file]

        # This is a Molecular Coulomb matrix sorted by row norm
        mol.generate_coulomb_matrix(size=23, sorting="row-norm")
        # mol.generate_bob()

        mols.append(mol)

    # Shuffle molecules

    # Make training and test sets
    n_train = int(len(mols) * 6400.0/7101)
    n_test  = len(mols) - n_train

    training = mols[:n_train]
    test  = mols[-n_test:]

    # List of representations
    X  = np.array([mol.representation for mol in training])
    Xs = np.array([mol.representation for mol in test])

    # List of properties
    Y = np.array([mol.properties for mol in training])
    Ys = np.array([mol.properties for mol in test])

    # Set hyper-parameters
    sigma = 1e4
    llambda = __LAMBDA__

    # Generate training Kernel
    K = laplacian_kernel(X, X, sigma)

    # Solve alpha
    C = deepcopy(K)
    C[np.diag_indices_from(C)] += llambda
    Ks = laplacian_kernel(Xs, X, sigma)
   

    cost = L1(K, Y)
    alpha1 = cho_solve(C,Y)
    # alpha1 = np.ones(len(Y))

    from scipy.optimize import minimize
    # from fista import Fista

    #K= np.ascontiguousarray(K)
    #Ks= np.ascontiguousarray(Ks)

    #  model = Fista(lambda_=0.5, loss="least-square", penalty="l11", recompute_Lipschitz_constant=True)

    # model.fit(K, Y, verbose=1)
    # Yss = model.predict(Ks)

    # np.set_printoptions(linewidth=100)
    # print(Yss - np.dot(Ks, alpha1))


    alpha = deepcopy(alpha1)
    print(alpha[0])
    res = minimize(cost.l1_norm, alpha, method="L-BFGS-B", callback = cost.output,
            options={"maxiter": 1000, "disp": True})

    print(alpha[0])
    alpha = deepcopy(res.x)
    print((alpha - alpha1)[:10])


    # alpha *= np.random.random(len(alpha))
    Yss = np.dot(Ks, alpha)
    Yss2 = np.dot(Ks, alpha1)

    # model = Lasso(alpha=1e-2, max_iter=1000, warm_start=True)
    # model = ElasticNet(l1_ratio=0.000001, warm_start=True)
    # model = Ridge(alpha=1e-7)
    # model = Ridge(alpha=1e-7)
    # model = LogisticRegression()
    # model = ElasticNet(alpha=1e-8, l1_ratio=0.8)
    # model.fit(K, Y)
    # Yss = model.predict(Ks)


    # Calculate prediction kernel

    #print(Yss)
    # print(Ys)

    mae = np.mean(np.abs(Ys - Yss))
    print(mae)

    mae = np.mean(np.abs(Ys - Yss2))
    print(mae)

    Y2 = np.dot(K, alpha)
    mae = np.mean(np.abs(Y - Y2))
    print(mae)
    print(np.sqrt(np.mean(np.square(Y - Y2))))

    Y2 = np.dot(K, alpha1)
    mae = np.mean(np.abs(Y - Y2))
    print(mae)

    print(np.sqrt(np.mean(np.square(Y - Y2))))

if __name__ == "__main__":
    test_krr_cmat()