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
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
test_linalg.py 文件源码
项目:PyDataLondon29-EmbarrassinglyParallelDAWithAWSLambda
作者: SignalMedia
项目源码
文件源码
阅读 24
收藏 0
点赞 0
评论 0
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)
test_linalg.py 文件源码
项目:PyDataLondon29-EmbarrassinglyParallelDAWithAWSLambda
作者: SignalMedia
项目源码
文件源码
阅读 46
收藏 0
点赞 0
评论 0
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]