我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用numpy.linalg()。
def logdet(C,eps=1e-6,safe=0): ''' Logarithm of the determinant of a matrix Works only with real-valued positive definite matrices ''' try: return 2.0*np.sum(np.log(np.diag(chol(C)))) except numpy.linalg.linalg.LinAlgError: if safe: C = check_covmat(C,eps=eps) w = np.linalg.eigh(C)[0] w = np.real(w) w[w<eps]=eps det = np.sum(np.log(w)) return det
def _get_skew(corners, board): """ Get skew for given checkerboard detection. Scaled to [0,1], which 0 = no skew, 1 = high skew Skew is proportional to the divergence of three outside corners from 90 degrees. """ # TODO Using three nearby interior corners might be more robust, outside corners occasionally # get mis-detected up_left, up_right, down_right, _ = _get_outside_corners(corners, board) def angle(a, b, c): """ Return angle between lines ab, bc """ ab = a - b cb = c - b return math.acos(numpy.dot(ab,cb) / (numpy.linalg.norm(ab) * numpy.linalg.norm(cb))) skew = min(1.0, 2. * abs((math.pi / 2.) - angle(up_left, up_right, down_right))) return skew
def get_ground_state(sparse_operator): """Compute lowest eigenvalue and eigenstate. Returns: eigenvalue: The lowest eigenvalue, a float. eigenstate: The lowest eigenstate in scipy.sparse csc format. """ if not is_hermitian(sparse_operator): raise ValueError('sparse_operator must be Hermitian.') values, vectors = scipy.sparse.linalg.eigsh( sparse_operator, 2, which='SA', maxiter=1e7) eigenstate = scipy.sparse.csc_matrix(vectors[:, 0]) eigenvalue = values[0] return eigenvalue, eigenstate.getH()
def __init__(self, mesh, transform_out=None, transform_in=None): transform_out = numpy.matrix(transform_out) if transform_out is not None else None transform_in = numpy.matrix(transform_in) if transform_in is not None else None if transform_in is None and transform_out is None: transform_in = numpy.identity(3) transform_out = numpy.identity(3) elif transform_in is None: try: transform_in = numpy.linalg.inv(transform_out) except: transform_in = None elif transform_out is None: try: transform_out = numpy.linalg.inv(transform_in) except: transform_out = None self.transform_out, self.transform_in = transform_out, transform_in super().__init__( mesh, warp_in=lambda vertex: self.transform_in.dot(vertex).tolist()[0][:3] if self.transform_in else None, warp_out=lambda vertex: self.transform_out.dot(vertex).tolist()[0][:3] if self.transform_out else None, )
def likelihood(x, m=None, Cinv=None, sigma=1, detC=None): """return likelihood of x for the normal density N(m, sigma**2 * Cinv**-1)""" # testing: MC integrate must be one: mean(p(x_i)) * volume(where x_i are uniformely sampled) # for i in xrange(3): print mean([cma.likelihood(20*r-10, dim * [0], # None, 3) for r in rand(10000,dim)]) * 20**dim if m is None: dx = x else: dx = x - m # array(x) - array(m) n = len(x) s2pi = (2 * np.pi)**(n / 2.) if Cinv is None: return exp(-sum(dx**2) / sigma**2 / 2) / s2pi / sigma**n if detC is None: detC = 1. / np.linalg.linalg.det(Cinv) return exp(-np.dot(dx, np.dot(Cinv, dx)) / sigma**2 / 2) / s2pi / abs(detC)**0.5 / sigma**n
def test_pseudoinverse_correctness(): rng = numpy.random.RandomState(utt.fetch_seed()) d1 = rng.randint(4) + 2 d2 = rng.randint(4) + 2 r = rng.randn(d1, d2).astype(theano.config.floatX) x = tensor.matrix() xi = pinv(x) ri = function([x], xi)(r) assert ri.shape[0] == r.shape[1] assert ri.shape[1] == r.shape[0] assert ri.dtype == r.dtype # Note that pseudoinverse can be quite unprecise so I prefer to compare # the result with what numpy.linalg returns assert _allclose(ri, numpy.linalg.pinv(r))
def test_numpy_compare(self): rng = numpy.random.RandomState(utt.fetch_seed()) M = tensor.matrix("A", dtype=theano.config.floatX) V = tensor.vector("V", dtype=theano.config.floatX) a = rng.rand(4, 4).astype(theano.config.floatX) b = rng.rand(4).astype(theano.config.floatX) A = ( [None, 'fro', 'inf', '-inf', 1, -1, None, 'inf', '-inf', 0, 1, -1, 2, -2], [M, M, M, M, M, M, V, V, V, V, V, V, V, V], [a, a, a, a, a, a, b, b, b, b, b, b, b, b], [None, 'fro', inf, -inf, 1, -1, None, inf, -inf, 0, 1, -1, 2, -2]) for i in range(0, 14): f = function([A[1][i]], norm(A[1][i], A[0][i])) t_n = f(A[2][i]) n_n = numpy.linalg.norm(A[2][i], A[3][i]) assert _allclose(n_n, t_n)
def test_eval(self): A = self.A Ai = tensorinv(A) n_ainv = numpy.linalg.tensorinv(self.a) tf_a = function([A], [Ai]) t_ainv = tf_a(self.a) assert _allclose(n_ainv, t_ainv) B = self.B Bi = tensorinv(B) Bi1 = tensorinv(B, ind=1) n_binv = numpy.linalg.tensorinv(self.b) n_binv1 = numpy.linalg.tensorinv(self.b1, ind=1) tf_b = function([B], [Bi]) tf_b1 = function([B], [Bi1]) t_binv = tf_b(self.b) t_binv1 = tf_b1(self.b1) assert _allclose(t_binv, n_binv) assert _allclose(t_binv1, n_binv1)
def test_perform(self): if not imported_scipy: raise SkipTest('kron tests need the scipy package to be installed') for shp0 in [(2,), (2, 3), (2, 3, 4), (2, 3, 4, 5)]: x = tensor.tensor(dtype='floatX', broadcastable=(False,) * len(shp0)) a = numpy.asarray(self.rng.rand(*shp0)).astype(config.floatX) for shp1 in [(6,), (6, 7), (6, 7, 8), (6, 7, 8, 9)]: if len(shp0) + len(shp1) == 2: continue y = tensor.tensor(dtype='floatX', broadcastable=(False,) * len(shp1)) f = function([x, y], kron(x, y)) b = self.rng.rand(*shp1).astype(config.floatX) out = f(a, b) # Newer versions of scipy want 4 dimensions at least, # so we have to add a dimension to a and flatten the result. if len(shp0) + len(shp1) == 3: scipy_val = scipy.linalg.kron( a[numpy.newaxis, :], b).flatten() else: scipy_val = scipy.linalg.kron(a, b) utt.assert_allclose(out, scipy_val)
def real_eig(M,eps=1e-9): ''' This code expects a real hermetian matrix and should throw a ValueError if not. This is probably redundant to the scipy eigh function. Do not use. ''' if not (type(M)==np.ndarray): raise ValueError("Expected array; type is %s"%type(M)) if np.any(np.abs(np.imag(M))>eps): raise ValueError("Matrix has imaginary values >%0.2e; will not extract real eigenvalues"%eps) M = np.real(M) w,v = np.linalg.eig(M) if np.any(abs(np.imag(w))>eps): raise ValueError('Eigenvalues with imaginary part >%0.2e; matrix has complex eigenvalues'%eps) w = np.real(w) order = np.argsort(w) w = w[order] v = v[:,order] return w,v
def _getAplus(A): ''' Please see the documentation for nearPDHigham ''' eigval, eigvec = np.linalg.eig(A) Q = np.matrix(eigvec) xdiag = np.matrix(np.diag(np.maximum(eigval, 0))) return Q*xdiag*Q.T
def bench_on(runner, sym, Ns, trials, dtype=None): global args, kernel, out, mkl_layer prepare = globals().get("prepare_"+sym, prepare_default) kernel = globals().get("kernel_"+sym, None) if not kernel: kernel = getattr(np.linalg, sym) out_lvl = runner.__doc__.split('.')[0].strip() func_s = kernel.__doc__.split('.')[0].strip() log.debug('Preparing input data for %s (%s).. ' % (sym, func_s)) args = [prepare(int(i)) for i in Ns] it = range(len(Ns)) # pprint(Ns) out = np.empty(shape=(len(Ns), trials)) b = body(trials) tic, toc = (0, 0) log.debug('Warming up %s (%s).. ' % (sym, func_s)) runner(range(1000), empty_work) kernel(*args[0]) runner(range(1000), empty_work) log.debug('Benchmarking %s on %s: ' % (func_s, out_lvl)) gc_old = gc.isenabled() # gc.disable() tic = time.time() runner(it, b) toc = time.time() - tic if gc_old: gc.enable() if 'reused_pool' in globals(): del globals()['reused_pool'] #calculate average time and min time and also keep track of outliers (max time in the loop) min_time = np.amin(out) max_time = np.amax(out) mean_time = np.mean(out) stdev_time = np.std(out) #print("Min = %.5f, Max = %.5f, Mean = %.5f, stdev = %.5f " % (min_time, max_time, mean_time, stdev_time)) #final_times = [min_time, max_time, mean_time, stdev_time] print('## %s: Outter:%s, Inner:%s, Wall seconds:%f\n' % (sym, out_lvl, mkl_layer, float(toc))) return out
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 get_precision(self): """Compute data precision matrix with the generative model. Equals the inverse of the covariance but computed with the matrix inversion lemma for efficiency. Returns ------- precision : array, shape=(n_features, n_features) Estimated precision of data. """ n_features = self.components_.shape[1] # handle corner cases first if self.n_components_ == 0: return np.eye(n_features) / self.noise_variance_ if self.n_components_ == n_features: return linalg.inv(self.get_covariance()) # Get precision using matrix inversion lemma components_ = self.components_ exp_var = self.explained_variance_ exp_var_diff = np.maximum(exp_var - self.noise_variance_, 0.) precision = np.dot(components_, components_.T) / self.noise_variance_ precision.flat[::len(precision) + 1] += 1. / exp_var_diff precision = np.dot(components_.T, np.dot(linalg.inv(precision), components_)) precision /= -(self.noise_variance_ ** 2) precision.flat[::len(precision) + 1] += 1. / self.noise_variance_ return precision
def __init__(self, dimension, lazy_update_gap=0, constant_trace='', randn=np.random.randn, eigenmethod=np.linalg.eigh): try: self.dimension = len(dimension) standard_deviations = np.asarray(dimension) except TypeError: self.dimension = dimension standard_deviations = np.ones(dimension) assert len(standard_deviations) == self.dimension # prevent equal eigenvals, a hack for np.linalg: self.C = np.diag(standard_deviations**2 * np.exp((1e-4 / self.dimension) * np.arange(self.dimension))) "covariance matrix" self.lazy_update_gap = lazy_update_gap self.constant_trace = constant_trace self.randn = randn self.eigenmethod = eigenmethod self.B = np.eye(self.dimension) "columns, B.T[i] == B[:, i], are eigenvectors of C" self.D = np.diag(self.C)**0.5 # we assume that C is yet diagonal idx = self.D.argsort() self.D = self.D[idx] self.B = self.B[:, idx] "axis lengths, roots of eigenvalues, sorted" self._inverse_root_C = None # see transform_inv... self.last_update = 0 self.count_tell = 0 self.count_eigen = 0
def rotation_from_matrix(matrix): """Return rotation angle and axis from rotation matrix. >>> angle = (random.random() - 0.5) * (2*math.pi) >>> direc = numpy.random.random(3) - 0.5 >>> point = numpy.random.random(3) - 0.5 >>> R0 = rotation_matrix(angle, direc, point) >>> angle, direc, point = rotation_from_matrix(R0) >>> R1 = rotation_matrix(angle, direc, point) >>> is_same_transform(R0, R1) True """ R = numpy.array(matrix, dtype=numpy.float64, copy=False) R33 = R[:3, :3] # direction: unit eigenvector of R33 corresponding to eigenvalue of 1 w, W = numpy.linalg.eig(R33.T) i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0] if not len(i): raise ValueError("no unit eigenvector corresponding to eigenvalue 1") direction = numpy.real(W[:, i[-1]]).squeeze() # point: unit eigenvector of R33 corresponding to eigenvalue of 1 w, Q = numpy.linalg.eig(R) i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0] if not len(i): raise ValueError("no unit eigenvector corresponding to eigenvalue 1") point = numpy.real(Q[:, i[-1]]).squeeze() point /= point[3] # rotation angle depending on direction cosa = (numpy.trace(R33) - 1.0) / 2.0 if abs(direction[2]) > 1e-8: sina = (R[1, 0] + (cosa-1.0)*direction[0]*direction[1]) / direction[2] elif abs(direction[1]) > 1e-8: sina = (R[0, 2] + (cosa-1.0)*direction[0]*direction[2]) / direction[1] else: sina = (R[2, 1] + (cosa-1.0)*direction[1]*direction[2]) / direction[0] angle = math.atan2(sina, cosa) return angle, direction, point # Function to translate handshape coding to degrees of rotation, adduction, flexion
def vector_norm(data, axis=None, out=None): """Return length, i.e. Euclidean norm, of ndarray along axis. >>> v = numpy.random.random(3) >>> n = vector_norm(v) >>> numpy.allclose(n, numpy.linalg.norm(v)) True >>> v = numpy.random.rand(6, 5, 3) >>> n = vector_norm(v, axis=-1) >>> numpy.allclose(n, numpy.sqrt(numpy.sum(v*v, axis=2))) True >>> n = vector_norm(v, axis=1) >>> numpy.allclose(n, numpy.sqrt(numpy.sum(v*v, axis=1))) True >>> v = numpy.random.rand(5, 4, 3) >>> n = numpy.empty((5, 3)) >>> vector_norm(v, axis=1, out=n) >>> numpy.allclose(n, numpy.sqrt(numpy.sum(v*v, axis=1))) True >>> vector_norm([]) 0.0 >>> vector_norm([1]) 1.0 """ data = numpy.array(data, dtype=numpy.float64, copy=True) if out is None: if data.ndim == 1: return math.sqrt(numpy.dot(data, data)) data *= data out = numpy.atleast_1d(numpy.sum(data, axis=axis)) numpy.sqrt(out, out) return out else: data *= data numpy.sum(data, axis=axis, out=out) numpy.sqrt(out, out)
def expms(A, eig=np.linalg.eigh): """matrix exponential for a symmetric matrix""" # TODO: check that this works reliably for low rank matrices # first: symmetrize A D, B = eig(A) return np.dot(B, (np.exp(D) * B).T)
def likelihood(x, m=None, Cinv=None, sigma=1, detC=None): """return likelihood of x for the normal density N(m, sigma**2 * Cinv**-1)""" # testing: MC integrate must be one: mean(p(x_i)) * volume(where x_i are uniformely sampled) # for i in xrange(3): print mean([cma.likelihood(20*r-10, dim * [0], None, 3) for r in rand(10000,dim)]) * 20**dim if m is None: dx = x else: dx = x - m # array(x) - array(m) n = len(x) s2pi = (2 * np.pi)**(n / 2.) if Cinv is None: return exp(-sum(dx**2) / sigma**2 / 2) / s2pi / sigma**n if detC is None: detC = 1. / np.linalg.linalg.det(Cinv) return exp(-np.dot(dx, np.dot(Cinv, dx)) / sigma**2 / 2) / s2pi / abs(detC)**0.5 / sigma**n
def rmsd(X, Y): """ Calculate the root mean squared deviation (RMSD) using Kabsch' formula. @param X: (n, d) input vector @type X: numpy array @param Y: (n, d) input vector @type Y: numpy array @return: rmsd value between the input vectors @rtype: float """ from numpy import sum, dot, sqrt, clip, average from numpy.linalg import svd, det X = X - X.mean(0) Y = Y - Y.mean(0) R_x = sum(X ** 2) R_y = sum(Y ** 2) V, L, U = svd(dot(Y.T, X)) if det(dot(V, U)) < 0.: L[-1] *= -1 return sqrt(clip(R_x + R_y - 2 * sum(L), 0., 1e300) / len(X))
def wrmsd(X, Y, w): """ Calculate the weighted root mean squared deviation (wRMSD) using Kabsch' formula. @param X: (n, d) input vector @type X: numpy array @param Y: (n, d) input vector @type Y: numpy array @param w: input weights @type w: numpy array @return: rmsd value between the input vectors @rtype: float """ from numpy import sum, dot, sqrt, clip, average from numpy.linalg import svd ## normalize weights w = w / w.sum() X = X - dot(w, X) Y = Y - dot(w, Y) R_x = sum(X.T ** 2 * w) R_y = sum(Y.T ** 2 * w) L = svd(dot(Y.T * w, X))[1] return sqrt(clip(R_x + R_y - 2 * sum(L), 0., 1e300))
def is_mirror_image(X, Y): """ Check if two configurations X and Y are mirror images (i.e. their optimal superposition involves a reflection). @param X: n x 3 input vector @type X: numpy array @param Y: n x 3 input vector @type Y: numpy array @rtype: bool """ from numpy.linalg import det, svd ## center configurations X = X - numpy.mean(X, 0) Y = Y - numpy.mean(Y, 0) ## SVD of correlation matrix V, L, U = svd(numpy.dot(numpy.transpose(X), Y)) #@UnusedVariable R = numpy.dot(V, U) return det(R) < 0
def eig(a): u,v = np.linalg.eig(a) return u.T
def sparse_eigenspectrum(sparse_operator): """Perform a dense diagonalization. Returns: eigenspectrum: The lowest eigenvalues in a numpy array. """ dense_operator = sparse_operator.todense() if is_hermitian(sparse_operator): eigenspectrum = numpy.linalg.eigvalsh(dense_operator) else: eigenspectrum = numpy.linalg.eigvals(dense_operator) return numpy.sort(eigenspectrum)
def get_gap(sparse_operator): """Compute gap between lowest eigenvalue and first excited state. Returns: A real float giving eigenvalue gap. """ if not is_hermitian(sparse_operator): raise ValueError('sparse_operator must be Hermitian.') values, _ = scipy.sparse.linalg.eigsh( sparse_operator, 2, which='SA', maxiter=1e7) gap = abs(values[1] - values[0]) return gap
def perpedndicular1(v): """calculate the perpendicular unit vector""" return numpy.array((-v[1], v[0])) / numpy.linalg.norm((-v[1], v[0]))
def normalize(v): return v / numpy.linalg.norm(v)
def test_qr_modes(): rng = numpy.random.RandomState(utt.fetch_seed()) A = tensor.matrix("A", dtype=theano.config.floatX) a = rng.rand(4, 4).astype(theano.config.floatX) f = function([A], qr(A)) t_qr = f(a) n_qr = numpy.linalg.qr(a) assert _allclose(n_qr, t_qr) for mode in ["reduced", "r", "raw"]: f = function([A], qr(A, mode)) t_qr = f(a) n_qr = numpy.linalg.qr(a, mode) if isinstance(n_qr, (list, tuple)): assert _allclose(n_qr[0], t_qr[0]) assert _allclose(n_qr[1], t_qr[1]) else: assert _allclose(n_qr, t_qr) try: n_qr = numpy.linalg.qr(a, "complete") f = function([A], qr(A, "complete")) t_qr = f(a) assert _allclose(n_qr, t_qr) except TypeError as e: assert "name 'complete' is not defined" in str(e)
def test_svd(): rng = numpy.random.RandomState(utt.fetch_seed()) A = tensor.matrix("A", dtype=theano.config.floatX) U, V, T = svd(A) fn = function([A], [U, V, T]) a = rng.rand(4, 4).astype(theano.config.floatX) n_u, n_v, n_t = numpy.linalg.svd(a) t_u, t_v, t_t = fn(a) assert _allclose(n_u, t_u) assert _allclose(n_v, t_v) assert _allclose(n_t, t_t)
def test_inverse_singular(): singular = numpy.array([[1, 0, 0]] + [[0, 1, 0]] * 2, dtype=theano.config.floatX) a = tensor.matrix() f = function([a], matrix_inverse(a)) try: f(singular) except numpy.linalg.LinAlgError: return assert False
def test_wrong_coefficient_matrix(self): x = tensor.vector() y = tensor.vector() z = tensor.scalar() b = theano.tensor.nlinalg.lstsq()(x, y, z) f = function([x, y, z], b) self.assertRaises(numpy.linalg.linalg.LinAlgError, f, [2, 1], [2, 1], 1)
def test_wrong_rcond_dimension(self): x = tensor.vector() y = tensor.vector() z = tensor.vector() b = theano.tensor.nlinalg.lstsq()(x, y, z) f = function([x, y, z], b) self.assertRaises(numpy.linalg.LinAlgError, f, [2, 1], [2, 1], [2, 1])
def test_numpy_compare(self): rng = numpy.random.RandomState(utt.fetch_seed()) A = tensor.matrix("A", dtype=theano.config.floatX) Q = matrix_power(A, 3) fn = function([A], [Q]) a = rng.rand(4, 4).astype(theano.config.floatX) n_p = numpy.linalg.matrix_power(a, 3) t_p = fn(a) assert numpy.allclose(n_p, t_p)
def test_eigvalsh_grad(): if not imported_scipy: raise SkipTest("Scipy needed for the geigvalsh op.") import scipy.linalg rng = numpy.random.RandomState(utt.fetch_seed()) a = rng.randn(5, 5) a = a + a.T b = 10 * numpy.eye(5, 5) + rng.randn(5, 5) tensor.verify_grad(lambda a, b: eigvalsh(a, b).dot([1, 2, 3, 4, 5]), [a, b], rng=numpy.random)
def test_solve_correctness(self): if not imported_scipy: raise SkipTest("Scipy needed for the Cholesky and Solve ops.") rng = numpy.random.RandomState(utt.fetch_seed()) A = theano.tensor.matrix() b = theano.tensor.matrix() y = self.op(A, b) gen_solve_func = theano.function([A, b], y) cholesky_lower = Cholesky(lower=True) L = cholesky_lower(A) y_lower = self.op(L, b) lower_solve_func = theano.function([L, b], y_lower) cholesky_upper = Cholesky(lower=False) U = cholesky_upper(A) y_upper = self.op(U, b) upper_solve_func = theano.function([U, b], y_upper) b_val = numpy.asarray(rng.rand(5, 1), dtype=config.floatX) # 1-test general case A_val = numpy.asarray(rng.rand(5, 5), dtype=config.floatX) # positive definite matrix: A_val = numpy.dot(A_val.transpose(), A_val) assert numpy.allclose(scipy.linalg.solve(A_val, b_val), gen_solve_func(A_val, b_val)) # 2-test lower traingular case L_val = scipy.linalg.cholesky(A_val, lower=True) assert numpy.allclose(scipy.linalg.solve_triangular(L_val, b_val, lower=True), lower_solve_func(L_val, b_val)) # 3-test upper traingular case U_val = scipy.linalg.cholesky(A_val, lower=False) assert numpy.allclose(scipy.linalg.solve_triangular(U_val, b_val, lower=False), upper_solve_func(U_val, b_val))
def test_expm(): if not imported_scipy: raise SkipTest("Scipy needed for the expm op.") rng = numpy.random.RandomState(utt.fetch_seed()) A = rng.randn(5, 5).astype(config.floatX) ref = scipy.linalg.expm(A) x = tensor.matrix() m = expm(x) expm_f = function([x], m) val = expm_f(A) numpy.testing.assert_array_almost_equal(val, ref)
def rcond(x): ''' Reciprocal condition number ''' return 1./np.linalg.cond(x)
def check_covmat_fast(C,N=None,eps=1e-6): ''' Verify that matrix M is a size NxN precision or covariance matirx ''' if not type(C)==np.ndarray: raise ValueError("Covariance matrix should be a 2D numpy array") if not len(C.shape)==2: raise ValueError("Covariance matrix should be a 2D numpy array") if N is None: N = C.shape[0] if not C.shape==(N,N): raise ValueError("Expected size %d x %d matrix"%(N,N)) if np.any(~np.isreal(C)): raise ValueError("Covariance matrices should not contain complex numbers") C = np.real(C) if np.any(~np.isfinite(C)): raise ValueError("Covariance matrix contains NaN or ±inf!") if not np.all(np.abs(C-C.T)<eps): raise ValueError("Covariance matrix is not symmetric up to precision %0.1e"%eps) try: ch = chol(C) except numpy.linalg.linalg.LinAlgError: # Check smallest eigenvalue if cholesky fails mine = np.real(scipy.linalg.decomp.eigh(C,eigvals=(0,0))[0][0]) if np.any(mine<-eps): raise ValueError('Covariance matrix contains eigenvalue %0.3e<%0.3e'%(mine,-eps)) if (mine<eps): C = C + np.eye(N)*(eps-mine) C = 0.5*(C+C.T) return C
def rsolve(a,b): ''' wraps solve, applies to right hand solution solves system x A = B ''' return scipy.linalg.solve(b.T,a.T).T
def lwlr(testPoint, xMat, yMat, k=1.0): m = np.shape(xMat)[0] weights = np.matrix(np.eye(m)) # ?????? for j in range(m): diffMat = testPoint-xMat[j, :] weights[j, j] = np.exp(diffMat*diffMat.T/(-2.0*k**2)) # ??? print weights xTx = xMat.T*(weights*xMat) if np.linalg.det(xTx) == 0.0: print 'This matrix is singular, cannot do inverse' return ws = xTx.I*(xMat.T*(weights*yMat)) return testPoint*ws
def ridgeRegres(xMat, yMat, lam=0.2): xTx = xMat.T*xMat denom = xTx+np.eye(np.shape(xMat)[1])*lam if np.linalg.det(denom) == 0.0: print 'This matrix is singular, cannot do inverse' return ws = denom.I*(xMat.T*yMat) return ws
