python类cholesky()的实例源码

def test_sample_independent_momentum(self):
        mom_resample_coeff = 1.
        dtype = np.float64
        for n_dim in [10, 100, 1000]:
            mass_matrix = self.prng.normal(size=(n_dim, n_dim))
            mass_matrix = mass_matrix.dot(mass_matrix.T)
            mass_matrix_chol = la.cholesky(mass_matrix, lower=True)
            sampler = uhmc.EuclideanMetricHmcSampler(
                energy_func=energy_func,
                mass_matrix=mass_matrix,
                energy_grad=energy_grad,
                prng=self.prng,
                mom_resample_coeff=mom_resample_coeff,
                dtype=dtype)
            pos = self.prng.normal(size=(n_dim,)).astype(dtype)
            mom = sampler.sample_independent_momentum_given_position(
                pos, cache={}
            )
            assert mom.ndim == pos.ndim and mom.shape[0] == pos.shape[0], (
                'Momentum sampling returning incorrect shaped array.'
            )
            assert mom.dtype == pos.dtype, (
                'Momentum sampling returning array with incorrect dtype.'
            )
            sum_std = sum(mass_matrix.diagonal()**0.5)
            assert abs(mom.mean()) < 5. * sum_std / n_dim**0.5, (
                'Mean of sampled momentum > 5 std. from expected value.'
            )

def __init__(self, energy_func, mass_matrix, energy_grad=None, prng=None,
                 mom_resample_coeff=1., dtype=np.float64):
        super(EuclideanMetricHmcSampler, self).__init__(
            energy_func, energy_grad, prng, mom_resample_coeff, dtype)
        self.mass_matrix = mass_matrix
        self.mass_matrix_chol = la.cholesky(mass_matrix, lower=True)

def matrixInverse(M):
    return chol2inv(spla.cholesky(M, lower=False))

def test_initialize(self,window):
        points = np.random.randn(self.n_points,self.dim)
        GM = GaussianMixture(self.n_components,window=window)
        GM.initialize(points)

        checking.verify_covariance(GM.get('cov'),self.n_components,self.dim)
        checking.verify_means(GM.get('means'),self.n_components,self.dim)
        checking.verify_log_pi(GM.get('log_weights'),self.n_components)

        cov_chol = np.empty_like(GM.get('cov'))
        for i in range(self.n_components):
            cov_chol[i] = linalg.cholesky(GM.get('cov')[i],lower=True)

        assert_almost_equal(cov_chol,GM.get('cov_chol'))
        assert GM.get('_is_initialized') == True

def test_update(self,window):
        points = np.random.randn(self.n_points,self.dim)
        GM = GaussianMixture(self.n_components,window=window,update=True)

        GM.initialize(points)
        GM.fit(points)

        expected_cov_chol = np.zeros((self.n_components,self.dim,self.dim))
        for i in range(self.n_components):
            expected_cov_chol[i] = linalg.cholesky(GM.get('cov')[i],lower=True)

        predected_cov_chol = GM.get('cov_chol')

        assert_almost_equal(expected_cov_chol,predected_cov_chol)

def test_compute_precisions_chol_full():
    n_components, n_features = 5,2

    cov = generate_covariance_matrices_full(n_components,n_features)
    expected_precisions_chol = np.empty((n_components,n_features,n_features))
    for i in range(n_components):
        cov_chol = linalg.cholesky(cov[i],lower=True)
        expected_precisions_chol[i] = np.linalg.inv(cov_chol).T

    predected_precisions_chol = _compute_precisions_chol(cov,'full')

    assert_almost_equal(expected_precisions_chol,predected_precisions_chol)

def factor(X, rho):
    """
    computes cholesky factorization of the kernel K = 1/rho*XX^T + I
    Input:
    X design matrix: n_s x n_f (we assume n_s << n_f)
    rho: regularizaer
    Output:
    L  lower triangular matrix
    U  upper triangular matrix
    """
    n_s, n_f = X.shape
    K = 1 / rho * scipy.dot(X, X.T) + scipy.eye(n_s)
    U = linalg.cholesky(K)
    return U

def _pull_from_cache_or_compute(self):
        if self.caching and len(self._cache_list) == self.num_states:
            chol  = self._cache_list[self.state]['chol']
            alpha = self._cache_list[self.state]['alpha']
        else:
            chol  = spla.cholesky(self.kernel.cov(self.inputs), lower=True)
            alpha = spla.cho_solve((chol, True), self.values - self.mean.value)

        return chol, alpha

def _prepare_cache(self):
        inputs_hash = hash(self.inputs.tostring())
        for i in xrange(self.num_states):
            self.set_state(i)
            chol  = spla.cholesky(self.kernel.cov(self.inputs), lower=True)
            alpha = spla.cho_solve((chol, True), self.values - self.mean.value)
            cache_dict = {
                'chol'  : chol,
                'alpha' : alpha
            }
            self._cache_list.append(cache_dict)

