Python linalg.cholesky方法代码示例

# 需要导入模块: from numpy import linalg [as 别名]
# 或者: from numpy.linalg import cholesky [as 别名]
def _t(u, rho, nu):
    d = u.shape[1]
    nu = float(nu)
    
    try:
        R = cholesky(rho)
    except LinAlgError:
        raise ValueError('Provided Rho matrix is not Positive Definite!')
    
    ticdf = t.ppf(u, nu)
    
    z = solve(R,ticdf.T)
    z = z.T
    logSqrtDetRho = np.sum(np.log(np.diag(R)))
    const = gammaln((nu+d)/2.0) + (d-1)*gammaln(nu/2.0) - d*gammaln((nu+1)/2.0) - logSqrtDetRho
    sq = np.power(z,2)
    summer = np.sum(np.power(z,2),axis=1)
    numer = -((nu+d)/2.0) * np.log(1.0 + np.sum(np.power(z,2),axis=1)/nu)
    denom = np.sum(-((nu+1)/2) * np.log(1 + (np.power(ticdf,2))/nu), axis=1)
    y = np.exp(const + numer - denom)
    
    return y 

# 需要导入模块: from numpy import linalg [as 别名]
# 或者: from numpy.linalg import cholesky [as 别名]
def test_basic_property(self):
        # Check A = L L^H
        shapes = [(1, 1), (2, 2), (3, 3), (50, 50), (3, 10, 10)]
        dtypes = (np.float32, np.float64, np.complex64, np.complex128)

        for shape, dtype in itertools.product(shapes, dtypes):
            np.random.seed(1)
            a = np.random.randn(*shape)
            if np.issubdtype(dtype, np.complexfloating):
                a = a + 1j*np.random.randn(*shape)

            t = list(range(len(shape)))
            t[-2:] = -1, -2

            a = np.matmul(a.transpose(t).conj(), a)
            a = np.asarray(a, dtype=dtype)

            c = np.linalg.cholesky(a)

            b = np.matmul(c, c.transpose(t).conj())
            assert_allclose(b, a,
                            err_msg="{} {}\n{}\n{}".format(shape, dtype, a, c),
                            atol=500 * a.shape[0] * np.finfo(dtype).eps) 

# 需要导入模块: from numpy import linalg [as 别名]
# 或者: from numpy.linalg import cholesky [as 别名]
def test_design_r(self):
        design = simple_mv_velocity_design(3)
        batch_design = design.for_batch(2, 1)

        cov = batch_design.R(0)[0]
        self.assertTupleEqual(cov.size(), (3, 3))

        self.assertTrue(cov.requires_grad)
        cholesky_log_diag = design.measure_covariance.param_dict()['cholesky_log_diag']
        cholesky_off_diag = design.measure_covariance.param_dict()['cholesky_off_diag']

        cov = cov.data.numpy()
        self.assertTrue(np.isclose(cov, cov.T).all(), msg="Covariance is not symmetric.")
        chol = cholesky(cov)

        for a, b in zip(torch.exp(cholesky_log_diag).tolist(), np.diag(chol).tolist()):
            self.assertAlmostEqual(a, b, places=4)

        for a, b in zip(cholesky_off_diag.tolist(), chol[np.tril_indices_from(chol, k=-1)].tolist()):
            self.assertAlmostEqual(a, b, places=4) 

# 需要导入模块: from numpy import linalg [as 别名]
# 或者: from numpy.linalg import cholesky [as 别名]
def dgp_dL_via_Sigma(self, L: np.ndarray, L_inv: np.ndarray, dsigma: np.ndarray) -> np.ndarray:
        """
        Partial derivatives of the gp posterior samples with respect to the cholesky of the posterior covariance matrix given the partial derivative values with respect to the posterior covariance matrix.
        
        :param s: Samples of observations from the posterior distribution of the model
        :param L: Cholesky decomposition(s) of the posterior covariance matrix (samples)
        :param L_inv: Inverse(s) of Cholesky decomposition(s) of the posterior covariance matrix (samples)
        :param dsigma: Partial derivatives with respect to the posterior covariance matrix
        :return: the derivative of the gp samples with respect to the choleskies
        """
        E = np.tril(2*np.ones((L.shape[1],L.shape[2])),-1 ) +np.eye(L.shape[2])
        dl = np.empty(dsigma.shape) #np.empty((N,b,b,b,d))
        for i in range(dsigma.shape[3]):
            for j in range(dsigma.shape[4]):
                tmp1 = np.matmul(L_inv, dsigma[:,:,:,i,j]) # N x b x b
                tmp2 = np.matmul(tmp1, np.swapaxes(L_inv,1,2)) # N x b x b
                tmp3 = tmp2 * E[None,:,:] # N x b x b
                dl[:,:,:,i,j] = 0.5 * np.matmul(L, tmp3) # N x b x b
        return dl # N x b x b 

