sgcrf.py

# -*- coding: utf-8 -*-

from __future__ import division

from collections import namedtuple

import numpy as np
from numpy import random as rng
from scipy.linalg import cho_factor, cho_solve, LinAlgError
from sklearn.base import BaseEstimator
from progressbar import ProgressBar

from copy import deepcopy


def soft_thresh(r, w):
    return np.sign(w) * np.max(np.abs(w)-r, 0)


def check_pd(A, lower=True):
    """
    Checks if A is PD.
    If so returns True and Cholesky decomposition,
    otherwise returns False and None
    """
    try:
        return True, np.tril(cho_factor(A, lower=lower)[0])
    except LinAlgError as err:
        if 'not positive definite' in str(err):
            return False, None


def chol_inv(B, lower=True):
    """
    Returns the inverse of matrix A, where A = B*B.T,
    ie B is the Cholesky decomposition of A.

    Solves Ax = I
    given B is the cholesky factorization of A.
    """
    return cho_solve((B, lower), np.eye(B.shape[0]))


def inv(A):
    """
    Inversion of a SPD matrix using Cholesky decomposition.
    """
    return chol_inv(check_pd(A)[1])


def log(x):
    return np.log(x) if x > 0 else -np.inf


class SparseGaussianCRF(BaseEstimator):
    """
    GCRF models conditional probability density of y in R^q given x in R^p as
    p(y|x, Λ, Θ) = exp(-y' * Λ * y - 2 * x' * Θ * y) / Z(x)

    where Z(x) = c * |Λ|^-1 * exp(x' * Θ * Λ^-1 * Θ' * x)

    This is equivalent to:
    p(y|x) = N(-Θ * Λ^-1 * x, Λ^-1)

    Parameters
    ----------
    learning_rate : float, default 1.0
        step size

    lamL : float, default 0.01
        l1 regularization for the Λ matrix

    lamT : float, default 0.01
        l1 regularization for the Θ matrix

    References
    ----------

    Wytock and Kolter 2013
    Probabilistic Forecasting using Sparse Gaussian CRFs
    http://www.zicokolter.com/wp-content/uploads/2015/10/wytock-cdc13.pdf

    Wytock and Kolter 2013
    Sparse Gaussian CRFs Algorithms Theory and Application
    https://www.cs.cmu.edu/~mwytock/papers/gcrf_full.pdf

    McCarter and Kim 2015
    Large-Scale Optimization Algorithms for Sparse CGGMs
    http://arxiv.org/pdf/1509.04681.pdf

    McCarter and Kim 2016
    On Sparse Gaussian Chain Graph Models
    (info on Multi-Layer Sparse Gaussian CRFs)
    http://papers.nips.cc/paper/5320-on-sparse-gaussian-chain-graph-models.pdf

    Klinger and Tomanek 2007
    Classical Probabilistic Models and Conditional Random Fields
    http://www.scai.fraunhofer.de/fileadmin/images/bio/data_mining/paper/crf_klinger_tomanek.pdf

    Tong Tong Wu and Kenneth Lange 2008
    Coordinate Descent Algorithms for Lasso Penalized Regression
    http://arxiv.org/pdf/0803.3876.pdf

    """
    def __init__(self, learning_rate=1.0, lamL=1, lamT=1, n_iter=1000):
        self.lamL = lamL
        self.lamT = lamT
        self.learning_rate = learning_rate
        self.n_iter = n_iter

        self.Lam = None
        self.Theta = None

        # stuff for line search
        self.beta = 0.5
        self.slack = 0.05


    def fit(self, X, Y):
        """TODO: Docstring for fit.

        Parameters
        ----------
        X : np.array, shape (n_samples, input_dimension)
        Y : np.array, shape (n_samples, output_dimension)

        Returns
        -------
        TODO

        """
        assert X.shape[0] == Y.shape[0], 'Inputs and Outputs must have the same number of observations'

        self.alt_newton_coord_descent(X=X, Y=Y)
        return self


    def loss(self, X, Y, Lam=None, Theta=None):
        if Lam is None:
            Lam = self.Lam
        if Theta is None:
            Theta = self.Theta

        n, p, q = self._problem_size(X, Y)
        FixedParams = namedtuple('FixedParams', ['Sxx', 'Sxy', 'Syy'])
        VariableParams = namedtuple('VariableParams', ['Sigma', 'Psi'])

        fixed = FixedParams(Sxx=np.dot(X.T, X) / n,
                            Syy=np.dot(Y.T, Y) / n,
                            Sxy=np.dot(X.T, Y) / n)
        Sigma = inv(Lam)
        R = np.dot(np.dot(X, self.Theta), Sigma) / np.sqrt(n)
        vary = VariableParams(Sigma=Sigma,
                              Psi=np.dot(R.T, R))
        return self.l1_neg_log_likelihood(Lam, Theta, fixed, vary)


    # def check_gradient(self, fixed, vary):
    #     grad_lam = np.zeros_like(self.Lam)
    #     grad_theta = np.zeros_like(self.Theta)
    #     for i in range(self.Lam.shape[0]):
    #         for j in range(self.Lam.shape[1]):
    #             L = self.Lam.copy()
    #             run = 1e-10
    #             L[i,j] += run
    #             rise = self.neg_log_likelihood(L, self.Theta, fixed, vary) - self.neg_log_likelihood(self.Lam, self.Theta, fixed, vary)
    #             grad_lam[i,j] = rise/run
    #
    #     for i in range(self.Theta.shape[0]):
    #         for j in range(self.Theta.shape[1]):
    #             T = self.Theta.copy()
    #             run = 1e-10
    #             T[i,j] += run
    #             rise = self.neg_log_likelihood(self.Lam, T, fixed, vary) - self.neg_log_likelihood(self.Lam, self.Theta, fixed, vary)
    #             grad_theta[i,j] = rise/run
    #     return grad_lam, grad_theta


    def neg_log_likelihood(self, Lam, Theta, fixed, vary):
        "compute the negative log-likelihood of the GCRF"
        return -log(np.linalg.det(Lam)) + \
                np.trace(np.dot(fixed.Syy, Lam) + \
                2*np.dot(fixed.Sxy.T, Theta) + \
                np.dot(vary.Psi, Lam))


    @staticmethod
    def l1_norm_off_diag(A):
        "convenience method for l1 norm, excluding diagonal"
        # let's speed this up later
        # assume A symmetric, sparse too
        return np.linalg.norm(A - np.diag(A.diagonal()), ord=1)


    def l1_neg_log_likelihood(self, Lam, Theta, fixed, vary):
        "regluarized negative log likelihood"
        return self.neg_log_likelihood( Lam, Theta, fixed, vary) + \
               self.lamL * self.l1_norm_off_diag(Lam) + \
               self.lamT * np.linalg.norm(Theta, ord=1)


    def neg_log_likelihood_wrt_Lam(self, Lam, fixed, vary):
        # compute the negative log-likelihood of the GCRF when Theta is fixed
        return -log(np.linalg.det(Lam)) + \
                np.trace(np.dot(fixed.Syy, Lam) + \
                np.dot(vary.Psi, Lam))


    def l1_neg_log_likelihood_wrt_Lam(self, Lam, fixed, vary):
        # regularized neg log loss
        return self.neg_log_likelihood_wrt_Lam(Lam, fixed, vary) + \
               self.lamL * np.linalg.norm(Lam - np.diag(Lam.diagonal()), ord=1)

    # def neg_log_likelihood_wrt_Theta(self, Theta, fixed, vary):
    #     # compute the negative log-likelihood of the GCRF when Lamba is fixed
    #     return 2*np.dot(fixed.Sxy.T, Theta) + np.dot(vary.Sigma, vary.Psi)

    def grad_wrt_Lam(self, fixed, vary):
        return fixed.Syy - vary.Sigma - vary.Psi


    def grad_wrt_Theta(self, fixed, vary):
        # TODO this is not avoiding the Gamma computation!!!
        # gamma = Sxx Theta Sigma
        return 2 * fixed.Sxy + 2 * np.dot(fixed.Sxx, np.dot(self.Theta, vary.Sigma))


    def active_set(self, fixed, vary):
        return (self.active_set_Lam(fixed, vary),
                self.active_set_Theta(fixed, vary))


    def active_set_Lam(self, fixed, vary):
        grad = self.grad_wrt_Lam(fixed, vary)
        assert np.allclose(grad, grad.T, 1e-3)
        return np.where((np.abs(np.triu(grad)) > self.lamL) | (self.Lam != 0))
        # return np.where((np.abs(grad) > self.lamL) | (~np.isclose(self.Lam, 0)))


    def active_set_Theta(self, fixed, vary):
        grad = self.grad_wrt_Theta(fixed, vary)
        return np.where((np.abs(grad) > self.lamT) | (self.Theta != 0))
        # return np.where((np.abs(grad) > self.lamT) | (~np.isclose(self.Theta, 0)))


    def _problem_size(self, X, Y):
        (n, p), q = X.shape, Y.shape[1]
        return n, p, q


    def check_descent(self, newton_lambda, alpha, fixed, vary):
        # check if we have made suffcient descent
        DLam = np.trace(np.dot(self.grad_wrt_Lam(fixed, vary), newton_lambda)) + \
               self.lamL * np.linalg.norm(self.Lam + newton_lambda, ord=1) - \
               self.lamL * np.linalg.norm(self.Lam, ord=1)

        nll_a = self.l1_neg_log_likelihood_wrt_Lam(self.Lam + alpha * newton_lambda, fixed, vary)
        nll_b = self.l1_neg_log_likelihood_wrt_Lam(self.Lam, fixed, vary) + alpha * self.slack * DLam
        return nll_a <= nll_b


    def check_descent2(self, newton_lambda, alpha, fixed, vary):

        lhs = self.l1_neg_log_likelihood(self.Lam + alpha*newton_lambda, self.Theta, fixed, vary)

        mu = np.trace(np.dot(self.grad_wrt_Lam(fixed, vary), newton_lambda)) + \
             self.lamL*self.l1_norm_off_diag(self.Lam + newton_lambda) +\
             self.lamT*np.linalg.norm(self.Theta, ord=1)

        rhs = self.neg_log_likelihood(self.Lam, self.Theta, fixed, vary) +\
              alpha * self.slack * mu
        return lhs <= rhs

    def line_search(self, newton_lambda, fixed, vary):
        # returns cholesky decomposition of Lambda and the learning rate
        alpha = self.learning_rate
        while True:
            pd, L = check_pd(self.Lam + alpha * newton_lambda)
            if pd and self.check_descent(newton_lambda, alpha, fixed, vary):
                # step is positive definite and we have sufficient descent
                break
                #TODO maybe want to return newt+alpha, to reuse computation
            alpha = alpha * self.beta
            # if alpha < 0.1:
            #   return L, alpha
        return L, alpha


    def lambda_newton_direction(self, active, fixed, vary, max_iter=1):
        # TODO we should be able to do a warm start...
        delta = np.zeros_like(vary.Sigma)
        U = np.zeros_like(vary.Sigma)

        for _ in range(max_iter):
            for i, j in rng.permutation(np.array(active).T):
                if i > j:
                    # seems ok since we look for upper triangular indices in active set
                    continue

                if i==j:
                    a = vary.Sigma[i,i] ** 2 + 2 * vary.Sigma[i,i] * vary.Psi[i,i]
                else:
                    a = (vary.Sigma[i, j] ** 2 + vary.Sigma[i, i] * vary.Sigma[j, j] +
                         vary.Sigma[i, i] * vary.Psi[j, j] + 2 * vary.Sigma[i, j] * vary.Psi[i, j] +
                         vary.Sigma[j, j] * vary.Psi[i, i])

                b = (fixed.Syy[i, j] - vary.Sigma[i, j] - vary.Psi[i, j] +
                     np.dot(vary.Sigma[i,:], U[:,j]) +
                     np.dot(vary.Psi[i,:], U[:,j]) +
                     np.dot(vary.Psi[j,:], U[:,i]))

                if i==j:
                    u = -b/a
                    delta[i, i] += u
                    U[i, :] +=  u * vary.Sigma[i, :]
                else:
                    c = self.Lam[i, j] + delta[i, j]
                    u = soft_thresh(self.lamL / a, c - b/a) - c
                    delta[j, i] += u
                    delta[i, j] += u
                    U[j, :] +=  u * vary.Sigma[i, :]
                    U[i, :] +=  u * vary.Sigma[j, :]

        return delta


    def theta_coordinate_descent(self, active, fixed, vary, max_iter=1):
        V = np.dot(self.Theta, vary.Sigma)
        for _ in range(max_iter):
            for i, j in np.array(active).T:
                a = 2 * vary.Sigma[j, j] * fixed.Sxx[i, i]
                b = 2 * fixed.Sxy[i, j] + 2 * np.dot(fixed.Sxx[i,:], V[:,j])
                c = self.Theta[i, j]

                u = soft_thresh(self.lamT / a, c - b/a) - c

                self.Theta[i, j] += u
                V[i,:] += u * vary.Sigma[j,:]

        return self.Theta


    def alt_newton_coord_descent(self, X, Y):
        """
        JK Trying to follow Calvin's algorithm, then will merge back into orig
        """
        n, p, q = self._problem_size(X, Y)

        FixedParams = namedtuple('FixedParams', ['Sxx', 'Sxy', 'Syy'])
        VariableParams = namedtuple('VariableParams', ['Sigma', 'Psi'])

        fixed = FixedParams(Sxx=np.dot(X.T, X) / n,
                            Syy=np.dot(Y.T, Y) / n,
                            Sxy=np.dot(X.T, Y) / n)

        # allow for continued fitting
        if self.Lam is None:
            self.Lam = np.eye(q)
            Sigma = np.eye(q)
        else:
            Sigma = inv(self.Lam) # use cholesky decomp then solve system of eqs
        if self.Theta is None:
            self.Theta = np.zeros((p, q))

        self.nll = []
        self.lnll = []
        self.lrs = []
        from progressbar import ProgressBar
        pbar = ProgressBar()
        for it in pbar(range(self.n_iter)):

            # update variable params
            R = np.dot(np.dot(X, self.Theta), Sigma) / np.sqrt(n)
            vary = VariableParams(Sigma=Sigma,
                                  Psi=np.dot(R.T, R))

            self.nll.append(self.neg_log_likelihood(self.Lam, self.Theta, fixed, vary))
            self.lnll.append(self.l1_neg_log_likelihood(self.Lam, self.Theta, fixed, vary))

            # determine active set
            active_Lam = self.active_set_Lam(fixed, vary)

            # solve D_lambda via coordinate descent
            newton_lambda = self.lambda_newton_direction(active_Lam, fixed, vary, max_iter=1)

            # line search for best step size
            learning_rate = self.learning_rate
            LL, learning_rate = self.line_search(newton_lambda, fixed, vary)
            self.lrs.append(learning_rate)
            self.Lam = self.Lam.copy() + learning_rate * newton_lambda

            # update variable params
            Sigma = chol_inv(LL) # use chol decomp from the backtracking
            vary = VariableParams(Sigma=Sigma,
                                  Psi=None) # dont need psi here

            # determine active set
            active_Theta = self.active_set_Theta(fixed, vary)

            # solve theta
            self.Theta = self.theta_coordinate_descent(active_Theta, fixed, vary, max_iter=1)


    def sample(self, X, n=1, verbose=True):
        """
        Draw samples from the conditional probability for y given x.

        Inference in  GCRF given by:
        y|x ~ N(-Θ * Λ^-1 * x, Λ^-1)

        This algorithm uses some clever accelerations to generate samples.
        Inspired by Calvin McCarter[1].

        The idea is to make the desired number of draws from a white multivariate
        normal distribution, and then multiply these draws by the cholesky
        decomposition of the desired covariance matrix. This adds the necessary
        color to the original draws.

        An equivalent thing to do is to multiply the white draws by the inverse of the
        cholesky decompostition of the desired precision matrix:

        If S is the desired covariance matrix and L, the precision matrix then:
        S = SL * SL.T # cholesky decomposition
        L = LL * LL.T # cholesky decomposition

        S = L^-1 # by asssumption, then
        SL * SL.T = (LL * LL.T)^-1 = LL.T ^-1 * LL ^-1

        Although SL != LL.T^-1, the effective coloring is the same.

        A final acceleration is to solve linear systems of equations instead of
        explicitly computing matrix inversions.

        [1] https://calvinmccarter.wordpress.com/2015/01/06/multivariate-normal-random-number-generation-in-matlab/
        """
        LL = np.linalg.cholesky(self.Lam)
        Sigma = chol_inv(LL)

        means = -np.dot(np.dot(Sigma, self.Theta.T), X.T)
        means = np.tile(np.atleast_2d(means), n)
        N = means.shape[1]

        z = rng.randn(self.Lam.shape[0], N)
        samples = np.linalg.solve(LL.T, z) + means

        return samples.squeeze().T


    def predict(self, X, Y=None):
        """
        Return the mean of y given x.

        Inference in  GCRF given by:
        y|x ~ N(-Θ * Λ^-1 * x, Λ^-1)

        so this method returns:
        -Θ * Λ^-1 * x
        """
        return -np.dot(np.dot(inv(self.Lam), self.Theta.T), X.T).T


    def get_params(self, deep=True):
        return {'lamL': self.lamL,
                'lamT': self.lamT,
                'learning_rate': self.learning_rate,
                'n_iter': self.n_iter}


    def set_params(self, **parameters):
        for parameter, value in parameters.items():
            self.setattr(parameter, value)