def log_likelihood(self):
        """
        GP Marginal likelihood

        Notes
        -----
        This is called by the samplers when fitting the hyperparameters.
        """
        cov   = self.kernel.cov(self.observed_inputs)
        chol  = spla.cholesky(cov, lower=True)
        solve = spla.cho_solve((chol, True), self.observed_values - self.mean.value)

        # Uses the identity that log det A = log prod diag chol A = sum log diag chol A
        return -np.sum(np.log(np.diag(chol)))-0.5*np.dot(self.observed_values - self.mean.value, solve)

def _compute_implied_y(self, model, nu):
        K_XX = model.noiseless_kernel.cov(model.inputs)
        L    = spla.cholesky(K_XX, lower=True) 

        return np.dot(L, nu) + model.mean.value

def sample_fun(self, model, **sampler_options):
        params_array = hyperparameter_utils.params_to_array(self.params)

        if model.has_data:
            K_XX      = model.noiseless_kernel.cov(model.inputs)
            current_L = spla.cholesky(K_XX, lower=True)
            nu        = spla.solve_triangular(current_L, model.latent_values.value-model.mean.value, lower=True)
        else:
            nu = None # if no data

        new_params, current_ll = slice_sample(params_array, self.logprob, model, nu, **sampler_options)

        new_latent_values = self._compute_implied_y(model, nu)

        return new_params, new_latent_values, current_ll

def chol_add(L, A):
    N = L.shape[0]
    S12 = spla.solve_triangular(L, A[:N,N:], lower=True)
    S22 = spla.cholesky(A[N:,N:] - S12.T.dot(S12)).T
    L_update = np.zeros(A.shape)
    L_update[:N,:N] = L
    L_update[N:,:N] = S12.T
    L_update[N:,N:] = S22
    return L_update

def _pull_from_cache_or_compute(self):
        if self.caching and len(self._cache_list) == self.num_states:
            chol  = self._cache_list[self.state]['chol']
            alpha = self._cache_list[self.state]['alpha']
        else:
            chol  = spla.cholesky(self.kernel.cov(self.inputs), lower=True)
            alpha = spla.cho_solve((chol, True), self.values - self.mean.value)

        return chol, alpha

def _prepare_cache(self):
        inputs_hash = hash(self.inputs.tostring())
        for i in xrange(self.num_states):
            self.set_state(i)
            chol  = spla.cholesky(self.kernel.cov(self.inputs), lower=True)
            alpha = spla.cho_solve((chol, True), self.values - self.mean.value)
            cache_dict = {
                'chol'  : chol,
                'alpha' : alpha
            }
            self._cache_list.append(cache_dict)

def log_likelihood(self):
        """
        GP Marginal likelihood

        Notes
        -----
        This is called by the samplers when fitting the hyperparameters.
        """
        cov   = self.kernel.cov(self.observed_inputs)
        chol  = spla.cholesky(cov, lower=True)
        solve = spla.cho_solve((chol, True), self.observed_values - self.mean.value)

        # Uses the identity that log det A = log prod diag chol A = sum log diag chol A
        return -np.sum(np.log(np.diag(chol)))-0.5*np.dot(self.observed_values - self.mean.value, solve)

def _compute_implied_y(self, model, nu):
        K_XX = model.noiseless_kernel.cov(model.inputs)
        L    = spla.cholesky(K_XX, lower=True) 

        return np.dot(L, nu) + model.mean.value

def sample(self, model):
        if not model.has_data:
            return np.zeros(0) # TODO this should be a sample from the prior...

        prior_cov      = model.noiseless_kernel.cov(model.inputs)
        prior_cov_chol = spla.cholesky(prior_cov, lower=True)
        # Here get the Cholesky from model

        params_array = hyperparameter_utils.params_to_array(self.params)
        for i in xrange(self.thinning + 1):
            params_array, current_ll = elliptical_slice(
                params_array,
                self.logprob,
                prior_cov_chol,
                model.mean.value,
                model,
                **self.sampler_options
            )
            hyperparameter_utils.set_params_from_array(self.params, params_array)
        self.current_ll = current_ll # for diagnostics


# xx: the initial point
# sample_nu: a function that samples from the multivariate Gaussian prior
# log_like_fn: a function that computes the log likelihood of an input
# cur_log_like (optional): the current log likelihood
# angle_range: not sure

def chol_add(L, A):
    N = L.shape[0]
    S12 = spla.solve_triangular(L, A[:N,N:], lower=True)
    S22 = spla.cholesky(A[N:,N:] - S12.T.dot(S12)).T
    L_update = np.zeros(A.shape)
    L_update[:N,:N] = L
    L_update[N:,:N] = S12.T
    L_update[N:,N:] = S22
    return L_update

def fit(self, X, y):
        """Fit the parameters of the probabilistic process based on the
        available training data.
        """
        # Store the training data (both the inputs and the targets).
        self.X, self.y = X, y
        # Compute the covariance matrix of the observed inputs.
        K = self.kernel.cov(self.X)
        # For a numerically stable algorithm, we use Cholesky decomposition.
        self.L = spla.cholesky(K, lower=True)
        self.alpha = spla.cho_solve((self.L, True), self.y).ravel()
        L_inv = spla.solve_triangular(self.L.T, np.eye(self.L.shape[0]))
        self.K_inv = L_inv.dot(L_inv.T)
        self.beta = self.y.dot(self.alpha)