# 需要导入模块: from numpy import linalg [as 别名]
# 或者: from numpy.linalg import cholesky [as 别名]
def update_item_params(self):
        x_bar = np.mean(self.item_features, 0).T
        x_bar = np.reshape(x_bar, (self.num_features, 1))
        S_bar = np.cov(self.item_features.T)
        norm_X_bar = self.mu0_item - x_bar

        WI_post = inv(inv(self.WI_item) + self.num_items * S_bar + \
            np.dot(norm_X_bar, norm_X_bar.T) * \
            (self.num_items * self.beta_item) / (self.beta_item + self.num_items))

        # Not sure why we need this...
        WI_post = (WI_post + WI_post.T) / 2.0

        # update alpha_item
        self.alpha_item = sample_wishart(WI_post, self.df_post_item)

        # update mu_item
        mu_temp = (self.beta_item * self.mu0_item + self.num_items * x_bar) / \
            (self.beta_item + self.num_items)
        lam = cholesky(inv((self.beta_item + self.num_items) * self.alpha_item))
        # lam = lam.T
        self.mu_item = mu_temp + np.dot(lam, NormalRandom.generate_matrix(self.num_features, 1))

        # raise ValueError('AAA') 

# 需要导入模块: from numpy import linalg [as 别名]
# 或者: from numpy.linalg import cholesky [as 别名]
def update_user_params(self):
        x_bar = np.mean(self.user_features, 0).T
        x_bar = np.reshape(x_bar, (self.num_features, 1))
        S_bar = np.cov(self.user_features.T)
        norm_X_bar = self.mu0_user - x_bar

        WI_post = inv(inv(self.WI_user) + self.num_users * S_bar + \
            np.dot(norm_X_bar, norm_X_bar.T) * \
            (self.num_users * self.beta_user) / (self.beta_user + self.num_users))

        # Not sure why we need this...
        WI_post = (WI_post + WI_post.T) / 2.0

        # update alpha_user
        self.alpha_user = sample_wishart(WI_post, self.df_post_user)

        # update mu_item
        mu_temp = (self.beta_user * self.mu0_user + self.num_users * x_bar) / \
            (self.beta_user + self.num_users)
        lam = cholesky(inv((self.beta_user + self.num_users) * self.alpha_user))
        self.mu_user = mu_temp + np.dot(lam, NormalRandom.generate_matrix(self.num_features, 1)) 

# 需要导入模块: from numpy import linalg [as 别名]
# 或者: from numpy.linalg import cholesky [as 别名]
def udpate_item_features(self):
        self.matrix = self.matrix.T
        # Gibbs sampling for item features
        for item_id in range(self.num_items):
            self.results_file.write('Item %d\n' % (item_id+1))
            users = self.matrix[:, item_id] > 0.0
            features = self.user_features[users, :]
            ratings = self.matrix[users, item_id] - self.mean_rating
            rating_len = len(ratings)
            ratings = np.reshape(ratings, (rating_len, 1))

            covar = inv(
                self.alpha_item + self.beta * np.dot(features.T, features))
            lam = cholesky(covar)
            temp = self.beta * \
                np.dot(features.T, ratings) + np.dot(
                    self.alpha_item, self.mu_item)
            mean = np.dot(covar, temp)
            temp_feature = mean + np.dot(lam, NormalRandom.generate_matrix(self.num_features, 1))
            temp_feature = np.reshape(temp_feature, (self.num_features,))
            self.item_features[item_id, :] = temp_feature 

# 需要导入模块: from numpy import linalg [as 别名]
# 或者: from numpy.linalg import cholesky [as 别名]
def transformer(self):
        """Computes a transformation matrix from the Mahalanobis matrix.
        ..math:: L = M^{1/2}

        Returns
        -------
        L : (d' x d) matrix, with d' <= d. It defines a projection. The distance can be calculated by
        ..math:: d(x,y) = \\|L(x-y)\\|_2.
        """

        if hasattr(self, 'L_'):
            return self.L_
        else:
            if hasattr(self, 'M_'):
                try:
                    L = cholesky(self.metric()).T
                    return L
                except:
                    L = metric_to_linear(self.metric())
                    return L
                # self.L_ = L
                return L
            else:
                raise NameError("Transformer was not defined. Algorithm was not fitted.") 

