我们从Python开源项目中，提取了以下40个代码示例，用于说明如何使用numpy.linalg.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 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 _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 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 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_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]
def has_constant(x): """ Parameters ---------- x: ndarray Array to be checked for a constant (n,k) Returns ------- const : bool Flag indicating whether x contains a constant or has column span with a constant loc : int Column location of constant """ if np.any(np.all(x == 1, axis=0)): loc = np.argwhere(np.all(x == 1, axis=0)) return True, int(loc) if np.any((np.ptp(x, axis=0) == 0) & ~np.all(x == 0, axis=0)): loc = np.any((np.ptp(x, axis=0) == 0) & ~np.all(x == 0, axis=0)) loc = np.argwhere(loc) return True, int(loc) n = x.shape[0] aug_rank = matrix_rank(np.c_[np.ones((n, 1)), x]) rank = matrix_rank(x) has_const = bool(aug_rank == rank) loc = None if has_const: out = np.linalg.lstsq(x, np.ones((n, 1))) beta = out[0].ravel() loc = np.argmax(np.abs(beta) * x.var(0)) return has_const, loc
def _validate_data(self): """Check input shape and remove missing""" y = self._y = self.dependent.values2d x = self._x = self.exog.values2d w = self._w = self.weights.values2d if y.shape[0] != x.shape[0]: raise ValueError('dependent and exog must have the same number of ' 'observations.') if y.shape[0] != w.shape[0]: raise ValueError('weights must have the same number of ' 'observations as dependent.') all_missing = np.any(np.isnan(y), axis=1) & np.all(np.isnan(x), axis=1) missing = (np.any(np.isnan(y), axis=1) | np.any(np.isnan(x), axis=1) | np.any(np.isnan(w), axis=1)) missing_warning(all_missing ^ missing) if np.any(missing): self.dependent.drop(missing) self.exog.drop(missing) self.weights.drop(missing) x = self.exog.values2d self._not_null = ~missing w = self.weights.dataframe if np.any(w.values <= 0): raise ValueError('weights must be strictly positive.') w = w / w.mean() self.weights = PanelData(w) if matrix_rank(x) < x.shape[1]: raise ValueError('exog does not have full column rank.') self._constant, self._constant_index = has_constant(x)
def train(x,y): """ Build the linear least weight vector W :param x: NxD matrix containing N attributes vectors for training :param y: NxK matrix containing N class vectors for training """ # D = Number of attributes D = x.shape[1] + 1 # K = Number of classes K = y.shape[1] # Build the sums of xi*xi' and xi*yi' sum1 = np.zeros((D,D)) # init placeholder sum2 = np.zeros((D,K)) i = 0 for x_i in x: # loop over all vectors x_i = np.append(1, x_i) # augment vector with a 1 y_i = y[i] sum1 += np.outer(x_i, x_i) # find xi*xi' sum2 += np.outer(x_i, y_i) # find xi*yi' i += 1 # Check that condition number is finite # and therefore sum1 is nonsingular (invertable) while matrix_rank(sum1) != D: # Naive choice of sigma. # Could cause inaccuracies when sum1 has small values # However, in most cases the matrix WILL be invertable sum1 = sum1 + 0.001 * np.eye(D) # Return weight vector # Weight vector multiplies sums and inverse of sum1 return np.dot(inv(sum1),sum2)
def ils(B, y, p=1): m, n = B.shape if rank(B) < n: raise ValueError("Matrix is rank deficient!") # Reduction R, Z_r, y_r = reduction(B.copy(), y.copy()) # Search _z_hat = search(R.copy(), y_r[:n].copy(), p) return dot(Z_r, _z_hat)
def mils(A, B, y, p=1): # x_hat,z_hat = mils(A,B,y,p) produces p pairs of optimal solutions to # the mixed integer least squares problem min_{x,z}||y-Ax-Bz||, # where x and z are real and integer vectors, respectively. # # Input arguments: # A - m by k real matrix # B - m by n real matrix # [A,B] has full column rank # y - m-dimensional real vector # p - the number of optimal solutions # # Output arguments: # x_hat - k by p real matrix # z_hat - n by p integer matrix (in double precision). # The pair {x_hat(:,j),z_hat(:,j)} is the j-th optimal solution # i.e., its residual is the j-th smallest, so # ||y-A*x_hat(:,1)-B*z_hat(:,1)||<=...<=||y-A*x_hat(:,p)-B*z_hat(:,p)|| m, k = A.shape m2, n = B.shape if m != m2 or m != len(y) or len(y[1]) != 1: raise ValueError("Input arguments have a matrix dimension error!") if rank(A) + rank(B) < k + n: raise ValueError("hmmm...") Q, R = qr(A, mode='complete') Q_A = Q[:, :k] Q_Abar = Q[:, k:] R_A = R[:k, :] # Compute the p optimal integer least squares solutions z_hat = ils(dot(Q_Abar.T, B), dot(Q_Abar.T, y), p) # Compute the corresponding real least squares solutions x_hat = lstsq(R_A, dot(Q_A.T, (dot(y, ones((1, p))) - dot(B, z_hat)))) return x_hat, z_hat
def qr_decomposition(X, job=1): """Performs the QR decomposition using LINPACK, BLAS and LAPACK Fortran subroutines. Parameters ---------- X : array_like, shape (n_samples, n_features) The matrix to decompose job : int, optional (default=1) Whether to perform pivoting. 0 is False, any other value will be coerced to 1 (True). Returns ------- X : np.ndarray, shape=(n_samples, n_features) The matrix rank : int The rank of the matrix qraux : np.ndarray, shape=(n_features,) Contains further information required to recover the orthogonal part of the decomposition. pivot : np.ndarray, shape=(n_features,) The pivot array, or None if not ``job`` """ X = check_array(X, dtype='numeric', order='F', copy=True) n, p = X.shape # check on size _validate_matrix_size(n, p) rank = matrix_rank(X) # validate job: job_ = 0 if not job else 1 qraux, pivot, work = (np.zeros(p, dtype=np.double, order='F'), # can't use arange, because need fortran order ('order' not kw in arange) np.array([i for i in range(1, p + 1)], dtype=np.int, order='F'), np.zeros(p, dtype=np.double, order='F')) # sanity checks assert qraux.shape[0] == p, 'expected qraux to be of length %i' % p assert pivot.shape[0] == p, 'expected pivot to be of length %i' % p assert work.shape[0] == p, 'expected work to be of length %i' % p # call the fortran module IN PLACE _safecall(dqrsl.dqrdc, X, n, n, p, qraux, pivot, work, job_) # do returns return (X, rank, qraux, (pivot - 1) if job_ else None) # subtract one because pivot started at 1 for the fortran
def _argrelmax(data, axis=0, order=1, mode='clip'): """Calculate the relative maxima of `data`. Parameters ---------- data : ndarray Array in which to find the relative maxima. axis : int, optional Axis over which to select from `data`. Default is 0. order : int, optional How many points on each side to use for the comparison to consider ``comparator(n, n+x)`` to be True. mode : str, optional How the edges of the vector are treated. Available options are 'wrap' (wrap around) or 'clip' (treat overflow as the same as the last (or first) element). Default 'clip'. See `numpy.take`. Returns ------- extrema : tuple of ndarrays Indices of the maxima in arrays of integers. ``extrema[k]`` is the array of indices of axis `k` of `data`. Note that the return value is a tuple even when `data` is one-dimensional. """ comparator = np.greater if((int(order) != order) or (order < 1)): raise ValueError('Order must be an int >= 1') datalen = data.shape[axis] locs = np.arange(0, datalen) results = np.ones(data.shape, dtype=bool) main = data.take(locs, axis=axis, mode=mode) for shift in xrange(1, order + 1): plus = data.take(locs + shift, axis=axis, mode=mode) minus = data.take(locs - shift, axis=axis, mode=mode) results &= comparator(main, plus) results &= comparator(main, minus) if(~results.any()): return results return np.where(results) ############################################################################### # Back porting matrix_rank for numpy < 1.7
def qrmcp(B, y): """ QR factorization of B. Input arguments: B - m by n real matrix to be factorized (np.array) y - m-dimensional real vector to be transformed to Q'y Output arguments: R - n by n real upper triangular matrix piv - n-vector storing the information of the permutation matrix P y - m-vector transformed from the input y by Q, i.e., y := Q'*y """ # Check input arguments m, n = B.shape if m < n: raise ValueError("Matrix to be factorized is column-rank deficient!") m2, n2 = y.shape if m != m2 and n2 != 1: raise ValueError("Input arguments have a dimension error!") # Initialization colnormB = zeros((2, n)) piv = arange(n) # Compute the 2-norm squared of each column of B for j in xrange(n): colnormB[0][j] = norm(B[:, j]) ** 2 for k in xrange(n): # Find the column with minimum 2-norm in B(k+1:m,k+1:n) tmp = colnormB[0, k:] - colnormB[1, k:] i = tmp.argmin() q = i + k # Column interchange if q > k: piv[k], piv[q] = piv[q].copy(), piv[k].copy() colnormB[:, k], colnormB[:, q] = colnormB[:, q].copy(), colnormB[:, k].copy() B[:, k], B[:, q] = B[:, q].copy(), B[:, k].copy() # Compute and apply the Householder transformation I-tau*v*v' if norm(B[k + 1:, k]) > 0: # otherwise no Householder transformation is needed v = (B[k:, k].copy()).reshape(m - k, 1) # v = column-vector rho = norm(v) # scalar if v[0] >= 0: rho = -rho v[0] -= rho # B(k,k)+sgn(B(k,k))*norm(B(k:n,k)) tao = -1 / (rho * v[0]) B[k, k] = rho B[k:, k + 1:] -= tao * dot(v, dot(v.T, B[k:, k + 1:])) # Update y by the Householder transformation y[k:] = y[k:, :] - tao * dot(v, dot(v.T, y[k:])) # Update colnormB(2,k+1:n) colnormB[1, k + 1:] += B[k, k + 1:] * B[k, k + 1:] return triu(B[:n, :n]), piv, y
