Python tensorflow.matrix_triangular_solve() Examples

def _log_prob(self, given):
        mean, cov_tril = (self.path_param(self.mean),
                          self.path_param(self.cov_tril))
        log_det = 2 * tf.reduce_sum(
            tf.log(tf.matrix_diag_part(cov_tril)), axis=-1)
        n_dim = tf.cast(self._n_dim, self.dtype)
        log_z = - n_dim / 2 * tf.log(
            2 * tf.constant(np.pi, dtype=self.dtype)) - log_det / 2
        # log_z.shape == batch_shape
        if self._check_numerics:
            log_z = tf.check_numerics(log_z, "log[det(Cov)]")
        # (given-mean)' Sigma^{-1} (given-mean) =
        # (g-m)' L^{-T} L^{-1} (g-m) = |x|^2, where Lx = g-m =: y.
        y = tf.expand_dims(given - mean, -1)
        L, _ = maybe_explicit_broadcast(
            cov_tril, y, 'MultivariateNormalCholesky.cov_tril',
            'expand_dims(given, -1)')
        x = tf.matrix_triangular_solve(L, y, lower=True)
        x = tf.squeeze(x, -1)
        stoc_dist = -0.5 * tf.reduce_sum(tf.square(x), axis=-1)
        return log_z + stoc_dist 

def conditional_ND(self, Xnew, full_cov=False):
        ## modified from GPR
        Kx = self.kern.K(self._X_mean, Xnew)
        K = self.kern.K(self._X_mean) + tf.eye(tf.shape(self._X_mean)[0], dtype=settings.float_type) * self._lik_variance
        L = tf.cholesky(K)
        A = tf.matrix_triangular_solve(L, Kx, lower=True)
        V = tf.matrix_triangular_solve(L, self._Y - self.mean_function(self._X_mean))
        fmean = tf.matmul(A, V, transpose_a=True) + self.mean_function(Xnew)
        if full_cov:
            fvar = self.kern.K(Xnew) - tf.matmul(A, A, transpose_a=True)
            shape = tf.stack([1, 1, tf.shape(self._Y)[1]])
            fvar = tf.tile(tf.expand_dims(fvar, 2), shape)
        else:
            fvar = self.kern.Kdiag(Xnew) - tf.reduce_sum(tf.square(A), 0)
            fvar = tf.tile(tf.reshape(fvar, (-1, 1)), [1, tf.shape(self._Y)[1]])
        return fmean, fvar 

def _build_predict(self, n_test, full_cov=False):
    with tf.name_scope("build_predict"):
      self.x_test_pl = tf.identity(
          tf.placeholder(tf.float64, [n_test, self.input_dim], name="x_test_pl")
      )

    tf.logging.info("Using pre-computed Kernel")
    self.k_data_test = self.kern.k_full(self.x_pl, self.x_test_pl)

    with tf.name_scope("build_predict"):
      a = tf.matrix_triangular_solve(self.l, self.k_data_test)
      fmean = tf.matmul(a, self.v, transpose_a=True)

      if full_cov:
        fvar = self.kern.k_full(self.x_test_pl) - tf.matmul(
            a, a, transpose_a=True)
        shape = [1, 1, self.y_pl.shape[1]]
        fvar = tf.tile(tf.expand_dims(fvar, 2), shape)
      else:
        fvar = self.kern.k_diag(self.x_test_pl) - tf.reduce_sum(tf.square(a), 0)
        fvar = tf.tile(tf.reshape(fvar, (-1, 1)), [1, self.output_y.shape[1]])

      self.fmean = fmean
      self.fvar = fvar 

def build_likelihood_terms(self):
        Kdiag = self.kern.Kdiag(self.X)
        Kuu = make_Kuu(self.kern, self.a, self.b, self.ms)
        sigma2 = self.likelihood.variance

        # Compute intermediate matrices
        P = self.KufKfu / sigma2 + Kuu.get()
        L = tf.cholesky(P)
        log_det_P = tf.reduce_sum(tf.log(tf.square(tf.diag_part(L))))
        c = tf.matrix_triangular_solve(L, self.KufY) / sigma2

        # compute log marginal bound
        ND = tf.cast(tf.size(self.Y), float_type)
        D = tf.cast(tf.shape(self.Y)[1], float_type)
        return (-0.5 * ND * tf.log(2 * np.pi * sigma2),
                -0.5 * D * log_det_P,
                0.5 * D * Kuu.logdet(),
                -0.5 * self.tr_YTY / sigma2,
                0.5 * tf.reduce_sum(tf.square(c)),
                -0.5 * tf.reduce_sum(Kdiag)/sigma2,
                0.5 * Kuu.trace_KiX(self.KufKfu) / sigma2) 

def log_det_kron_sum(L1, L2):
    """
    L1 is a list of lower triangular arrays.
    L2 is a list of lower triangular arrays.

    if S1 = kron(L1) * kron(L1).T, and S2 similarly,
    this function computes the log determinant of S1 + S2
    """
    L1_logdets = [tf.reduce_sum(tf.log(tf.square(tf.diag_part(L)))) for L in L1]
    total_size = reduce(tf.mul, [L.shape[0] for L in L1])
    N_other = [total_size / tf.shape(L)[0] for L in L1]
    L1_logdet = reduce(tf.add, [s*ld for s, ld in zip(N_other, L1_logdets)])
    LiL = [tf.matrix_triangular_solve(L, R) for L, R in zip(L1, L2)]
    eigvals = [tf.self_adjoint_eigvals(tf.matmul(mat, mat.T)) for mat in LiL]
    eigvals_kronned = kron_vec_mul([tf.reshape(e, [-1, 1]) for e in eigvals], tf.ones([1, 1], tf.float64))
    return tf.reduce_sum(tf.log(1 + eigvals_kronned)) + L1_logdet 

def _define_full_covariance_probs(self, shard_id, shard):
    """Defines the full covariance probabilties per example in a class.

    Updates a matrix with dimension num_examples X num_classes.

    Args:
      shard_id: id of the current shard.
      shard: current data shard, 1 X num_examples X dimensions.
    """
    diff = shard - self._means
    cholesky = tf.cholesky(self._covs + self._min_var)
    log_det_covs = 2.0 * tf.reduce_sum(tf.log(tf.matrix_diag_part(cholesky)), 1)
    x_mu_cov = tf.square(
        tf.matrix_triangular_solve(
            cholesky, tf.transpose(
                diff, perm=[0, 2, 1]), lower=True))
    diag_m = tf.transpose(tf.reduce_sum(x_mu_cov, 1))
    self._probs[shard_id] = -0.5 * (
        diag_m + tf.to_float(self._dimensions) * tf.log(2 * np.pi) +
        log_det_covs) 

def __init__(self, prec_mean, prec, d=None):

        prec_mean = tf.convert_to_tensor(prec_mean)
        prec = tf.convert_to_tensor(prec)
        
        try:
            d1, = util.extract_shape(prec_mean)
            prec_mean = tf.reshape(prec_mean, (d1,1))
        except:
            d1,k = util.extract_shape(prec_mean)
            assert(k == 1)

            
        d2,_ = util.extract_shape(prec)
        assert(d1==d2)
        if d is None:
            d = d1
        else:
            assert(d==d1)

        super(MVGaussianNatural, self).__init__(d=d)

        self._prec_mean = prec_mean
        self._prec = prec
        
        self._L_prec = tf.cholesky(prec)
        self._entropy = util.dists.multivariate_gaussian_entropy(L_prec=self._L_prec)

        # want to solve prec * mean = prec_mean for mean.
        # this is equiv to (LL') * mean = prec_mean.
        # since tf doesn't have a cholSolve shortcut, just
        # do it directly:
        #   solve L y = prec_mean
        # to get y = (L' * mean), then
        #   solve L' mean = y
        y = tf.matrix_triangular_solve(self._L_prec, self._prec_mean, lower=True, adjoint=False)
        self._mean = tf.matrix_triangular_solve(self._L_prec, y, lower=True, adjoint=True)

        L_cov_transpose = util.triangular_inv(self._L_prec)
        self._L_cov = tf.transpose(L_cov_transpose)
        self._cov = tf.matmul(self._L_cov, L_cov_transpose) 