# 需要导入模块: from numpy import linalg [as 别名]
# 或者: from numpy.linalg import cholesky [as 别名]
def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.float64)
        # for documentation purpose:
        assert_(isinstance(res, np.ndarray))

        a = np.zeros((1, 0, 0), dtype=np.complex64).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.complex64)
        assert_(isinstance(res, np.ndarray)) 

# 需要导入模块: from numpy import linalg [as 别名]
# 或者: from numpy.linalg import cholesky [as 别名]
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 

# 需要导入模块: from numpy import linalg [as 别名]
# 或者: from numpy.linalg import cholesky [as 别名]
def CovarianceMaatrixAdaptionEvolutionStrategyF(task, epsilon=1e-20, rnd=rand):
	lam, alpha_mu, hs, sigma0 = (4 + round(3 * log(task.D))) * 10, 2, 0, 0.3 * task.bcRange()
	mu = int(round(lam / 2))
	w = log(mu + 0.5) - log(range(1, mu + 1))
	w = w / sum(w)
	mueff = 1 / sum(w ** 2)
	cs = (mueff + 2) / (task.D + mueff + 5)
	ds = 1 + cs + 2 * max(sqrt((mueff - 1) / (task.D + 1)) - 1, 0)
	ENN = sqrt(task.D) * (1 - 1 / (4 * task.D) + 1 / (21 * task.D ** 2))
	cc, c1 = (4 + mueff / task.D) / (4 + task.D + 2 * mueff / task.D), 2 / ((task.D + 1.3) ** 2 + mueff)
	cmu, hth = min(1 - c1, alpha_mu * (mueff - 2 + 1 / mueff) / ((task.D + 2) ** 2 + alpha_mu * mueff / 2)), (1.4 + 2 / (task.D + 1)) * ENN
	ps, pc, C, sigma, M = full(task.D, 0.0), full(task.D, 0.0), eye(task.D), sigma0, full(task.D, 0.0)
	x = rnd.uniform(task.bcLower(), task.bcUpper())
	x_f = task.eval(x)
	while not task.stopCondI():
		pop_step = asarray([rnd.multivariate_normal(full(task.D, 0.0), C) for _ in range(int(lam))])
		pop = asarray([task.repair(x + sigma * ps, rnd) for ps in pop_step])
		pop_f = apply_along_axis(task.eval, 1, pop)
		isort = argsort(pop_f)
		pop, pop_f, pop_step = pop[isort[:mu]], pop_f[isort[:mu]], pop_step[isort[:mu]]
		if pop_f[0] < x_f: x, x_f = pop[0], pop_f[0]
		M = sum(w * pop_step.T, axis=1)
		ps = solve(chol(C).conj() + epsilon, ((1 - cs) * ps + sqrt(cs * (2 - cs) * mueff) * M + epsilon).T)[0].T
		sigma *= exp(cs / ds * (norm(ps) / ENN - 1)) ** 0.3
		ifix = where(sigma == inf)
		if any(ifix): sigma[ifix] = sigma0
		if norm(ps) / sqrt(1 - (1 - cs) ** (2 * (task.Iters + 1))) < hth: hs = 1
		else: hs = 0
		delta = (1 - hs) * cc * (2 - cc)
		pc = (1 - cc) * pc + hs * sqrt(cc * (2 - cc) * mueff) * M
		C = (1 - c1 - cmu) * C + c1 * (tile(pc, [len(pc), 1]) * tile(pc.reshape([len(pc), 1]), [1, len(pc)]) + delta * C)
		for i in range(mu): C += cmu * w[i] * tile(pop_step[i], [len(pop_step[i]), 1]) * tile(pop_step[i].reshape([len(pop_step[i]), 1]), [1, len(pop_step[i])])
		E, V = eig(C)
		if any(E < epsilon):
			E = fmax(E, 0)
			C = lstsq(V.T, dot(V, diag(E)).T, rcond=None)[0].T
	return x, x_f 

# 需要导入模块: from numpy import linalg [as 别名]
# 或者: from numpy.linalg import cholesky [as 别名]
def cholesky_decomposition(a):
    return linalg.cholesky(a)

# Eigenvalues 

# 需要导入模块: from numpy import linalg [as 别名]
# 或者: from numpy.linalg import cholesky [as 别名]
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367 ],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp)