Python linalg.solve方法代码示例

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import solve [as 别名]
def compute_mtx_obj(GtG_lst, Tbeta_lst, Rc0, num_bands, K):
    """
    compute the matrix (M) in the objective function:
        min   c^H M c
        s.t.  c0^H c = 1

    :param GtG_lst: list of G^H * G
    :param Tbeta_lst: list of Teoplitz matrices for beta-s
    :param Rc0: right dual matrix for the annihilating filter (same of each block -> not a list)
    :return:
    """
    mtx = np.zeros((K + 1, K + 1), dtype=float)  # <= assume G, Tbeta and Rc0 are real-valued
    for loop in range(num_bands):
        Tbeta_loop = Tbeta_lst[loop]
        GtG_loop = GtG_lst[loop]
        mtx += np.dot(Tbeta_loop.T,
                      linalg.solve(np.dot(Rc0, linalg.solve(GtG_loop, Rc0.T)),
                                   Tbeta_loop)
                      )
    return mtx 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import solve [as 别名]
def mkk2ab(mk, k):
    """
    Transforms MK and M TD matrices into A and B matrices.
    Args:
        mk (numpy.ndarray): TD MK-matrix;
        k (numpy.ndarray): TD K-matrix;

    Returns:
        A and B submatrices.
    """
    if numpy.iscomplexobj(mk) or numpy.iscomplexobj(k):
        raise ValueError("MK- and/or K-matrixes are complex-valued: no transform is possible")
    m = solve(k.T, mk.T).T
    a = 0.5 * (m + k)
    b = 0.5 * (m - k)
    return a, b 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import solve [as 别名]
def compute_finite_difference_coefficients(derivative_order, grid_size):

    # from http://www.scientificpython.net/pyblog/uniform-finite-differences-all-orders

    n = 2*grid_size -1
    A = np.tile(np.arange(grid_size), (n,1)).T
    B = np.tile(np.arange(1-grid_size,grid_size), (grid_size,1))
    M = (B**A)/gamma(A+1)

    r = np.zeros(grid_size)
    r[derivative_order] = 1

    D = np.zeros((grid_size, grid_size))
    for k in range(grid_size):
        indexes = k + np.arange(grid_size)
        D[:,k] = solve(M[:,indexes], r)

    return D

#################################################################################################### 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import solve [as 别名]
def solve(self, f, tol=0):
        dx = -self.alpha*f

        n = len(self.dx)
        if n == 0:
            return dx

        df_f = np.empty(n, dtype=f.dtype)
        for k in xrange(n):
            df_f[k] = vdot(self.df[k], f)

        try:
            gamma = solve(self.a, df_f)
        except LinAlgError:
            # singular; reset the Jacobian approximation
            del self.dx[:]
            del self.df[:]
            return dx

        for m in xrange(n):
            dx += gamma[m]*(self.dx[m] + self.alpha*self.df[m])
        return dx 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import solve [as 别名]
def log_multivariate_normal_density(X, means, covars, min_covar=1.e-7):
    """Log probability for full covariance matrices. """
    if hasattr(linalg, 'solve_triangular'):
        # only in scipy since 0.9
        solve_triangular = linalg.solve_triangular
    else:
        # slower, but works
        solve_triangular = linalg.solve
    n_samples, n_dim = X.shape
    nmix = len(means)
    log_prob = np.empty((n_samples, nmix))
    for c, (mu, cv) in enumerate(zip(means, covars)):
        try:
            cv_chol = linalg.cholesky(cv, lower=True)
        except linalg.LinAlgError:
            # The model is most probabily stuck in a component with too
            # few observations, we need to reinitialize this components
            cv_chol = linalg.cholesky(cv + min_covar * np.eye(n_dim),
                                      lower=True)
        cv_log_det = 2 * np.sum(np.log(np.diagonal(cv_chol)))
        cv_sol = solve_triangular(cv_chol, (X - mu).T, lower=True).T
        log_prob[:, c] = - .5 * (np.sum(cv_sol ** 2, axis=1) + \
                                     n_dim * np.log(2 * np.pi) + cv_log_det)

    return log_prob 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import solve [as 别名]
def optimize_restpose(Tr, W, C):
    num_bones, num_verts = W.shape
    # reshape Tr according to section 4.3 into 
    # a matrix of d * P of 4-vectors
    Tr1 = Tr.reshape((-1, num_bones, 4))
    new_verts0 = []
    for j in xrange(num_verts):
        # multiply with weights and sum over num_bones axis
        Lambda = (Tr1 * W[np.newaxis,:,j,np.newaxis]).sum(axis=1)
        # convert to system matrix - last column must be one, 
        # so subtract from rhs
        gamma = C[:,j] - Lambda[:,3]
        Lambda = Lambda[:,:3]
        v = linalg.solve(np.dot(Lambda.T, Lambda), 
                         np.dot(Lambda.T, gamma))
        new_verts0.append(v)
    return np.array(new_verts0) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import solve [as 别名]
def _solve_cholesky(X, y, alpha):
    # w = inv(X^t X + alpha*Id) * X.T y
    n_samples, n_features = X.shape
    n_targets = y.shape[1]

    A = safe_sparse_dot(X.T, X, dense_output=True)
    Xy = safe_sparse_dot(X.T, y, dense_output=True)

    one_alpha = np.array_equal(alpha, len(alpha) * [alpha[0]])

    if one_alpha:
        A.flat[::n_features + 1] += alpha[0]
        return linalg.solve(A, Xy, sym_pos=True,
                            overwrite_a=True).T
    else:
        coefs = np.empty([n_targets, n_features], dtype=X.dtype)
        for coef, target, current_alpha in zip(coefs, Xy.T, alpha):
            A.flat[::n_features + 1] += current_alpha
            coef[:] = linalg.solve(A, target, sym_pos=True,
                                   overwrite_a=False).ravel()
            A.flat[::n_features + 1] -= current_alpha
        return coefs 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import solve [as 别名]
def matvec(self, f):
        dx = -f/self.alpha

        n = len(self.dx)
        if n == 0:
            return dx

        df_f = np.empty(n, dtype=f.dtype)
        for k in xrange(n):
            df_f[k] = vdot(self.df[k], f)

        b = np.empty((n, n), dtype=f.dtype)
        for i in xrange(n):
            for j in xrange(n):
                b[i,j] = vdot(self.df[i], self.dx[j])
                if i == j and self.w0 != 0:
                    b[i,j] -= vdot(self.df[i], self.df[i])*self.w0**2*self.alpha
        gamma = solve(b, df_f)

        for m in xrange(n):
            dx += gamma[m]*(self.df[m] + self.dx[m]/self.alpha)
        return dx 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import solve [as 别名]
def test_solve_equivalence():
    # For toeplitz matrices, solve_toeplitz() should be equivalent to solve().
    random = np.random.RandomState(1234)
    for n in (1, 2, 3, 10):
        c = random.randn(n)
        if random.rand() < 0.5:
            c = c + 1j * random.randn(n)
        r = random.randn(n)
        if random.rand() < 0.5:
            r = r + 1j * random.randn(n)
        y = random.randn(n)
        if random.rand() < 0.5:
            y = y + 1j * random.randn(n)

        # Check equivalence when both the column and row are provided.
        actual = solve_toeplitz((c,r), y)
        desired = solve(toeplitz(c, r=r), y)
        assert_allclose(actual, desired)

        # Check equivalence when the column is provided but not the row.
        actual = solve_toeplitz(c, b=y)
        desired = solve(toeplitz(c), y)
        assert_allclose(actual, desired) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import solve [as 别名]
