python类matrix_rank()的实例源码

def test_matrix_rank(self):
        # Full rank matrix
        yield assert_equal, 4, matrix_rank(np.eye(4))
        # rank deficient matrix
        I = np.eye(4)
        I[-1, -1] = 0.
        yield assert_equal, matrix_rank(I), 3
        # All zeros - zero rank
        yield assert_equal, matrix_rank(np.zeros((4, 4))), 0
        # 1 dimension - rank 1 unless all 0
        yield assert_equal, matrix_rank([1, 0, 0, 0]), 1
        yield assert_equal, matrix_rank(np.zeros((4,))), 0
        # accepts array-like
        yield assert_equal, matrix_rank([1]), 1
        # greater than 2 dimensions raises error
        yield assert_raises, TypeError, matrix_rank, np.zeros((2, 2, 2))
        # works on scalar
        yield assert_equal, matrix_rank(1), 1

def test_matrix_rank(self):
        # Full rank matrix
        yield assert_equal, 4, matrix_rank(np.eye(4))
        # rank deficient matrix
        I = np.eye(4)
        I[-1, -1] = 0.
        yield assert_equal, matrix_rank(I), 3
        # All zeros - zero rank
        yield assert_equal, matrix_rank(np.zeros((4, 4))), 0
        # 1 dimension - rank 1 unless all 0
        yield assert_equal, matrix_rank([1, 0, 0, 0]), 1
        yield assert_equal, matrix_rank(np.zeros((4,))), 0
        # accepts array-like
        yield assert_equal, matrix_rank([1]), 1
        # greater than 2 dimensions raises error
        yield assert_raises, TypeError, matrix_rank, np.zeros((2, 2, 2))
        # works on scalar
        yield assert_equal, matrix_rank(1), 1

def test_matrix_rank(self):
        # Full rank matrix
        yield assert_equal, 4, matrix_rank(np.eye(4))
        # rank deficient matrix
        I = np.eye(4)
        I[-1, -1] = 0.
        yield assert_equal, matrix_rank(I), 3
        # All zeros - zero rank
        yield assert_equal, matrix_rank(np.zeros((4, 4))), 0
        # 1 dimension - rank 1 unless all 0
        yield assert_equal, matrix_rank([1, 0, 0, 0]), 1
        yield assert_equal, matrix_rank(np.zeros((4,))), 0
        # accepts array-like
        yield assert_equal, matrix_rank([1]), 1
        # greater than 2 dimensions raises error
        yield assert_raises, TypeError, matrix_rank, np.zeros((2, 2, 2))
        # works on scalar
        yield assert_equal, matrix_rank(1), 1

def test_matrix_rank(self):
        # Full rank matrix
        yield assert_equal, 4, matrix_rank(np.eye(4))
        # rank deficient matrix
        I = np.eye(4)
        I[-1, -1] = 0.
        yield assert_equal, matrix_rank(I), 3
        # All zeros - zero rank
        yield assert_equal, matrix_rank(np.zeros((4, 4))), 0
        # 1 dimension - rank 1 unless all 0
        yield assert_equal, matrix_rank([1, 0, 0, 0]), 1
        yield assert_equal, matrix_rank(np.zeros((4,))), 0
        # accepts array-like
        yield assert_equal, matrix_rank([1]), 1
        # greater than 2 dimensions raises error
        yield assert_raises, TypeError, matrix_rank, np.zeros((2, 2, 2))
        # works on scalar
        yield assert_equal, matrix_rank(1), 1

def test_matrix_rank(self):
        # Full rank matrix
        yield assert_equal, 4, matrix_rank(np.eye(4))
        # rank deficient matrix
        I = np.eye(4)
        I[-1, -1] = 0.
        yield assert_equal, matrix_rank(I), 3
        # All zeros - zero rank
        yield assert_equal, matrix_rank(np.zeros((4, 4))), 0
        # 1 dimension - rank 1 unless all 0
        yield assert_equal, matrix_rank([1, 0, 0, 0]), 1
        yield assert_equal, matrix_rank(np.zeros((4,))), 0
        # accepts array-like
        yield assert_equal, matrix_rank([1]), 1
        # greater than 2 dimensions raises error
        yield assert_raises, TypeError, matrix_rank, np.zeros((2, 2, 2))
        # works on scalar
        yield assert_equal, matrix_rank(1), 1

def qr_fast_givens(A):
    '''
    QR factorization via fast givens transformation
    '''
    m = A.shape[0]
    n = A.shape[1]
    d = np.ones(m)
    for j in range(n):
        for i in range(m-1, j, -1):
            alpha, beta, ty = givens.fast_givens(A[i-1:i+1, j], d[i-1:i+1])
            if ty == 1:
                A[i-1:i+1, j:n+1] = np.array([[beta, 1],
                                              [1, alpha]]).dot(A[i-1:i+1, j:n+1])
            else:
                A[i-1:i+1, j:n+1] = np.array([[1, alpha],
                                              [beta, 1]]).dot(A[i-1:i+1, j:n+1])

    return A

# least square using QR (A must be full column rank)

def qr_ls(A, b):
    '''
    least square using QR (A must be full column rank)
    '''
    m = A.shape[0]
    n = A.shape[1]
    if rank(A) < n:
        raise Exception('Rank deficient')

    A = qr_householder(A)
    for j in range(n):
        v = np.hstack((1, A[j+1:, j]))
        A[j+1:, j] = 0
        b[j:] = (np.eye(m - j + 1) - 2 * np.outer(v, v) / la.norm(v, 2)).dot(b[j:])

    x_ls = la.solve(A[:n, :n], b[:n])

    return x_ls

def fullrank(X):
    '''
    return full-rank decomposition of X = FG^T
    '''
    rankX = rank(X)

    U, eigvals, Vh = la.svd(X)

    #construct a r-rank sigma-square-root matrix

    sigma = np.eye(rankX)
    for i in range(sigma.shape[0]):
        sigma[i, i] = np.sqrt(eigvals[i])

    F = U.dot(np.vstack((sigma, np.zeros((X.shape[0] - rankX, rankX)))))
    Gh = np.hstack((sigma, np.zeros((rankX, X.shape[1] - rankX)))).dot(Vh)

    return F, Gh

def hausdorff_distance(X, Y, n):
    '''
    distace between subspaces by Hausdorff's definition
    '''
    if rank(X) != X.shape[1] & rank(Y) != Y.shape[1]:
        raise Exception('Please provide subspaces with full COLUMN rank')

    inner = 0

    for i in range(X.shape[1]):
        for j in range(Y.shape[1]):
            inner = inter + np.square(X[:, i].conjugate().T.dot(Y[:, j]))

    distance = np.sqrt(np.max(rank(X), rank(Y)) - inner)

    return distance

# distance with inner-product

def kernel_distance(X, Y, n):
    '''
    distance with inner-product
    '''
    if rank(X) != X.shape[1] & rank(Y) != Y.shape[1]:
        raise Exception('Please provide subspaces with full COLUMN rank')

    inner = 0

    for i in range(X.shape[1]):
        for j in range(Y.shape[1]):
            inter = inter + np.square(X[:, i].conjugate().T.dot(Y[:, j]))

    distance = np.sqrt(inner)


# return the dimension of the intersection of two subspaces

def ls_qr(A, b):
    '''
    least square using QR (A must be full column rank)
    '''
    m = A.shape[0]
    n = A.shape[1]
    if rank(A) < n:
        raise Exception('Rank deficient')

    A = qr.qr_householder(A)
    for j in range(n):
        v = np.hstack((1, A[j+1:, j]))
        A[j+1:, j] = 0
        b[j:] = (np.eye(m - j + 1) - 2 * np.outer(v, v) / la.norm(v, 2)).dot(b[j:])

    x_ls = la.solve(A, b)

    return x_ls0

