Python linalg.fractional_matrix_power函数代码示例

    def test_type_preservation_and_conversion(self):
        # The fractional_matrix_power matrix function should preserve
        # the type of a matrix whose eigenvalues
        # are positive with zero imaginary part.
        # Test this preservation for variously structured matrices.
        complex_dtype_chars = ('F', 'D', 'G')
        for matrix_as_list in (
                [[1, 0], [0, 1]],
                [[1, 0], [1, 1]],
                [[2, 1], [1, 1]],
                [[2, 3], [1, 2]]):

            # check that the spectrum has the expected properties
            W = scipy.linalg.eigvals(matrix_as_list)
            assert_(not any(w.imag or w.real < 0 for w in W))

            # Check various positive and negative powers
            # with absolute values bigger and smaller than 1.
            for p in (-2.4, -0.9, 0.2, 3.3):

                # check float type preservation
                A = np.array(matrix_as_list, dtype=float)
                A_power = fractional_matrix_power(A, p)
                assert_(A_power.dtype.char not in complex_dtype_chars)

                # check complex type preservation
                A = np.array(matrix_as_list, dtype=complex)
                A_power = fractional_matrix_power(A, p)
                assert_(A_power.dtype.char in complex_dtype_chars)

                # check float->complex for the matrix negation
                A = -np.array(matrix_as_list, dtype=float)
                A_power = fractional_matrix_power(A, p)
                assert_(A_power.dtype.char in complex_dtype_chars)

	def fidelity(self, m, n):
		"""Compute the fidelity between the density matrces m and n.

		:param numpy_array m: Density matrix
		:param numpy_array n: Density matrix

		:return: The fideltiy between m and n (:math:`\mathrm{Tr}(\sqrt{\sqrt{m}n\sqrt{m}})^2`).
		:rtype: complex
		"""
		sqrt_m = fractional_matrix_power(m, 0.5)
		F = np.trace(fractional_matrix_power(np.dot(sqrt_m,np.dot(n, sqrt_m)), 0.5))**2.0

		return F

    def get_mean_and_covariance(self, w2v, num_of_occurences):
        """ get mean and covariance of words vectors over the training set of word2vec model
            w2v -- word2vec model (in matrix form)
            num_of_occurences -- array that specifies weights for averaging over words
        """
        weights = num_of_occurences/np.sum(num_of_occurences)
        try:
            w2v_temp = np.multiply(w2v, weights)
        except MemoryError:
            w2v_temp = np.copy(w2v)
            for wn in range(np.shape(w2v)[1]):
                w2v_temp[:, wn] *= weights[wn]
        self.Mean = np.sum(w2v_temp, 1)
        try:
            w2v_except0 = w2v - np.reshape(self.Mean, (len(self.Mean), 1))
        except MemoryError:
            w2v_except0 = w2v_temp  # just to set the right shape (to avoid memoryError)
            for wn in range(np.shape(w2v)[1]):
                w2v_except0[:, wn] = w2v[:, wn] - self.Mean

        try:
            w2v_normalized = np.multiply(w2v_except0, np.power(weights, 0.5))
        except MemoryError:
            w2v_normalized = w2v_except0
            for wn in range(np.shape(w2v_except0)[1]):
                w2v_normalized[:, wn] *= weights[wn]**0.5

        self.Cov = np.dot(w2v_normalized, np.transpose(w2v_normalized))
        self.Cov = self.Cov/np.shape(w2v)[1]
        self.SqrtCov = fractional_matrix_power(self.Cov, -0.5)

 def test_larger_abs_fractional_matrix_powers(self):
     np.random.seed(1234)
     for n in (2, 3, 5):
         for i in range(10):
             M = np.random.randn(n, n) + 1j * np.random.randn(n, n)
             M_one_fifth = fractional_matrix_power(M, 0.2)
             # Test the round trip.
             M_round_trip = np.linalg.matrix_power(M_one_fifth, 5)
             assert_allclose(M, M_round_trip)
             # Test a large abs fractional power.
             X = fractional_matrix_power(M, -5.4)
             Y = np.linalg.matrix_power(M_one_fifth, -27)
             assert_allclose(X, Y)
             # Test another large abs fractional power.
             X = fractional_matrix_power(M, 3.8)
             Y = np.linalg.matrix_power(M_one_fifth, 19)
             assert_allclose(X, Y)

    def test_singular(self):
        # Negative fractional powers do not work with singular matrices.
        for matrix_as_list in (
                [[0, 0], [0, 0]],
                [[1, 1], [1, 1]],
                [[1, 2], [3, 6]],
                [[0, 0, 0], [0, 1, 1], [0, -1, 1]]):

            # Check fractional powers both for float and for complex types.
            for newtype in (float, complex):
                A = np.array(matrix_as_list, dtype=newtype)
                for p in (-0.7, -0.9, -2.4, -1.3):
                    A_power = fractional_matrix_power(A, p)
                    assert_(np.isnan(A_power).all())
                for p in (0.2, 1.43):
                    A_power = fractional_matrix_power(A, p)
                    A_round_trip = fractional_matrix_power(A_power, 1/p)
                    assert_allclose(A_round_trip, A)

 def test_round_trip_random_complex(self):
     np.random.seed(1234)
     for p in range(1, 5):
         for n in range(1, 5):
             M_unscaled = np.random.randn(n, n) + 1j * np.random.randn(n, n)
             for scale in np.logspace(-4, 4, 9):
                 M = M_unscaled * scale
                 M_root = fractional_matrix_power(M, 1/p)
                 M_round_trip = np.linalg.matrix_power(M_root, p)
                 assert_allclose(M_round_trip, M)

    def test_al_mohy_higham_2012_experiment_1(self):
        # Fractional powers of a tricky upper triangular matrix.
        A = _get_al_mohy_higham_2012_experiment_1()

        # Test remainder matrix power.
        A_funm_sqrt, info = funm(A, np.sqrt, disp=False)
        A_sqrtm, info = sqrtm(A, disp=False)
        A_rem_power = _matfuncs_inv_ssq._remainder_matrix_power(A, 0.5)
        A_power = fractional_matrix_power(A, 0.5)
        assert_array_equal(A_rem_power, A_power)
        assert_allclose(A_sqrtm, A_power)
        assert_allclose(A_sqrtm, A_funm_sqrt)

        # Test more fractional powers.
        for p in (1/2, 5/3):
            A_power = fractional_matrix_power(A, p)
            A_round_trip = fractional_matrix_power(A_power, 1/p)
            assert_allclose(A_round_trip, A, rtol=1e-2)
            assert_allclose(np.tril(A_round_trip, 1), np.tril(A, 1))

    def test_type_conversion_mixed_sign_or_complex_spectrum(self):
        complex_dtype_chars = ("F", "D", "G")
        for matrix_as_list in ([[1, 0], [0, -1]], [[0, 1], [1, 0]], [[0, 1, 0], [0, 0, 1], [1, 0, 0]]):

            # check that the spectrum has the expected properties
            W = scipy.linalg.eigvals(matrix_as_list)
            assert_(any(w.imag or w.real < 0 for w in W))

            # Check various positive and negative powers
            # with absolute values bigger and smaller than 1.
            for p in (-2.4, -0.9, 0.2, 3.3):

                # check complex->complex
                A = np.array(matrix_as_list, dtype=complex)
                A_power = fractional_matrix_power(A, p)
                assert_(A_power.dtype.char in complex_dtype_chars)

                # check float->complex
                A = np.array(matrix_as_list, dtype=float)
                A_power = fractional_matrix_power(A, p)
                assert_(A_power.dtype.char in complex_dtype_chars)

 def fit(self,X,Y):
     self._compute_covariance(X,Y)
     S_11_ = fractional_matrix_power(self.S_11,-0.5)
     
     S_22_ = fractional_matrix_power(self.S_22,-0.5)
     
     T = np.dot(np.dot(S_11_,self.S_12),S_22_)
     
     U ,S ,V = np.linalg.svd(T)
     
     self.U = U[:,:self.k]
     self.S = S[:self.k]
     self.V = V[:self.k,:]
     
     self.A = np.dot(S_11_,U)
     self.B = np.dot(S_22_,V.T)
     self.A = self.A[:,:self.k]
     self.B = self.B[:,:self.k]
     
     self.coeff_ = np.dot(self.A,self.V)
     
     return self

 def test_round_trip_random_float(self):
     # This test is more annoying because it can hit the branch cut;
     # this happens when the matrix has an eigenvalue
     # with no imaginary component and with a real negative component,
     # and it means that the principal branch does not exist.
     np.random.seed(1234)
     for p in range(1, 5):
         for n in range(1, 5):
             M_unscaled = np.random.randn(n, n)
             for scale in np.logspace(-4, 4, 9):
                 M = M_unscaled * scale
                 M_root = fractional_matrix_power(M, 1/p)
                 M_round_trip = np.linalg.matrix_power(M_root, p)
                 assert_allclose(M_round_trip, M)

    def test_singular(self):
        # Negative fractional powers do not work with singular matrices.
        # Neither do non-integer fractional powers,
        # because the scaling and squaring cannot deal with it.
        for matrix_as_list in (
                [[0, 0], [0, 0]],
                [[1, 1], [1, 1]],
                [[1, 2], [3, 6]],
                [[0, 0, 0], [0, 1, 1], [0, -1, 1]]):

            # check that the spectrum has the expected properties
            W = scipy.linalg.eigvals(matrix_as_list)
            assert_(np.count_nonzero(W) < len(W))

            # check fractional powers both for float and for complex types
            for newtype in (float, complex):
                A = np.array(matrix_as_list, dtype=newtype)
                for p in (-0.7, -0.9, -2.4, -1.3):
                    A_power = fractional_matrix_power(A, p)
                    assert_(np.isnan(A_power).all())
                for p in (0.2, 1.43):
                    A_power = fractional_matrix_power(A, p)
                    A_round_trip = fractional_matrix_power(A_power, 1/p)
                    assert_allclose(A_round_trip, A)

    def test_random_matrices_and_powers(self):
        # Each independent iteration of this fuzz test picks random parameters.
        # It tries to hit some edge cases.
        np.random.seed(1234)
        nsamples = 20
        for i in range(nsamples):
            # Sample a matrix size and a random real power.
            n = random.randrange(1, 5)
            p = np.random.randn()

            # Sample a random real or complex matrix.
            matrix_scale = np.exp(random.randrange(-4, 5))
            A = np.random.randn(n, n)
            if random.choice((True, False)):
                A = A + 1j * np.random.randn(n, n)
            A = A * matrix_scale

            # Check a couple of analytically equivalent ways
            # to compute the fractional matrix power.
            # These can be compared because they both use the principal branch.
            A_power = fractional_matrix_power(A, p)
            A_logm, info = logm(A, disp=False)
            A_power_expm_logm = expm(A_logm * p)
            assert_allclose(A_power, A_power_expm_logm)

