def do(self, a, b):
# note that eigenvalue arrays returned by eig must be sorted since
# their order isn't guaranteed.
ev, evc = linalg.eigh(a)
evalues, evectors = linalg.eig(a)
evalues.sort(axis=-1)
assert_almost_equal(ev, evalues)
assert_allclose(dot_generalized(a, evc),
np.asarray(ev)[..., None, :] * np.asarray(evc),
rtol=get_rtol(ev.dtype))
ev2, evc2 = linalg.eigh(a, 'U')
assert_almost_equal(ev2, evalues)
assert_allclose(dot_generalized(a, evc2),
np.asarray(ev2)[..., None, :] * np.asarray(evc2),
rtol=get_rtol(ev.dtype), err_msg=repr(a))
python类eigh()的实例源码
def test_UPLO(self):
Klo = np.array([[0, 0], [1, 0]], dtype=np.double)
Kup = np.array([[0, 1], [0, 0]], dtype=np.double)
tgt = np.array([-1, 1], dtype=np.double)
rtol = get_rtol(np.double)
# Check default is 'L'
w, v = np.linalg.eigh(Klo)
assert_allclose(w, tgt, rtol=rtol)
# Check 'L'
w, v = np.linalg.eigh(Klo, UPLO='L')
assert_allclose(w, tgt, rtol=rtol)
# Check 'l'
w, v = np.linalg.eigh(Klo, UPLO='l')
assert_allclose(w, tgt, rtol=rtol)
# Check 'U'
w, v = np.linalg.eigh(Kup, UPLO='U')
assert_allclose(w, tgt, rtol=rtol)
# Check 'u'
w, v = np.linalg.eigh(Kup, UPLO='u')
assert_allclose(w, tgt, rtol=rtol)
def GetBestFitPlane(pts, weights=None):
if weights is None:
wSum = len(pts)
origin = np.sum(pts, 0)
origin /= wSum
sums = np.zeros((3, 3), np.double)
for pt in pts:
dp = pt - origin
for i in range(3):
sums[i, i] += dp[i] * dp[i]
for j in range(i + 1, 3):
sums[i, j] += dp[i] * dp[j]
sums[j, i] += dp[i] * dp[j]
sums /= wSum
vals, vects = linalg.eigh(sums)
order = np.argsort(vals)
normal = vects[:, order[0]]
plane = np.zeros((4, ), np.double)
plane[:3] = normal
plane[3] = -1 * normal.dot(origin)
return plane
def do(self, a, b):
# note that eigenvalue arrays returned by eig must be sorted since
# their order isn't guaranteed.
ev, evc = linalg.eigh(a)
evalues, evectors = linalg.eig(a)
evalues.sort(axis=-1)
assert_almost_equal(ev, evalues)
assert_allclose(dot_generalized(a, evc),
np.asarray(ev)[..., None, :] * np.asarray(evc),
rtol=get_rtol(ev.dtype))
ev2, evc2 = linalg.eigh(a, 'U')
assert_almost_equal(ev2, evalues)
assert_allclose(dot_generalized(a, evc2),
np.asarray(ev2)[..., None, :] * np.asarray(evc2),
rtol=get_rtol(ev.dtype), err_msg=repr(a))
def test_UPLO(self):
Klo = np.array([[0, 0], [1, 0]], dtype=np.double)
Kup = np.array([[0, 1], [0, 0]], dtype=np.double)
tgt = np.array([-1, 1], dtype=np.double)
rtol = get_rtol(np.double)
# Check default is 'L'
w, v = np.linalg.eigh(Klo)
assert_allclose(w, tgt, rtol=rtol)
# Check 'L'
w, v = np.linalg.eigh(Klo, UPLO='L')
assert_allclose(w, tgt, rtol=rtol)
# Check 'l'
w, v = np.linalg.eigh(Klo, UPLO='l')
assert_allclose(w, tgt, rtol=rtol)
# Check 'U'
w, v = np.linalg.eigh(Kup, UPLO='U')
assert_allclose(w, tgt, rtol=rtol)
# Check 'u'
w, v = np.linalg.eigh(Kup, UPLO='u')
assert_allclose(w, tgt, rtol=rtol)
test_linalg.py 文件源码
项目:PyDataLondon29-EmbarrassinglyParallelDAWithAWSLambda
作者: SignalMedia
项目源码
文件源码
阅读 27
收藏 0
点赞 0
评论 0
def do(self, a, b):
# note that eigenvalue arrays returned by eig must be sorted since
# their order isn't guaranteed.
ev, evc = linalg.eigh(a)
evalues, evectors = linalg.eig(a)
evalues.sort(axis=-1)
assert_almost_equal(ev, evalues)
assert_allclose(dot_generalized(a, evc),
np.asarray(ev)[..., None, :] * np.asarray(evc),
rtol=get_rtol(ev.dtype))
ev2, evc2 = linalg.eigh(a, 'U')
assert_almost_equal(ev2, evalues)
assert_allclose(dot_generalized(a, evc2),
np.asarray(ev2)[..., None, :] * np.asarray(evc2),
rtol=get_rtol(ev.dtype), err_msg=repr(a))
test_linalg.py 文件源码
项目:PyDataLondon29-EmbarrassinglyParallelDAWithAWSLambda
作者: SignalMedia
项目源码
文件源码
阅读 26
收藏 0
点赞 0
评论 0
def test_UPLO(self):
Klo = np.array([[0, 0], [1, 0]], dtype=np.double)
Kup = np.array([[0, 1], [0, 0]], dtype=np.double)
tgt = np.array([-1, 1], dtype=np.double)
rtol = get_rtol(np.double)
# Check default is 'L'
w, v = np.linalg.eigh(Klo)
assert_allclose(w, tgt, rtol=rtol)
# Check 'L'
w, v = np.linalg.eigh(Klo, UPLO='L')
assert_allclose(w, tgt, rtol=rtol)
# Check 'l'
w, v = np.linalg.eigh(Klo, UPLO='l')
assert_allclose(w, tgt, rtol=rtol)
# Check 'U'
w, v = np.linalg.eigh(Kup, UPLO='U')
assert_allclose(w, tgt, rtol=rtol)
# Check 'u'
w, v = np.linalg.eigh(Kup, UPLO='u')
assert_allclose(w, tgt, rtol=rtol)
def __newlambdagammaC(self, theta, l):
""" Apply Eqns. 25-27 in Vidal 2003 to update lambda^C and gamma^C
(lambda and gamma of this qbit).
"""
gamma_ket = self.coefs[l+1].lam
gamma_bra = np.conjugate(gamma_ket)
Gamma_star = np.conjugate(self.coefs[l+1].gamma)
inputs = [Gamma_star, theta, gamma_bra, gamma_ket]
Gamma_star_idx = [1, -3, -2]
theta_idx = [-1, 1, -4, -5]
gamma_bra_idx = [-6]
gamma_ket_idx = [-7]
idx = [Gamma_star_idx, theta_idx, gamma_bra_idx, gamma_ket_idx]
contract_me = scon(inputs, idx)
svd_me = np.einsum('agibggg', contract_me)
evals, evecs = la.eigh(svd_me)
return evals, evecs
def do(self, a, b):
# note that eigenvalue arrays returned by eig must be sorted since
# their order isn't guaranteed.
ev, evc = linalg.eigh(a)
evalues, evectors = linalg.eig(a)
evalues.sort(axis=-1)
assert_almost_equal(ev, evalues)
assert_allclose(dot_generalized(a, evc),
np.asarray(ev)[..., None, :] * np.asarray(evc),
rtol=get_rtol(ev.dtype))
ev2, evc2 = linalg.eigh(a, 'U')
assert_almost_equal(ev2, evalues)
assert_allclose(dot_generalized(a, evc2),
np.asarray(ev2)[..., None, :] * np.asarray(evc2),
rtol=get_rtol(ev.dtype), err_msg=repr(a))
def test_UPLO(self):
Klo = np.array([[0, 0], [1, 0]], dtype=np.double)
Kup = np.array([[0, 1], [0, 0]], dtype=np.double)
tgt = np.array([-1, 1], dtype=np.double)
rtol = get_rtol(np.double)
# Check default is 'L'
w, v = np.linalg.eigh(Klo)
assert_allclose(w, tgt, rtol=rtol)
# Check 'L'
w, v = np.linalg.eigh(Klo, UPLO='L')
assert_allclose(w, tgt, rtol=rtol)
# Check 'l'
w, v = np.linalg.eigh(Klo, UPLO='l')
assert_allclose(w, tgt, rtol=rtol)
# Check 'U'
w, v = np.linalg.eigh(Kup, UPLO='U')
assert_allclose(w, tgt, rtol=rtol)
# Check 'u'
w, v = np.linalg.eigh(Kup, UPLO='u')
assert_allclose(w, tgt, rtol=rtol)
def plot_ellipses(self):
q = 0.60
r = ops.get_ellipse_rad(q)
H = ops.get_hessian(self.X)
eigv, rotation = la.eigh(H)
center = self.Bh
u = np.linspace(0.0, 2.0 * np.pi, 100)
v = np.linspace(0.0, np.pi, 100)
x = (r/math.sqrt(eigv[0]))* np.outer(np.cos(u), np.sin(v))
y = (r/math.sqrt(eigv[1])) * np.outer(np.sin(u), np.sin(v))
z = (r/math.sqrt(eigv[2])) * np.outer(np.ones_like(u), np.cos(v))
for i in range(len(x)):
for j in range(len(x)):
[x[i,j],y[i,j],z[i,j]] = np.dot([x[i,j],y[i,j],z[i,j]], rotation) + center
plt.plot(z,x)
plt.axis('equal')
plt.show()
def spectral_embedding(mat, target_dim, gramian=True, discard_first=True):
if discard_first:
last = -1
first = target_dim - last
else:
first = target_dim
last = None
if not gramian:
mat = mat.dot(mat.T)
eigvals, eigvecs = eigh(mat)
sl = slice(-first, last)
eigvecs = eigvecs[:, sl]
eigvals_crop = eigvals[sl]
Y = eigvecs.dot(np.diag(np.sqrt(eigvals_crop)))
Y = Y[:, ::-1]
variance_explaned(eigvals, eigvals_crop)
return Y
def do(self, a, b):
# note that eigenvalue arrays returned by eig must be sorted since
# their order isn't guaranteed.
ev, evc = linalg.eigh(a)
evalues, evectors = linalg.eig(a)
evalues.sort(axis=-1)
assert_almost_equal(ev, evalues)
assert_allclose(dot_generalized(a, evc),
np.asarray(ev)[..., None, :] * np.asarray(evc),
rtol=get_rtol(ev.dtype))
ev2, evc2 = linalg.eigh(a, 'U')
assert_almost_equal(ev2, evalues)
assert_allclose(dot_generalized(a, evc2),
np.asarray(ev2)[..., None, :] * np.asarray(evc2),
rtol=get_rtol(ev.dtype), err_msg=repr(a))
def test_UPLO(self):
Klo = np.array([[0, 0], [1, 0]], dtype=np.double)
Kup = np.array([[0, 1], [0, 0]], dtype=np.double)
tgt = np.array([-1, 1], dtype=np.double)
rtol = get_rtol(np.double)
# Check default is 'L'
w, v = np.linalg.eigh(Klo)
assert_allclose(w, tgt, rtol=rtol)
# Check 'L'
w, v = np.linalg.eigh(Klo, UPLO='L')
assert_allclose(w, tgt, rtol=rtol)
# Check 'l'
w, v = np.linalg.eigh(Klo, UPLO='l')
assert_allclose(w, tgt, rtol=rtol)
# Check 'U'
w, v = np.linalg.eigh(Kup, UPLO='U')
assert_allclose(w, tgt, rtol=rtol)
# Check 'u'
w, v = np.linalg.eigh(Kup, UPLO='u')
assert_allclose(w, tgt, rtol=rtol)
def make_eigvals_positive(am, targetprod):
"""For the symmetric square matrix `am`, increase any zero eigenvalues
such that the total product of eigenvalues is greater or equal to
`targetprod`. Returns a (possibly) new, non-singular matrix."""
w, v = linalg.eigh(am) # use eigh since a is symmetric
mask = w < 1.e-10
if np.any(mask):
nzprod = np.product(w[~mask]) # product of nonzero eigenvalues
nzeros = mask.sum() # number of zero eigenvalues
new_val = max(1.e-10, (targetprod / nzprod) ** (1. / nzeros))
w[mask] = new_val # adjust zero eigvals
am_new = np.dot(np.dot(v, np.diag(w)), linalg.inv(v)) # re-form cov
else:
am_new = am
return am_new
def do(self, a, b):
# note that eigenvalue arrays returned by eig must be sorted since
# their order isn't guaranteed.
ev, evc = linalg.eigh(a)
evalues, evectors = linalg.eig(a)
evalues.sort(axis=-1)
assert_almost_equal(ev, evalues)
assert_allclose(dot_generalized(a, evc),
np.asarray(ev)[..., None, :] * np.asarray(evc),
rtol=get_rtol(ev.dtype))
ev2, evc2 = linalg.eigh(a, 'U')
assert_almost_equal(ev2, evalues)
assert_allclose(dot_generalized(a, evc2),
np.asarray(ev2)[..., None, :] * np.asarray(evc2),
rtol=get_rtol(ev.dtype), err_msg=repr(a))
def test_UPLO(self):
Klo = np.array([[0, 0], [1, 0]], dtype=np.double)
Kup = np.array([[0, 1], [0, 0]], dtype=np.double)
tgt = np.array([-1, 1], dtype=np.double)
rtol = get_rtol(np.double)
# Check default is 'L'
w, v = np.linalg.eigh(Klo)
assert_allclose(w, tgt, rtol=rtol)
# Check 'L'
w, v = np.linalg.eigh(Klo, UPLO='L')
assert_allclose(w, tgt, rtol=rtol)
# Check 'l'
w, v = np.linalg.eigh(Klo, UPLO='l')
assert_allclose(w, tgt, rtol=rtol)
# Check 'U'
w, v = np.linalg.eigh(Kup, UPLO='U')
assert_allclose(w, tgt, rtol=rtol)
# Check 'u'
w, v = np.linalg.eigh(Kup, UPLO='u')
assert_allclose(w, tgt, rtol=rtol)
def inv_sqrth(x):
"""
Matrix inverse square root
Parameters
----------
x : ndarray
Real, symmetric matrix
Returns
-------
invsqrt : ndarray
Input to the power -1/2
"""
vals, vecs = eigh(x)
return vecs @ diag(1 / sqrt(vals)) @ vecs.T
def do(self, a, b):
# note that eigenvalue arrays returned by eig must be sorted since
# their order isn't guaranteed.
ev, evc = linalg.eigh(a)
evalues, evectors = linalg.eig(a)
evalues.sort(axis=-1)
assert_almost_equal(ev, evalues)
assert_allclose(dot_generalized(a, evc),
np.asarray(ev)[..., None, :] * np.asarray(evc),
rtol=get_rtol(ev.dtype))
ev2, evc2 = linalg.eigh(a, 'U')
assert_almost_equal(ev2, evalues)
assert_allclose(dot_generalized(a, evc2),
np.asarray(ev2)[..., None, :] * np.asarray(evc2),
rtol=get_rtol(ev.dtype), err_msg=repr(a))
def test_UPLO(self):
Klo = np.array([[0, 0], [1, 0]], dtype=np.double)
Kup = np.array([[0, 1], [0, 0]], dtype=np.double)
tgt = np.array([-1, 1], dtype=np.double)
rtol = get_rtol(np.double)
# Check default is 'L'
w, v = np.linalg.eigh(Klo)
assert_allclose(w, tgt, rtol=rtol)
# Check 'L'
w, v = np.linalg.eigh(Klo, UPLO='L')
assert_allclose(w, tgt, rtol=rtol)
# Check 'l'
w, v = np.linalg.eigh(Klo, UPLO='l')
assert_allclose(w, tgt, rtol=rtol)
# Check 'U'
w, v = np.linalg.eigh(Kup, UPLO='U')
assert_allclose(w, tgt, rtol=rtol)
# Check 'u'
w, v = np.linalg.eigh(Kup, UPLO='u')
assert_allclose(w, tgt, rtol=rtol)
def get_whitening_matrix(X, fudge=1E-18):
from numpy.linalg import eigh
Xcov = numpy.dot(X.T, X)/X.shape[0]
d,V = eigh(Xcov)
D = numpy.diag(1./numpy.sqrt(d+fudge))
W = numpy.dot(numpy.dot(V,D), V.T)
return W
def __init__(self, P, q, r, relop=None):
self.P, self.q, self.r = P, q, r
self.qarray = np.squeeze(np.asarray(q.todense()))
self.relop = relop
self.eigh = None # for ADMM
# Evalutes f with a numpy array x.
def dc_split(self, use_eigen_split=False):
n = self.P.shape[0]
if self.P.nnz == 0: # P is zero
P1, P2 = sp.csr_matrix((n, n)), sp.csr_matrix((n, n))
if use_eigen_split:
lmb, Q = LA.eigh(self.P.todense())
P1 = sum([Q[:, i]*lmb[i]*Q[:, i].T for i in range(n) if lmb[i] > 0])
P2 = sum([-Q[:, i]*lmb[i]*Q[:, i].T for i in range(n) if lmb[i] < 0])
assert abs(np.sum(P1 - P2 - self.P)) < 1e-8
else:
lmb_min = np.min(LA.eigh(self.P.todense())[0])
if lmb_min < 0:
P1 = self.P + (1-lmb_min)*sp.identity(n)
P2 = (1-lmb_min)*sp.identity(n)
else:
P1 = self.P
P2 = sp.csr_matrix((n, n))
f1 = QuadraticFunction(P1, self.q, self.r)
f2 = QuadraticFunction(P2, sp.csc_matrix((n, 1)), 0)
return (f1, f2)
# Returns the one-variable function when regarding f(x)
# as a quadratic expression in x[k].
# f is an instance of QuadraticFunction
# return value is an instance of OneVarQuadraticFunction
# TODO: speedup
def improve_admm(x0, prob, *args, **kwargs):
num_iters = kwargs.get('num_iters', 1000)
viol_lim = kwargs.get('viol_lim', 1e4)
tol = kwargs.get('tol', 1e-2)
rho = kwargs.get('rho', None)
phase1 = kwargs.get('phase1', True)
if rho is not None:
lmb0, P0Q = map(np.asmatrix, LA.eigh(prob.f0.P.todense()))
lmb_min = np.min(lmb0)
if lmb_min + prob.m*rho < 0:
logging.error("rho parameter is too small, z-update not convex.")
logging.error("Minimum possible value of rho: %.3f\n", -lmb_min/prob.m)
logging.error("Given value of rho: %.3f\n", rho)
raise Exception("rho parameter is too small, need at least %.3f." % rho)
# TODO: find a reasonable auto parameter
if rho is None:
lmb0, P0Q = map(np.asmatrix, LA.eigh(prob.f0.P.todense()))
lmb_min = np.min(lmb0)
lmb_max = np.max(lmb0)
if lmb_min < 0: rho = 2.*(1.-lmb_min)/prob.m
else: rho = 1./prob.m
rho *= 50.
logging.warning("Automatically setting rho to %.3f", rho)
if phase1:
x1 = prob.better(x0, admm_phase1(x0, prob, tol, num_iters))
else:
x1 = x0
x2 = prob.better(x1, admm_phase2(x1, prob, rho, tol, num_iters, viol_lim))
return x2
def test_eigh_build(self, level=rlevel):
# Ticket 662.
rvals = [68.60568999, 89.57756725, 106.67185574]
cov = array([[77.70273908, 3.51489954, 15.64602427],
[3.51489954, 88.97013878, -1.07431931],
[15.64602427, -1.07431931, 98.18223512]])
vals, vecs = linalg.eigh(cov)
assert_array_almost_equal(vals, rvals)
def test_types(self):
def check(dtype):
x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
w, v = np.linalg.eigh(x)
assert_equal(w.dtype, get_real_dtype(dtype))
assert_equal(v.dtype, dtype)
for dtype in [single, double, csingle, cdouble]:
yield check, dtype
def test_invalid(self):
x = np.array([[1, 0.5], [0.5, 1]], dtype=np.float32)
assert_raises(ValueError, np.linalg.eigh, x, UPLO="lrong")
assert_raises(ValueError, np.linalg.eigh, x, "lower")
assert_raises(ValueError, np.linalg.eigh, x, "upper")
def getPCAVideo(I):
ICov = I.dot(I.T)
[lam, V] = linalg.eigh(ICov)
lam[lam < 0] = 0
V = V*np.sqrt(lam[None, :])
return V
def average_structure(X):
"""
Calculate an average structure from an ensemble of structures
(i.e. X is a rank-3 tensor: X[i] is a (N,3) configuration matrix).
@param X: m x n x 3 input vector
@type X: numpy array
@return: average structure
@rtype: (n,3) numpy.array
"""
from numpy.linalg import eigh
B = csb.numeric.gower_matrix(X)
v, U = eigh(B)
if numpy.iscomplex(v).any():
v = v.real
if numpy.iscomplex(U).any():
U = U.real
indices = numpy.argsort(v)[-3:]
v = numpy.take(v, indices, 0)
U = numpy.take(U, indices, 1)
x = U * numpy.sqrt(v)
i = 0
while is_mirror_image(x, X[0]) and i < 2:
x[:, i] *= -1
i += 1
return x
def test_eigh_build(self, level=rlevel):
# Ticket 662.
rvals = [68.60568999, 89.57756725, 106.67185574]
cov = array([[77.70273908, 3.51489954, 15.64602427],
[3.51489954, 88.97013878, -1.07431931],
[15.64602427, -1.07431931, 98.18223512]])
vals, vecs = linalg.eigh(cov)
assert_array_almost_equal(vals, rvals)