def compute_b(G_lst, GtG_lst, beta_lst, Rc0, num_bands, a_ri):
    """
    compute the uniform sinusoidal samples b from the updated annihilating
    filter coeffiients.
    :param GtG_lst: list of G^H G for different subbands
    :param beta_lst: list of beta-s for different subbands
    :param Rc0: right-dual matrix, here it is the convolution matrix associated with c
    :param num_bands: number of bands
    :param L: size of b: L by 1
    :param a_ri: a 2D numpy array. each column corresponds to the measurements within a subband
    :return:
    """
    b_lst = []
    a_Gb_lst = []
    for loop in range(num_bands):
        GtG_loop = GtG_lst[loop]
        beta_loop = beta_lst[loop]
        b_loop = beta_loop - \
                 linalg.solve(GtG_loop,
                              np.dot(Rc0.T,
                                     linalg.solve(np.dot(Rc0, linalg.solve(GtG_loop, Rc0.T)),
                                                  np.dot(Rc0, beta_loop)))
                              )

        b_lst.append(b_loop)
        a_Gb_lst.append(a_ri[:, loop] - np.dot(G_lst[loop], b_loop))

    return np.column_stack(b_lst), linalg.norm(np.concatenate(a_Gb_lst)) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import solve [as 别名]
def test_care_g(self):
        A = array([[-2, -1],[-1, -1]])
        Q = array([[0, 0],[0, 1]])
        B = array([[1, 0],[0, 4]])
        R = array([[2, 0],[0, 1]])
        S = array([[0, 0],[0, 0]])
        E = array([[2, 1],[1, 2]])

        X,L,G = care(A,B,Q,R,S,E)
        # print("The solution obtained is", X)
        Gref = solve(R, B.T.dot(X).dot(E) + S.T)
        assert_array_almost_equal(Gref, G)
        assert_array_almost_equal(
            A.T.dot(X).dot(E) + E.T.dot(X).dot(A)
            - (E.T.dot(X).dot(B) + S).dot(Gref) + Q,
            zeros((2,2)))

        A = array([[-2, -1],[-1, -1]])
        Q = array([[0, 0],[0, 1]])
        B = array([[1],[0]])
        R = 1
        S = array([[1],[0]])
        E = array([[2, 1],[1, 2]])

        X,L,G = care(A,B,Q,R,S,E)
        # print("The solution obtained is", X)
        Gref = 1/R * (B.T.dot(X).dot(E) + S.T)
        assert_array_almost_equal(
            A.T.dot(X).dot(E) + E.T.dot(X).dot(A)
            - (E.T.dot(X).dot(B) + S).dot(Gref) + Q ,
            zeros((2,2)))
        assert_array_almost_equal(Gref , G) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import solve [as 别名]
def test_dare(self):
        A = array([[-0.6, 0],[-0.1, -0.4]])
        Q = array([[2, 1],[1, 0]])
        B = array([[2, 1],[0, 1]])
        R = array([[1, 0],[0, 1]])

        X,L,G = dare(A,B,Q,R)
        # print("The solution obtained is", X)
        Gref = solve(B.T.dot(X).dot(B) + R, B.T.dot(X).dot(A))
        assert_array_almost_equal(Gref, G)
        assert_array_almost_equal(
            A.T.dot(X).dot(A) - X -
            A.T.dot(X).dot(B).dot(Gref) + Q,
            zeros((2,2)))
        # check for stable closed loop
        lam = eigvals(A - B.dot(G))
        assert_array_less(abs(lam), 1.0)

        A = array([[1, 0],[-1, 1]])
        Q = array([[0, 1],[1, 1]])
        B = array([[1],[0]])
        R = 2

        X,L,G = dare(A,B,Q,R)
        # print("The solution obtained is", X)
        assert_array_almost_equal(
            A.T.dot(X).dot(A) - X -
            A.T.dot(X).dot(B) * solve(B.T.dot(X).dot(B) + R, B.T.dot(X).dot(A)) + Q, zeros((2,2)))
        assert_array_almost_equal(B.T.dot(X).dot(A) / (B.T.dot(X).dot(B) + R), G)
        # check for stable closed loop
        lam = eigvals(A - B.dot(G))
        assert_array_less(abs(lam), 1.0) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import solve [as 别名]
def test_dare_g(self):
        A = array([[-0.6, 0],[-0.1, -0.4]])
        Q = array([[2, 1],[1, 3]])
        B = array([[1, 5],[2, 4]])
        R = array([[1, 0],[0, 1]])
        S = array([[1, 0],[2, 0]])
        E = array([[2, 1],[1, 2]])

        X,L,G = dare(A,B,Q,R,S,E)
        # print("The solution obtained is", X)
        Gref = solve(B.T.dot(X).dot(B) + R, B.T.dot(X).dot(A) + S.T)
        assert_array_almost_equal(Gref,G)
        assert_array_almost_equal(
            A.T.dot(X).dot(A) - E.T.dot(X).dot(E)
            - (A.T.dot(X).dot(B) + S).dot(Gref) + Q,
            zeros((2,2)) )
        # check for stable closed loop
        lam = eigvals(A - B.dot(G), E)
        assert_array_less(abs(lam), 1.0)

        A = array([[-0.6, 0],[-0.1, -0.4]])
        Q = array([[2, 1],[1, 3]])
        B = array([[1],[2]])
        R = 1
        S = array([[1],[2]])
        E = array([[2, 1],[1, 2]])

        X,L,G = dare(A,B,Q,R,S,E)
        # print("The solution obtained is", X)
        assert_array_almost_equal(
            A.T.dot(X).dot(A) - E.T.dot(X).dot(E) -
            (A.T.dot(X).dot(B) + S).dot(solve(B.T.dot(X).dot(B) + R, B.T.dot(X).dot(A) + S.T)) + Q,
            zeros((2,2)) )
        assert_array_almost_equal((B.T.dot(X).dot(A) + S.T) / (B.T.dot(X).dot(B) + R), G)
        # check for stable closed loop
        lam = eigvals(A - B.dot(G), E)
        assert_array_less(abs(lam), 1.0) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import solve [as 别名]
def dare(A, B, Q, R, S=None, E=None, stabilizing=True):
    """ (X,L,G) = dare(A,B,Q,R) solves the discrete-time algebraic Riccati
    equation

        :math:`A^T X A - X - A^T X B (B^T X B + R)^{-1} B^T X A + Q = 0`

    where A and Q are square matrices of the same dimension. Further, Q
    is a symmetric matrix. The function returns the solution X, the gain
    matrix G = (B^T X B + R)^-1 B^T X A and the closed loop eigenvalues L,
    i.e., the eigenvalues of A - B G.

    (X,L,G) = dare(A,B,Q,R,S,E) solves the generalized discrete-time algebraic
    Riccati equation

        :math:`A^T X A - E^T X E - (A^T X B + S) (B^T X B + R)^{-1} (B^T X A + S^T) + Q = 0`

    where A, Q and E are square matrices of the same dimension. Further, Q and
    R are symmetric matrices. The function returns the solution X, the gain
    matrix :math:`G = (B^T X B + R)^{-1} (B^T X A + S^T)` and the closed loop
    eigenvalues L, i.e., the eigenvalues of A - B G , E.
    """
    if S is not None or E is not None or not stabilizing:
        return dare_old(A, B, Q, R, S, E, stabilizing)
    else:
        Rmat = _ssmatrix(R)
        Qmat = _ssmatrix(Q)
        X = solve_discrete_are(A, B, Qmat, Rmat)
        G = solve(B.T.dot(X).dot(B) + Rmat, B.T.dot(X).dot(A))
        L = eigvals(A - B.dot(G))
        return _ssmatrix(X), L, _ssmatrix(G) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import solve [as 别名]
def get_multivariate_ll(contexts, gmpe, imt):
    """
    Returns the multivariate loglikelihood, as described om equation 7 of
    Mak et al. (2017)
    """
    observations, v_mat, expected_mat, neqs, nrecs = _build_matrices(
        contexts, gmpe, imt)
    sign, logdetv = np.linalg.slogdet(v_mat)
    b_mat = observations - expected_mat
    return (float(nrecs) * np.log(2.0 * np.pi) + logdetv +
            (b_mat.T.dot(solve(v_mat, b_mat)))) / 2.