def random_orthogonal_matrix(n):
"""Returns a random, orthogonal matrix of n by n."""
a = np.random.randn(n, n)
U, _, _ = linalg.svd(a)
assert U.shape == (n, n)
return U
python类svd()的实例源码
def fit(self, X,
augment=False,
rounds=1,
seed=None):
'''Required for featurewise_center, featurewise_std_normalization
and zca_whitening.
# Arguments
X: Numpy array, the data to fit on.
augment: whether to fit on randomly augmented samples
rounds: if `augment`,
how many augmentation passes to do over the data
seed: random seed.
'''
if seed is not None:
np.random.seed(seed)
X = np.copy(X)
if augment:
aX = np.zeros(tuple([rounds * X.shape[0]] + list(X.shape)[1:]))
for r in range(rounds):
for i in range(X.shape[0]):
aX[i + r * X.shape[0]] = self.random_transform(X[i])
X = aX
if self.featurewise_center:
self.mean = np.mean(X, axis=0)
X -= self.mean
if self.featurewise_std_normalization:
self.std = np.std(X, axis=0)
X /= (self.std + 1e-7)
if self.zca_whitening:
flatX = np.reshape(X, (X.shape[0], X.shape[1] * X.shape[2] * X.shape[3]))
sigma = np.dot(flatX.T, flatX) / flatX.shape[0]
U, S, V = linalg.svd(sigma)
self.principal_components = np.dot(np.dot(U, np.diag(1. / np.sqrt(S + 10e-7))), U.T)
def svd_wrapper(X, rank = None):
"""
Computes the (possibly partial) SVD of a matrix. Handles the case where
X is either dense or sparse.
Parameters
----------
X: either dense or sparse
rank: rank of the desired SVD (required for sparse matrices)
Output
------
U, D, V
the columns of U are the left singular vectors
the COLUMNS of V are the left singular vectors
"""
if isinstance(X, LinearOperator):
scipy_svds = svds(convert2scipy(X), rank)
U, D, V = fix_scipy_svds(scipy_svds)
V = V.T
elif issparse(X):
scipy_svds = svds(X, rank)
U, D, V = fix_scipy_svds(scipy_svds)
V = V.T
else:
# TODO: implement partial SVD
U, D, V = full_svd(X, full_matrices=False)
V = V.T
if rank:
U = U[:, :rank]
D = D[:rank]
V = V[:, :rank]
return U, D, V
def low_rank_repr(X, n_dim):
U, S, V = svd(X.T,full_matrices=False)
mask = S > 1
V = V[mask]
S = S[mask]
R = (V.T * (1 - S**-2)).dot(V)
return R
def _hrchy(DB, classes, src_idx, trg_idx):
from hierarchy import load
def simplex(num_vecs):
from scipy.linalg import svd
cw = np.eye(num_vecs)
cw -= cw.mean(axis=1, keepdims=True)
cw, _, _ = svd(cw)
return cw[:, :num_vecs-1]
hrchy_fn = 'data/%s/wordnet_hierarchy.txt' % DB
hrchy = load(hrchy_fn)
constrains = OrderedDict()
k = 0
for q in hrchy.nodes:
if len(q.fpointer) <= 1:
continue
labs = len(q.fpointer)*np.ones((len(classes)))
for s, cls in enumerate(classes):
if cls in q.name['classes']:
labs[s] = [cls in hrchy[cid].name['classes'] for cid in q.fpointer].index(True)
assigns = np.zeros((len(q.fpointer)+1, len(classes)))
for c in range(len(q.fpointer)+1):
assigns[c, labs == c] = 1
constrains[q.identifier] = {'codewords': simplex(len(q.fpointer)+1).T,
'idxDim': np.arange(k, k+len(q.fpointer)),
'labels': {'source': assigns[:, src_idx],
'target': assigns[:, trg_idx]}}
k += len(q.fpointer)
source_deg = [key for key in constrains.keys() if constrains[key]['labels']['source'].sum(axis=1).max() == len(src_idx)]
[constrains.pop(key) for key in source_deg]
return constrains
def load_caltech101_30(folder=CALTECH101_30_DIR, tiny_problem=False):
caltech = scio.loadmat(folder + '/caltech101-30.matlab')
k_train, k_test = caltech['Ktrain'], caltech['Ktest']
label_tr, label_te = caltech['tr_label'], caltech['te_label']
file_tr, file_te = caltech['tr_files'], caltech['te_files']
if tiny_problem:
pattern_step = 5
fraction_limit = 0.2
k_train = k_train[:int(len(label_tr) * fraction_limit):pattern_step,
:int(len(label_tr) * fraction_limit):pattern_step]
label_tr = label_tr[:int(len(label_tr) * fraction_limit):pattern_step]
U, s, Vh = linalg.svd(k_train)
S_sqrt = linalg.diagsvd(s ** 0.5, len(s), len(s))
X = np.dot(U, S_sqrt) # examples in rows
train_x, val_x, test_x = X[0:len(X):3, :], X[1:len(X):3, :], X[2:len(X):3, :]
label_tr_enc = to_one_hot_enc(np.array(label_tr) - 1)
train_y, val_y, test_y = label_tr_enc[0:len(X):3, :], label_tr_enc[1:len(X):3, :], label_tr_enc[2:len(X):3, :]
train_file, val_file, test_file = file_tr[0:len(X):3], file_tr[1:len(X):3], file_tr[2:len(X):3]
test_dataset = Dataset(data=test_x, target=test_y, info={'files': test_file})
validation_dataset = Dataset(data=val_x, target=val_y, info={'files': val_file})
training_dataset = Dataset(data=train_x, target=train_y, info={'files': train_file})
return Datasets(train=training_dataset, validation=validation_dataset, test=test_dataset)
def test():
data = np.random.random((5000, 100))
u, s, v = linalg.svd(data, full_matrices=False)
pca = np.dot(u[:, :10].T, data)
results = fastica(pca.T, whiten=False)
def test():
data = np.random.random((5000, 100))
u, s, v = linalg.svd(data)
pca = np.dot(u[:, :10].T, data)
results = fastica(pca.T, whiten=False)
def main():
parser = argparse.ArgumentParser('Transforms the space by applying the SVD')
parser.add_argument('--input', '-i', help='Input file')
parser.add_argument('--output', '-o', help='Output file')
args = parser.parse_args()
space = load_mikolov_text(args.input)
U, s, Vh = svd(space.matrix, full_matrices=False)
transformed = U.dot(np.diag(s))
newspace = VectorSpace(transformed, space.vocab)
newspace.save_mikolov_text(args.output)
def fit(self, x):
s = x.shape
x = x.copy().reshape((s[0],np.prod(s[1:])))
m = np.mean(x, axis=0)
x -= m
sigma = np.dot(x.T,x) / x.shape[0]
U, S, V = linalg.svd(sigma)
tmp = np.dot(U, np.diag(1./np.sqrt(S+self.regularization)))
tmp2 = np.dot(U, np.diag(np.sqrt(S+self.regularization)))
self.ZCA_mat = th.shared(np.dot(tmp, U.T).astype(th.config.floatX))
self.inv_ZCA_mat = th.shared(np.dot(tmp2, U.T).astype(th.config.floatX))
self.mean = th.shared(m.astype(th.config.floatX))
def fit(self, X,
augment=False, # fit on randomly augmented samples
rounds=1, # if augment, how many augmentation passes over the data do we use
seed=None
):
'''
Required for featurewise_center, featurewise_std_normalization and zca_whitening.
'''
X = np.copy(X)
if augment:
aX = np.zeros(tuple([rounds*X.shape[0]]+list(X.shape)[1:]))
for r in range(rounds):
for i in range(X.shape[0]):
img = array_to_img(X[i])
img = self.random_transform(img)
aX[i+r*X.shape[0]] = img_to_array(img)
X = aX
if self.featurewise_center:
self.mean = np.mean(X, axis=0)
X -= self.mean
if self.featurewise_std_normalization:
self.std = np.std(X, axis=0)
X /= self.std
if self.zca_whitening:
flatX = np.reshape(X, (X.shape[0], X.shape[1]*X.shape[2]*X.shape[3]))
fudge = 10e-6
sigma = np.dot(flatX.T, flatX) / flatX.shape[1]
U, S, V = linalg.svd(sigma)
self.principal_components = np.dot(np.dot(U, np.diag(1. / np.sqrt(S + fudge))), U.T)
def condensed_svd(self, M, tol=1e-3, store_singular_vals=False):
U, S, Vt = linalg.svd(M, full_matrices=False)
if store_singular_vals:
self.singular_vals = S
#want tolerance on fraction of variance in singular value
#when not norm_covariance, need to normalize singular values
S_norm = 1.0 if self.norm_covariance else np.sum(S)
rank = np.sum( (S/S_norm) > tol )
return U[:,:rank], S[:rank], Vt[:rank,:], S_norm
def condensed_svd(self, M, tol=1e-3, store_singular_vals=False):
U, S, Vt = linalg.svd(M, full_matrices=False)
if store_singular_vals:
self.singular_vals = S
#want tolerance on fraction of variance in singular value
#when not norm_covariance, need to normalize singular values
S_norm = np.sum(S)
rank = np.sum( (S/S_norm) > tol )
return U[:,:rank], S[:rank], Vt[:rank,:], S_norm
def make_orthogonal(self):
# orthogonalize the columns
U, s, V = svd(self.Tm, full_matrices=False)
self.Tm = np.diag(s).dot(V)
def ComputeCCA(X, Y):
assert X.shape[0] == Y.shape[0], (X.shape, Y.shape, "Unequal number of rows")
assert X.shape[0] > 1, (X.shape, "Must have more than 1 row")
X = NormCenterMatrix(X)
Y = NormCenterMatrix(Y)
X_q, _, _ = decomp_qr.qr(X, overwrite_a=True, mode='economic', pivoting=True)
Y_q, _, _ = decomp_qr.qr(Y, overwrite_a=True, mode='economic', pivoting=True)
C = np.dot(X_q.T, Y_q)
r = linalg.svd(C, full_matrices=False, compute_uv=False)
d = min(X.shape[1], Y.shape[1])
r = r[:d]
r = np.minimum(np.maximum(r, 0.0), 1.0) # remove roundoff errs
return r.mean()
def pseudo_inverse(mat, eps=1e-10):
"""Computes pseudo-inverse of mat, treating eigenvalues below eps as 0."""
s, u, v = tf.svd(mat)
eps = 1e-10 # zero threshold for eigenvalues
si = tf.where(tf.less(s, eps), s, 1./s)
return u @ tf.diag(si) @ tf.transpose(v)
def symsqrt(mat, eps=1e-7):
"""Symmetric square root."""
s, u, v = tf.svd(mat)
# sqrt is unstable around 0, just use 0 in such case
print("Warning, cutting off at eps")
si = tf.where(tf.less(s, eps), s, tf.sqrt(s))
return u @ tf.diag(si) @ tf.transpose(v)
def pseudo_inverse_sqrt(mat, eps=1e-7):
"""half pseduo-inverse"""
s, u, v = tf.svd(mat)
# zero threshold for eigenvalues
si = tf.where(tf.less(s, eps), s, 1./tf.sqrt(s))
return u @ tf.diag(si) @ tf.transpose(v)
def pseudo_inverse_sqrt2(svd, eps=1e-7):
"""half pseduo-inverse, accepting existing values"""
# zero threshold for eigenvalues
if svd.__class__.__name__=='SvdTuple':
(s, u, v) = (svd.s, svd.u, svd.v)
elif svd.__class__.__name__=='SvdWrapper':
(s, u, v) = (svd.s, svd.u, svd.v)
else:
assert False, "Unknown type"
si = tf.where(tf.less(s, eps), s, 1./tf.sqrt(s))
return u @ tf.diag(si) @ tf.transpose(v)
def pseudo_inverse2(svd, eps=1e-7):
"""pseudo-inverse, accepting existing values"""
# use float32 machine precision as cut-off (works for MKL)
# https://www.wolframcloud.com/objects/927b2aa5-de9c-46f5-89fe-c4a58aa4c04b
if svd.__class__.__name__=='SvdTuple':
(s, u, v) = (svd.s, svd.u, svd.v)
elif svd.__class__.__name__=='SvdWrapper':
(s, u, v) = (svd.s, svd.u, svd.v)
else:
assert False, "Unknown type"
max_eigen = tf.reduce_max(s)
si = tf.where(s/max_eigen<eps, 0.*s, 1./s)
return u @ tf.diag(si) @ tf.transpose(v)