def test_matrix_rank(self):
        # Full rank matrix
        yield assert_equal, 4, matrix_rank(np.eye(4))
        # rank deficient matrix
        I = np.eye(4)
        I[-1, -1] = 0.
        yield assert_equal, matrix_rank(I), 3
        # All zeros - zero rank
        yield assert_equal, matrix_rank(np.zeros((4, 4))), 0
        # 1 dimension - rank 1 unless all 0
        yield assert_equal, matrix_rank([1, 0, 0, 0]), 1
        yield assert_equal, matrix_rank(np.zeros((4,))), 0
        # accepts array-like
        yield assert_equal, matrix_rank([1]), 1
        # greater than 2 dimensions raises error
        yield assert_raises, TypeError, matrix_rank, np.zeros((2, 2, 2))
        # works on scalar
        yield assert_equal, matrix_rank(1), 1

def _validate_inputs(self):
        x, z = self._x, self._z
        if x.shape[1] == 0:
            raise ValueError('Model must contain at least one regressor.')
        if self.instruments.shape[1] < self.endog.shape[1]:
            raise ValueError('The number of instruments ({0}) must be at least '
                             'as large as the number of endogenous regressors'
                             ' ({1}).'.format(self.instruments.shape[1],
                                              self.endog.shape[1]))
        if matrix_rank(x) < x.shape[1]:
            raise ValueError('regressors [exog endog] do not have full '
                             'column rank')
        if matrix_rank(z) < z.shape[1]:
            raise ValueError('instruments [exog instruments]  do not have '
                             'full column rank')
        self._has_constant, self._const_loc = has_constant(x)

def gen_report(M, name, obj, err, L_test, S_test, L_true, S_true,
               y_true):

    lam = 1.0/np.sqrt(np.max(M.shape))
    nn = norm(L_test, 'nuc')
    on = np.sum(np.abs(S_test))
    o = nn + lam*on
    print('Rank = %d, NNZs = %d' % (matrix_rank(L_test),
                                    np.count_nonzero(S_test)))
    print('Nuclear Norm = %e' % nn)
    print('One Norm = %e' % on)
    print('Objective = %e' % o)
    if L_true is not None:
        print('Recovery Error = %e' %
              (norm(L_test - L_true, 'fro')/norm(L_true, 'fro')))

    if y_true is not None:
        y_test = np.linalg.norm(S_test, axis=1)
        tp, fp, _ = metrics.roc_curve(y_true, y_test)
        score = metrics.roc_auc_score(y_true, y_test)
        auc_ax.plot(tp, fp, label=name + ' AUC=' + str(score))
    obj_ax.plot(obj, label=name + ' Objective')

def test_matrix_rank(self):
        # Full rank matrix
        yield assert_equal, 4, matrix_rank(np.eye(4))
        # rank deficient matrix
        I = np.eye(4)
        I[-1, -1] = 0.
        yield assert_equal, matrix_rank(I), 3
        # All zeros - zero rank
        yield assert_equal, matrix_rank(np.zeros((4, 4))), 0
        # 1 dimension - rank 1 unless all 0
        yield assert_equal, matrix_rank([1, 0, 0, 0]), 1
        yield assert_equal, matrix_rank(np.zeros((4,))), 0
        # accepts array-like
        yield assert_equal, matrix_rank([1]), 1
        # greater than 2 dimensions raises error
        yield assert_raises, TypeError, matrix_rank, np.zeros((2, 2, 2))
        # works on scalar
        yield assert_equal, matrix_rank(1), 1

def get_R_rank(self):
        """Get the rank of the R matrix.

        Returns
        -------

        rank : int
            The rank of the R matrix
        """
        return matrix_rank(self.get_R())

def test_reduced_rank():
    # Test matrices with reduced rank
    rng = np.random.RandomState(20120714)
    for i in range(100):
        # Make a rank deficient matrix
        X = rng.normal(size=(40, 10))
        X[:, 0] = X[:, 1] + X[:, 2]
        # Assert that matrix_rank detected deficiency
        assert_equal(matrix_rank(X), 9)
        X[:, 3] = X[:, 4] + X[:, 5]
        assert_equal(matrix_rank(X), 8)

def test_reduced_rank():
    # Test matrices with reduced rank
    rng = np.random.RandomState(20120714)
    for i in range(100):
        # Make a rank deficient matrix
        X = rng.normal(size=(40, 10))
        X[:, 0] = X[:, 1] + X[:, 2]
        # Assert that matrix_rank detected deficiency
        assert_equal(matrix_rank(X), 9)
        X[:, 3] = X[:, 4] + X[:, 5]
        assert_equal(matrix_rank(X), 8)

def test_reduced_rank():
    # Test matrices with reduced rank
    rng = np.random.RandomState(20120714)
    for i in range(100):
        # Make a rank deficient matrix
        X = rng.normal(size=(40, 10))
        X[:, 0] = X[:, 1] + X[:, 2]
        # Assert that matrix_rank detected deficiency
        assert_equal(matrix_rank(X), 9)
        X[:, 3] = X[:, 4] + X[:, 5]
        assert_equal(matrix_rank(X), 8)

def test_reduced_rank():
    # Test matrices with reduced rank
    rng = np.random.RandomState(20120714)
    for i in range(100):
        # Make a rank deficient matrix
        X = rng.normal(size=(40, 10))
        X[:, 0] = X[:, 1] + X[:, 2]
        # Assert that matrix_rank detected deficiency
        assert_equal(matrix_rank(X), 9)
        X[:, 3] = X[:, 4] + X[:, 5]
        assert_equal(matrix_rank(X), 8)

def assert_invertible(cls, A):
        rank = matrix_rank(A)
        is_invertible = rank != 0

        cls.assertTrue(is_invertible)

def assert_invertible(self, A):
        rank = matrix_rank(A)
        is_invertible = rank != 0

        self.assertTrue(is_invertible)

def test_cp_tensor():
    """test for random.cp_tensor"""
    shape = (10, 11, 12)
    rank = 4

    tensor = cp_tensor(shape, rank, full=True)
    for i in range(T.ndim(tensor)):
        T.assert_equal(matrix_rank(T.to_numpy(unfold(tensor, i))), rank)

    factors = cp_tensor(shape, rank, full=False)
    for i, factor in enumerate(factors):
        T.assert_equal(factor.shape, (shape[i], rank),
                err_msg=('{}-th factor has shape {}, expected {}'.format(
                     i, factor.shape, (shape[i], rank))))

def test_tucker_tensor():
    """test for random.tucker_tensor"""
    shape = (10, 11, 12)
    rank = 4

    tensor = tucker_tensor(shape, rank, full=True)
    for i in range(T.ndim(tensor)):
        T.assert_equal(matrix_rank(T.to_numpy(unfold(tensor, i))), rank)

    core, factors = tucker_tensor(shape, rank, full=False)
    for i, factor in enumerate(factors):
        T.assert_equal(factor.shape, (shape[i], rank),
                err_msg=('{}-th factor has shape {}, expected {}'.format(
                     i, factor.shape, (shape[i], rank))))

    shape = (10, 11, 12)
    rank = (6, 4, 5)
    tensor = tucker_tensor(shape, rank, full=True)
    for i in range(T.ndim(tensor)):
        T.assert_equal(matrix_rank(T.to_numpy(unfold(tensor, i))),  min(shape[i], rank[i]))

    core, factors = tucker_tensor(shape, rank, full=False)
    for i, factor in enumerate(factors):
        T.assert_equal(factor.shape, (shape[i], rank[i]),
                err_msg=('{}-th factor has shape {}, expected {}.'.format(
                     i, factor.shape, (shape[i], rank[i]))))
    T.assert_equal(core.shape, rank, err_msg='core has shape {}, expected {}.'.format(
                                     core.shape, rank))
    for factor in factors:
        T.assert_array_almost_equal(T.dot(T.transpose(factor), factor), T.tensor(np.eye(factor.shape[1])))
    tensor = tucker_to_tensor(core, factors)
    reconstructed = multi_mode_dot(tensor, factors, transpose=True)
    T.assert_array_almost_equal(core, reconstructed)

    with T.assert_raises(ValueError):
        tucker_tensor((3, 4, 5), (3, 6, 3))

def test_reduced_rank():
    # Test matrices with reduced rank
    rng = np.random.RandomState(20120714)
    for i in range(100):
        # Make a rank deficient matrix
        X = rng.normal(size=(40, 10))
        X[:, 0] = X[:, 1] + X[:, 2]
        # Assert that matrix_rank detected deficiency
        assert_equal(matrix_rank(X), 9)
        X[:, 3] = X[:, 4] + X[:, 5]
        assert_equal(matrix_rank(X), 8)

def test_grams():
    sys = 0.6*Alpha(0.01) + 0.4*Lowpass(0.05)

    A, B, C, D = sys2ss(sys)

    P = control_gram(sys)
    assert np.allclose(np.dot(A, P) + np.dot(P, A.T), -np.dot(B, B.T))
    assert matrix_rank(P) == len(P)  # controllable

    Q = observe_gram(sys)
    assert np.allclose(np.dot(A.T, Q) + np.dot(Q, A), -np.dot(C.T, C))
    assert matrix_rank(Q) == len(Q)  # observable

def projection(X, Y):

    rankX = rank(X)
    rankY = rank(Y)

    # rank, or dimension, or the original space
    rankO = rankX + rankY

    # check if two subspaces have the same shapes
    if X.shape != Y.shape:
        raise Exception('The two subspaces do not have the same shapes')

    # check if O is singular
    if la.det(np.hstack((X, Y))) == 0:
        raise Exception('X + Y is not the direct sum of the original space')

    # check whether each subspace is of full column/row rank
    if rankX < min(X.shape):
        raise Exception('subspace X is not of full rank')

    elif rankY < min(Y.shape):
        raise Exception('subspace Y is not of full rank')
    # X and Y are of full column rank
    elif rankX == X.shape[1] & rankY == Y.shape[1]:
        return np.hstack((X, np.zeros((X.shape[0], rankO - rankX)))).dot(la.inv(np.hstack((X, Y))))
    # X and Y are of full row rank
    elif rankX == X.shape[0] & rankY == Y.shape[0]:
        return np.vstack((X, np.zeros((rankO - rankX, X.shape[1])))).dot(la.inv(np.vstack(X, Y)))

# orthogonal projection matrix

def orthoProjection(X, n):

    # check if X is a subspace of the original space
    if rank(X) < n:
        P = X.dot(la.inv(X.conjugate().T.dot(X))).dot(X.conjugate().T)

        # return: orthogonal projection onto subspace X, orthogonal projection X's orthogonal complement subspace
        return P, np.eye(P.shape[0]) - P
    else:
        raise Exception('not a subspace')

# projection transformation from original space onto its subspace X along its completement Y

def subspace_intersection(X, Y, n):
    '''
    return the dimension of the intersection of two subspaces
    '''
    U = principal_angles(X, Y, n)[1]
    V = principal_angles(X, Y, n)[2]

    return rank(np.hstack(U, V))

# distance between A and any lower rank matrix

def lowRank_distance(A, k):
    '''
    distance between A and any lower rank matrix
    '''
    if rank(A) >= k:
        raise Exception('Please provide a lower rank k')

    sigma = la.svdvals(A)
    # return the k+1'th singular value
    return sigma[k]