Python numpy.linalg() 使用实例

The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.

Example 1

```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 ```

Example 2

```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 ```

Example 3

```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() ```

Example 4

```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,
) ```

Example 5

```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 ```

Example 6

```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)) ```

Example 7

```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) ```

Example 8

```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) ```

Example 9

```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',
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',
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) ```

Example 10

```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 ```

Example 11

```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 ```

Example 12

```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 ```

Example 13

```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 ```

Example 14

```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 ```

Example 15

```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 ```

Example 16

```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 ```

Example 17

```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) ```

Example 18

```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 ```

Example 19

```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) ```

Example 20

```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 ```

Example 21

```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)) ```

Example 22

```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)) ```

Example 23

```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 ```

Example 24

```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) ```

Example 25

```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 ```

Example 26

```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) ```

Example 27

```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 ```

Example 28

```def eig(a):
u,v = np.linalg.eig(a)
return u.T ```

Example 29

```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) ```

Example 30

```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 ```

Example 31

```def perpedndicular1(v):
"""calculate the perpendicular unit vector"""
return numpy.array((-v[1], v[0])) / numpy.linalg.norm((-v[1], v[0])) ```

Example 32

```def normalize(v):
return v / numpy.linalg.norm(v) ```

Example 33

```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) ```

Example 34

```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) ```

Example 35

```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) ```

Example 36

```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 ```

Example 37

```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) ```

Example 38

```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]) ```

Example 39

```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) ```

Example 40

```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) ```

Example 41

```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)) ```

Example 42

```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) ```

Example 43

```def rcond(x):
'''
Reciprocal condition number
'''
return 1./np.linalg.cond(x) ```

Example 44

```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 ```

Example 45

```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 ```

Example 46

```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) ```

Example 47

```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 ```

Example 48

```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 ```

Example 49

```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 ```

Example 50

```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) ```