temp = np.array(xrange(1829))+1
adDup = np.hstack((temp[:,np.newaxis],adMat))
temp = np.array(xrange(1830))
adDup = np.vstack((temp,adDup))
np.savetxt("aff.csv", adDup, delimiter=';')

# social similarity matrix
socialsim = sklearn.metrics.pairwise.pairwise_kernels(social,metric='rbf',gamma=0.8)

adMat = adMat+socialsim

# Make D matrix
D = np.diag(adMat.sum(1))

# Make laplacian
Dinv = la.fractional_matrix_power(D,-0.5)
L = np.eye(numProj) - np.dot(np.dot(Dinv,adMat),Dinv)

#L = D - adMat

# Cal eigenvector of laplacian
eigval,eigvec = np.linalg.eigh(L)

# Sort from smallest to largest
idx = np.argsort(eigval)
eigval = eigval[idx]
eigvec = eigvec[:,idx]

 Perform PCA
uXPca,sXPca,wXPca = la.svd(eigvec, full_matrices=False)
K=500

def sq_IG_distance(cov_1, cov_2) :
    
    cov_1_pow = fractional_matrix_power(cov_1, -0.5)
    return norm(logm(np.linalg.multi_dot([cov_1_pow, cov_2, cov_1_pow])), ord='fro') ** 2

#adjMat[2,9] = 1
#adjMat[5,0] = 1
#adjMat[0,5] = 1
#adjMat[1,8] = 1
#adjMat[8,1] = 1

#print adjMat
#print(np.sum(adjMat))
#show_graph(adjMat)

# Create diagonal matrix
diagMat = np.diag(adjMat.sum(axis=0))
#print diagMat

# Create Laplacian matrix
diagMatinv = la.fractional_matrix_power(diagMat,-0.5)

lapMat = np.eye(2*N+N2) - np.dot(np.dot(diagMatinv,adjMat),diagMatinv)

# Unnormalized
#lapMat = diagMat - adjMat

#print lapMat

# Cal eigenvector of laplacian
eigval,eigvec = np.linalg.eig(lapMat)

# Sort from smallest to largest
idx = np.argsort(eigval)
eigval = eigval[idx]
eigvec = eigvec[:,idx]