def _build_predict(self, Xnew, full_cov=False):
        Kuu = make_Kuu(self.kern, self.a, self.b, self.ms)
        sigma2 = self.likelihood.variance

        # Compute intermediate matrices
        P = self.KufKfu / sigma2 + Kuu.get()
        L = tf.cholesky(P)
        c = tf.matrix_triangular_solve(L, self.KufY) / sigma2

        Kus = make_Kuf(self.kern, Xnew, self.a, self.b, self.ms)
        tmp = tf.matrix_triangular_solve(L, Kus)
        mean = tf.matmul(tf.transpose(tmp), c)
        KiKus = Kuu.solve(Kus)
        if full_cov:
            var = self.kern.k(Xnew) + \
               tf.matmul(tf.transpose(tmp), tmp) - \
               tf.matmul(tf.transpose(KiKus), Kus)
            shape = tf.stack([1, 1, tf.shape(self.y)[1]])
            var = tf.tile(tf.expand_dims(var, 2), shape)
        else:
            var = self.kern.Kdiag(Xnew)
            var += tf.reduce_sum(tf.square(tmp), 0)
            var -= tf.reduce_sum(KiKus * Kus, 0)
            shape = tf.stack([1, tf.shape(self.Y)[1]])
            var = tf.tile(tf.expand_dims(var, 1), shape)
        return mean, var 

def _build_likelihood(self):
        num_data = tf.shape(self.Y)[0]
        output_dim = tf.shape(self.Y)[1]

        total_variance = reduce(tf.add, [k.variance for k in self.kerns])
        Kuu = [make_Kuu(k, ai, bi, self.ms) for k, ai, bi in zip(self.kerns, self.a, self.b)]
        Kuu = BlockDiagMat_many([mat for k in Kuu for mat in [k.A, k.B]])
        sigma2 = self.likelihood.variance

        # Compute intermediate matrices
        P = self.KufKfu / sigma2 + Kuu.get()
        L = tf.cholesky(P)
        log_det_P = tf.reduce_sum(tf.log(tf.square(tf.diag_part(L))))
        c = tf.matrix_triangular_solve(L, self.KufY) / sigma2

        # compute log marginal bound
        ND = tf.cast(num_data * output_dim, float_type)
        D = tf.cast(output_dim, float_type)
        bound = -0.5 * ND * tf.log(2 * np.pi * sigma2)
        bound += -0.5 * D * log_det_P
        bound += 0.5 * D * Kuu.logdet()
        bound += -0.5 * self.tr_YTY / sigma2
        bound += 0.5 * tf.reduce_sum(tf.square(c))
        bound += -0.5 * ND * total_variance / sigma2
        bound += 0.5 * D * Kuu.trace_KiX(self.KufKfu) / sigma2

        return bound 

def predict_components(self, Xnew):
        """
        Here, Xnew should be a Nnew x 1 array of points at which to test each function
        """
        Kuu = [make_Kuu(k, ai, bi, self.ms) for k, ai, bi in zip(self.kerns, self.a, self.b)]
        Kuu = BlockDiagMat_many([mat for k in Kuu for mat in [k.A, k.B]])
        sigma2 = self.likelihood.variance

        # Compute intermediate matrices
        P = self.KufKfu / sigma2 + Kuu.get()
        L = tf.cholesky(P)
        c = tf.matrix_triangular_solve(L, self.KufY) / sigma2

        Kus_blocks = [make_Kuf(k, Xnew, a, b, self.ms)
                      for i, (k, a, b) in enumerate(zip(self.kerns, self.a, self.b))]
        Kus = []
        start = tf.constant(0, tf.int32)
        for i, b in enumerate(Kus_blocks):
            zeros_above = tf.zeros(tf.stack([start, tf.shape(b)[1]]), float_type)
            zeros_below = tf.zeros(tf.stack([tf.shape(L)[0] - start - tf.shape(b)[0], tf.shape(b)[1]]), float_type)
            Kus.append(tf.concat([zeros_above, b, zeros_below], axis=0))
            start = start + tf.shape(b)[0]

        tmp = [tf.matrix_triangular_solve(L, Kus_i) for Kus_i in Kus]
        mean = [tf.matmul(tf.transpose(tmp_i), c) for tmp_i in tmp]
        KiKus = [Kuu.solve(Kus_i) for Kus_i in Kus]
        var = [k.Kdiag(Xnew[:, i:i+1]) for i, k in enumerate(self.kerns)]
        var = [v + tf.reduce_sum(tf.square(tmp_i), 0) for v, tmp_i in zip(var, tmp)]
        var = [v - tf.reduce_sum(KiKus_i * Kus_i, 0) for v, KiKus_i, Kus_i in zip(var, KiKus, Kus)]
        var = [tf.expand_dims(v, 1) for v in var]
        return tf.concat(mean, axis=1), tf.concat(var, axis=1) 

def _build_likelihood_terms(self):
        Kdiag = reduce(tf.multiply, [k.Kdiag(self.X[:, i:i+1]) for i, k in enumerate(self.kerns)])
        Kuu = [make_Kuu(k, a, b, self.ms) for k, a, b, in zip(self.kerns, self.a, self.b)]
        Kuu_solid = kron([Kuu_d.get() for Kuu_d in Kuu])
        Kuu_inv_solid = kron([Kuu_d.inv().get() for Kuu_d in Kuu])
        sigma2 = self.likelihood.variance

        # Compute intermediate matrices
        P = self.KufKfu / sigma2 + Kuu_solid
        L = tf.cholesky(P)
        log_det_P = tf.reduce_sum(tf.log(tf.square(tf.diag_part(L))))
        c = tf.matrix_triangular_solve(L, self.KufY) / sigma2

        Kuu_logdets = [K.logdet() for K in Kuu]
        N_others = [float(np.prod(self.Ms)) / M for M in self.Ms]
        Kuu_logdet = reduce(tf.add, [N*logdet for N, logdet in zip(N_others, Kuu_logdets)])

        # compute log marginal bound
        ND = tf.cast(tf.size(self.Y), float_type)
        D = tf.cast(tf.shape(self.Y)[1], float_type)
        return (-0.5 * ND * tf.log(2 * np.pi * sigma2),
                -0.5 * D * log_det_P,
                0.5 * D * Kuu_logdet,
                -0.5 * self.tr_YTY / sigma2,
                0.5 * tf.reduce_sum(tf.square(c)),
                -0.5 * tf.reduce_sum(Kdiag)/sigma2,
                0.5 * tf.reduce_sum(Kuu_inv_solid * self.KufKfu) / sigma2) 

def _build_predict(self, Xnew, full_cov=False):
        Kuu = [make_Kuu(k, a, b, self.ms) for k, a, b, in zip(self.kerns, self.a, self.b)]
        Kuu_solid = kron([Kuu_d.get() for Kuu_d in Kuu])
        Kuu_inv_solid = kron([Kuu_d.inv().get() for Kuu_d in Kuu])
        sigma2 = self.likelihood.variance

        # Compute intermediate matrices
        P = self.KufKfu / sigma2 + Kuu_solid
        L = tf.cholesky(P)
        c = tf.matrix_triangular_solve(L, self.KufY) / sigma2

        Kus = [make_Kuf(k, Xnew[:, i:i+1], a, b, self.ms)
               for i, (k, a, b) in enumerate(zip(self.kerns, self.a, self.b))]
        Kus = tf.transpose(make_kvs([tf.transpose(Kus_d) for Kus_d in Kus]))
        tmp = tf.matrix_triangular_solve(L, Kus)
        mean = tf.matmul(tf.transpose(tmp), c)
        KiKus = tf.matmul(Kuu_inv_solid, Kus)
        if full_cov:
            raise NotImplementedError
        else:
            var = reduce(tf.multiply, [k.Kdiag(Xnew[:, i:i+1]) for i, k in enumerate(self.kerns)])
            var += tf.reduce_sum(tf.square(tmp), 0)
            var -= tf.reduce_sum(KiKus * Kus, 0)
            shape = tf.stack([1, tf.shape(self.Y)[1]])
            var = tf.tile(tf.expand_dims(var, 1), shape)
        return mean, var 

def _build_likelihood(self):
        num_inducing = tf.size(self.ms)
        num_data = tf.shape(self.Y)[0]
        output_dim = tf.shape(self.Y)[1]

        err = self.Y - self.mean_function(self.X)
        Kdiag = self.kern.Kdiag(self.X)
        Kuf = make_Kuf(self.kern, self.X, self.a, self.b, self.ms)
        Kuu = make_Kuu(self.kern, self.a, self.b, self.ms)
        Kuu = Kuu.get()
        sigma = tf.sqrt(self.likelihood.variance)

        # Compute intermediate matrices
        L = tf.cholesky(Kuu)
        A = tf.matrix_triangular_solve(L, Kuf) / sigma
        AAT = tf.matmul(A, tf.transpose(A))

        B = AAT + tf.eye(num_inducing * 2 - 1, dtype=float_type)
        LB = tf.cholesky(B)
        log_det_B = 2. * tf.reduce_sum(tf.log(tf.diag_part(LB)))
        c = tf.matrix_triangular_solve(LB, tf.matmul(A, err)) / sigma

        # compute log marginal bound
        ND = tf.cast(num_data * output_dim, float_type)
        D = tf.cast(output_dim, float_type)
        bound = -0.5 * ND * tf.log(2 * np.pi * self.likelihood.variance)
        bound += -0.5 * D * log_det_B
        bound += -0.5 * tf.reduce_sum(tf.square(err))/self.likelihood.variance
        bound += 0.5 * tf.reduce_sum(tf.square(c))
        bound += -0.5 * tf.reduce_sum(Kdiag)/self.likelihood.variance
        bound += 0.5 * tf.reduce_sum(tf.diag_part(AAT))

        return bound 

def _build_predict(self, Xnew, full_cov=False):
        num_inducing = tf.size(self.ms)

        err = self.Y - self.mean_function(self.X)
        Kuf = make_Kuf(self.kern, self.X, self.a, self.b, self.ms)
        Kuu = make_Kuu(self.kern, self.a, self.b, self.ms)
        Kuu = Kuu.get()
        sigma = tf.sqrt(self.likelihood.variance)

        # Compute intermediate matrices
        L = tf.cholesky(Kuu)
        A = tf.matrix_triangular_solve(L, Kuf) / sigma
        AAT = tf.matmul(A, tf.transpose(A))

        B = AAT + tf.eye(num_inducing * 2 - 1, dtype=float_type)
        LB = tf.cholesky(B)
        c = tf.matrix_triangular_solve(LB, tf.matmul(A, err)) / sigma

        Kus = make_Kuf(self.kern, Xnew, self.a, self.b, self.ms)
        tmp1 = tf.matrix_triangular_solve(L, Kus, lower=True)
        tmp2 = tf.matrix_triangular_solve(LB, tmp1, lower=True)
        mean = tf.matmul(tf.transpose(tmp2), c)
        if full_cov:
            var = self.kern.K(Xnew) + \
                tf.matmul(tf.transpose(tmp2), tmp2) - \
                tf.matmul(tf.transpose(tmp1), tmp1)
            var = var[:, :, None] * tf.ones(self.Y.shape[1], dtype=float_type)
        else:
            var = self.kern.Kdiag(Xnew) + \
                tf.reduce_sum(tf.square(tmp2), 0) - \
                tf.reduce_sum(tf.square(tmp1), 0)
            var = var[:, None]#  * tf.ones(self.Y.shape[1], dtype=float_type)
        return mean + self.mean_function(Xnew), var 

def _elbo_helper(self):
        # compute kernel stuff
        kern = self.kern; f = self.q_mu;
        num_data = tf.shape(self.Z)[0]  # M
        num_func = tf.shape(f)[1]  # K
        Kmn = kern.Kzx(self.Z, self.X)
        Kmm = kern.Kzz(self.Z) + tf.eye(num_data, dtype=float_type) * settings.numerics.jitter_level
        Lm = tf.cholesky(Kmm)
        # Compute the projection matrix A
        A = tf.matrix_triangular_solve(Lm, Kmn, lower=True)
        # compute the covariance due to the conditioning
        fvar = kern.Kdiag_tr() - tf.reduce_sum(tf.square(A), 0)
        shape = tf.stack([num_func, 1])
        fvar = tf.tile(tf.expand_dims(fvar, 0), shape)  # K x N x N or K x N
        # another backsubstitution in the unwhitened case
        if not self.whiten:
            A = tf.matrix_triangular_solve(tf.transpose(Lm), A, lower=False)
        # construct the conditional mean
        fmean = tf.matmul(A, f, transpose_a=True)
        if self.q_sqrt is not None:
            if self.q_sqrt.get_shape().ndims == 2:
                LTA = A * tf.expand_dims(tf.transpose(self.q_sqrt), 2)  # K x M x N
            elif self.q_sqrt.get_shape().ndims == 3:
                L = tf.matrix_band_part(tf.transpose(self.q_sqrt, (2, 0, 1)), -1, 0)  # K x M x M
                A_tiled = tf.tile(tf.expand_dims(A, 0), tf.stack([num_func, 1, 1]))
                LTA = tf.matmul(L, A_tiled, transpose_a=True)  # K x M x N
            else:  # pragma: no cover
                raise ValueError("Bad dimension for q_sqrt: %s" %
                                 str(self.q_sqrt.get_shape().ndims))
            fvar = fvar + tf.reduce_sum(tf.square(LTA), 1)  # K x N
        fvar = tf.transpose(fvar)  # N x K or N x N x K
        return fmean + self.mean_function(self.X), fvar 

def _sqrt_solve(self, rhs):
    return tf.matrix_triangular_solve(self._chol, rhs, lower=True) 

def _build_inverse_projection(self, X, W, mu, std):
        n, d_latent = self.shape
        variance = tf.square(std)

        M = tf.matmul(tf.transpose(W), W) + tf.diag(tf.ones((d_latent))*variance)
        L = tf.cholesky(M)
        
        # pred_Z = M^-1 W' r' as per (6) in tipping & bishop JRSS
        # https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/bishop-ppca-jrss.pdf
        r = X-mu
        WTr = tf.transpose(tf.matmul(r, W)) # W' (t - mu)  
        tmp = tf.matrix_triangular_solve(L, WTr, lower=True)
        pred_z = tf.transpose(tf.matrix_triangular_solve(tf.transpose(L), tmp, lower=False))        
        
        return pred_z, L, std 

def _sample_and_entropy(self, **input_samples):
        # more efficient to re-use the same network for sampled value and entropy
        pred_z, L, std = self._build_inverse_projection(**input_samples)
        eps = tf.random_normal(shape=self.shape)
        unit_entropy = multivariate_gaussian_entropy(L_prec=L/std)
        
        self.canonical_pred = pred_z
        self.canonical_L = L        

        n, d_latent = self.shape
        sample = pred_z + tf.transpose(tf.matrix_triangular_solve(tf.transpose(L/std), tf.transpose(eps), lower=False))     
        entropy = n*unit_entropy        
        return sample, entropy 

def _sample(self, X, W, mu, std):
        pred_z, L, std = self._build_inverse_projection(X, W, mu, std)
        eps = tf.random_normal(shape=self.shape)
        
        # sample from N(pred_z, M^-1 * std**2)
        # where LL' = M
        return pred_z + tf.transpose(tf.matrix_triangular_solve(tf.transpose(L/std), tf.transpose(eps), lower=False)) 

def triangular_inv(L):
    eye = tf.diag(tf.ones_like(tf.diag_part(L)))
    invL = tf.matrix_triangular_solve(L, eye)
    return invL 

def cholesky_inverse(A):
  """Compute the inverse of `A` using Choselky decomposition. NOTE: `A` must be
  symmetric positive definite. This method of inversion is not completely stable since
  tf.cholesky is not always stable. Might raise `tf.errors.InvalidArgumentError`
  """
  N     = tf.shape(A)[0]
  L     = tf.cholesky(A)
  L_inv = tf.matrix_triangular_solve(L, tf.eye(N))
  A_inv = tf.matmul(L_inv, L_inv, transpose_a=True)
  return A_inv 

def _log_prob(self, given):
        mean, u_tril, v_tril = (self.path_param(self.mean),
                                self.path_param(self.u_tril),
                                self.path_param(self.v_tril))

        log_det_u = 2 * tf.reduce_sum(
            tf.log(tf.matrix_diag_part(u_tril)), axis=-1)
        log_det_v = 2 * tf.reduce_sum(
            tf.log(tf.matrix_diag_part(v_tril)), axis=-1)
        n_row = tf.cast(self._n_row, self.dtype)
        n_col = tf.cast(self._n_col, self.dtype)
        logZ = - (n_row * n_col) / 2. * \
            tf.log(2. * tf.constant(np.pi, dtype=self.dtype)) - \
            n_row / 2. * log_det_v - n_col / 2. * log_det_u
        # logZ.shape == batch_shape
        if self._check_numerics:
            logZ = tf.check_numerics(logZ, "log[det(Cov)]")

        y = given - mean
        y_with_last_dim_changed = tf.expand_dims(tf.ones(tf.shape(y)[:-1]), -1)
        Lu, _ = maybe_explicit_broadcast(
            u_tril, y_with_last_dim_changed,
            'MatrixVariateNormalCholesky.u_tril', 'expand_dims(given, -1)')
        y_with_sec_last_dim_changed = tf.expand_dims(tf.ones(
            tf.concat([tf.shape(y)[:-2], tf.shape(y)[-1:]], axis=0)), -1)
        Lv, _ = maybe_explicit_broadcast(
            v_tril, y_with_sec_last_dim_changed,
            'MatrixVariateNormalCholesky.v_tril',
            'expand_dims(given, -1)')
        x_Lb_inv_t = tf.matrix_triangular_solve(Lu, y, lower=True)
        x_t = tf.matrix_triangular_solve(Lv, tf.matrix_transpose(x_Lb_inv_t),
                                         lower=True)
        stoc_dist = -0.5 * tf.reduce_sum(tf.square(x_t), [-1, -2])
        return logZ + stoc_dist 

def _build_cholesky(self):
    tf.logging.info("Computing Kernel")
    self.k_data_data_reg = self.k_data_data + tf.eye(
        self.input_x.shape[0], dtype=tf.float64) * self.stability_eps
    if FLAGS.print_kernel:
      self.k_data_data_reg = tf.Print(
          self.k_data_data_reg, [self.k_data_data_reg],
          message="K_DD = ", summarize=100)
    self.l = tf.cholesky(self.k_data_data_reg)
    self.v = tf.matrix_triangular_solve(self.l, self.y_pl) 

def KL(self):
        """
        The KL divergence from the variational distribution to the prior

        :return: KL divergence from N(q_mu, q_sqrt) to N(0, I), independently for each GP
        """
        # if self.white:
        #     return gauss_kl(self.q_mu, self.q_sqrt)
        # else:
        #     return gauss_kl(self.q_mu, self.q_sqrt, self.Ku)

        self.build_cholesky_if_needed()

        KL = -0.5 * self.num_outputs * self.num_inducing
        KL -= 0.5 * tf.reduce_sum(tf.log(tf.matrix_diag_part(self.q_sqrt) ** 2))

        if not self.white:
            KL += tf.reduce_sum(tf.log(tf.matrix_diag_part(self.Lu))) * self.num_outputs
            KL += 0.5 * tf.reduce_sum(tf.square(tf.matrix_triangular_solve(self.Lu_tiled, self.q_sqrt, lower=True)))
            Kinv_m = tf.cholesky_solve(self.Lu, self.q_mu)
            KL += 0.5 * tf.reduce_sum(self.q_mu * Kinv_m)
        else:
            KL += 0.5 * tf.reduce_sum(tf.square(self.q_sqrt))
            KL += 0.5 * tf.reduce_sum(self.q_mu**2)

        return KL 

def gp_conditional(z, fz, x, full_cov, kernel, Kzz_chol=None):
    '''
    GP gp_conditional f(x) | f(z)==fz
    :param z: shape [n_z, n_covariates]
    :param fz: shape [n_particles, n_z]
    :param x: shape [n_x, n_covariates]
    :return: a distribution with shape [n_particles, n_x]
    '''
    n_z = int(z.shape[0])
    n_particles = tf.shape(fz)[0]

    if Kzz_chol is None:
        Kzz_chol = tf.cholesky(kernel(z, z))

    # Mean[fx|fz] = Kxz @ inv(Kzz) @ fz; Cov[fx|z] = Kxx - Kxz @ inv(Kzz) @ Kzx
    # With ill-conditioned Kzz, the inverse is often asymmetric, which
    # breaks further cholesky decomposition. We compute a symmetric one.
    Kzz_chol_inv = tf.matrix_triangular_solve(Kzz_chol, tf.eye(n_z))
    Kzz_inv = tf.matmul(tf.transpose(Kzz_chol_inv), Kzz_chol_inv)
    Kxz = kernel(x, z)  # [n_x, n_z]
    Kxziz = tf.matmul(Kxz, Kzz_inv)
    mean_fx_given_fz = tf.matmul(fz, tf.matrix_transpose(Kxziz))

    if full_cov:
        cov_fx_given_fz = kernel(x, x) - tf.matmul(Kxziz, tf.transpose(Kxz))
        cov_fx_given_fz = tf.tile(
            tf.expand_dims(tf.cholesky(cov_fx_given_fz), 0),
            [n_particles, 1, 1])
        fx_given_fz = zs.distributions.MultivariateNormalCholesky(
            mean_fx_given_fz, cov_fx_given_fz)
    else:
        # diag(AA^T) = sum(A**2, axis=-1)
        var = kernel.Kdiag(x) - \
            tf.reduce_sum(tf.matmul(
                Kxz, tf.matrix_transpose(Kzz_chol_inv)) ** 2, axis=-1)
        std = tf.sqrt(var)
        fx_given_fz = zs.distributions.Normal(
            mean=mean_fx_given_fz, std=std, group_ndims=1)
    return fx_given_fz 

def build_likelihood(self):

        # w = w./repmat(ell',[m,1]);                                              % scaled model angular frequencies
        w = self.omega / self.kern.lengthscales
        m = tf.shape(self.omega)[0]
        m_float = tf.cast(m, tf.float64)

        # phi = x_tr*w';
        phi = tf.matmul(self.X, tf.transpose(w))
        # phi = [cos(phi) sin(phi)];                                              % design matrix
        phi = tf.concat([tf.cos(phi), tf.sin(phi)], axis=1)

        # R = chol((sf2/m)*(phi'*phi) + sn2*eye(2*m));                            % calculate some often-used constants
        A = (self.kern.variance / m_float) * tf.matmul(tf.transpose(phi), phi)\
            + self.likelihood.variance * gpflow.tf_wraps.eye(2*m)
        RT = tf.cholesky(A)
        R = tf.transpose(RT)

        # PhiRiphi/R;
        # RtiPhit = PhiRi';
        RtiPhit = tf.matrix_triangular_solve(RT, tf.transpose(phi))
        # Rtiphity=RtiPhit*y_tr;
        Rtiphity = tf.matmul(RtiPhit, self.Y)

        # % output NLML
        # out1=0.5/sn2*(sum(y_tr.^2)-sf2/m*sum(Rtiphity.^2))+ ...
        out = 0.5/self.likelihood.variance*(tf.reduce_sum(tf.square(self.Y)) -
                                            self.kern.variance/m_float*tf.reduce_sum(tf.square(Rtiphity)))
        # +sum(log(diag(R)))+(n/2-m)*log(sn2)+n/2*log(2*pi);
        n = tf.cast(tf.shape(self.X)[0], tf.float64)
        out += tf.reduce_sum(tf.log(tf.diag_part(R)))\
            + (n/2.-m_float) * tf.log(self.likelihood.variance)\
            + n/2*np.log(2*np.pi)
        return -out 

def inv(self):
        di = tf.reciprocal(self.d)
        d_col = tf.expand_dims(self.d, 1)
        DiW = self.W / d_col
        M = tf.eye(tf.shape(self.W)[1], float_type) + tf.matmul(tf.transpose(DiW), self.W)
        L = tf.cholesky(M)
        v = tf.transpose(tf.matrix_triangular_solve(L, tf.transpose(DiW), lower=True))
        return LowRankMatNeg(di, V) 

def solve(self, B):
        d_col = tf.expand_dims(self.d, 1)
        DiB = B / d_col
        DiW = self.W / d_col
        WTDiB = tf.matmul(tf.transpose(DiW), B)
        M = tf.eye(tf.shape(self.W)[1], float_type) + tf.matmul(tf.transpose(DiW), self.W)
        L = tf.cholesky(M)
        tmp1 = tf.matrix_triangular_solve(L, WTDiB, lower=True)
        tmp2 = tf.matrix_triangular_solve(tf.transpose(L), tmp1, lower=False)
        return DiB - tf.matmul(DiW, tmp2) 

def _verifySolve(self, x, y, lower=True, adjoint=False, batch_dims=None, use_gpu=False):
    for np_type in [np.float32, np.float64]:
      a = x.astype(np_type)
      b = y.astype(np_type)
      # For numpy.solve we have to explicitly zero out the strictly
      # upper or lower triangle.
      if lower and a.size > 0:
        a_np = np.tril(a)
      elif a.size > 0:
        a_np = np.triu(a)
      else:
        a_np = a
      if adjoint:
        a_np = np.conj(np.transpose(a_np))

      if batch_dims is not None:
        a = np.tile(a, batch_dims + [1, 1])
        a_np = np.tile(a_np, batch_dims + [1, 1])
        b = np.tile(b, batch_dims + [1, 1])

      with self.test_session(use_gpu=use_gpu):
        tf_ans = tf.matrix_triangular_solve(a, b, lower=lower, adjoint=adjoint)
        out = tf_ans.eval()
        np_ans = np.linalg.solve(a_np, b)
        self.assertEqual(np_ans.shape, tf_ans.get_shape())
        self.assertEqual(np_ans.shape, out.shape)
        self.assertAllClose(np_ans, out) 

def test_MatrixTriangularSolve(self):
        t = tf.matrix_triangular_solve(*self.random((2, 3, 3, 3), (2, 3, 3, 1)), adjoint=False, lower=False)
        self.check(t)
        t = tf.matrix_triangular_solve(*self.random((2, 3, 3, 3), (2, 3, 3, 1)), adjoint=True, lower=False)
        self.check(t)
        t = tf.matrix_triangular_solve(*self.random((2, 3, 3, 3), (2, 3, 3, 1)), adjoint=False, lower=True)
        self.check(t)