Python numpy.linalg() Examples

The following are 30 code examples for showing how to use numpy.linalg(). These examples are extracted from open source projects. You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example.

You may check out the related API usage on the sidebar.

You may also want to check out all available functions/classes of the module numpy , or try the search function .

Example 1
Project: pyscf   Author: pyscf   File: 001-remove_linear_dep.py    License: Apache License 2.0 6 votes vote down vote up
def eig(h, s):
    from pyscf import symm
    nirrep = len(mol.symm_orb)
    h = symm.symmetrize_matrix(h, mol.symm_orb)
    s = symm.symmetrize_matrix(s, mol.symm_orb)
    cs = []
    es = []
#
# Linear dependency are removed by looping over different symmetry irreps.
#
    for ir in range(nirrep):
        d, t = numpy.linalg.eigh(s[ir])
        x = t[:,d>1e-8] / numpy.sqrt(d[d>1e-8])
        xhx = reduce(numpy.dot, (x.T, h[ir], x))
        e, c = numpy.linalg.eigh(xhx)
        cs.append(reduce(numpy.dot, (mol.symm_orb[ir], x, c)))
        es.append(e)
    e = numpy.hstack(es)
    c = numpy.hstack(cs)
    return e, c 
Example 2
Project: pyscf   Author: pyscf   File: 42-remove_linear_dep.py    License: Apache License 2.0 6 votes vote down vote up
def eig(h, s):
    from pyscf import symm
    nirrep = len(mol.symm_orb)
    h = symm.symmetrize_matrix(h, mol.symm_orb)
    s = symm.symmetrize_matrix(s, mol.symm_orb)
    cs = []
    es = []
#
# Linear dependency are removed by looping over different symmetry irreps.
#
    for ir in range(nirrep):
        d, t = numpy.linalg.eigh(s[ir])
        x = t[:,d>1e-8] / numpy.sqrt(d[d>1e-8])
        xhx = reduce(numpy.dot, (x.T, h[ir], x))
        e, c = numpy.linalg.eigh(xhx)
        cs.append(reduce(numpy.dot, (mol.symm_orb[ir], x, c)))
        es.append(e)
    e = numpy.hstack(es)
    c = numpy.hstack(cs)
    return e, c 
Example 3
Project: SparseSC   Author: microsoft   File: fit_fast.py    License: MIT License 6 votes vote down vote up
def _weights(V , X_treated, X_control, w_pen):
    V = np.diag(V) #make square
    #weights = np.zeros((X_control.shape[0], X_treated.shape[0]))
    w_pen_mat = 2 * w_pen * np.diag(np.ones(X_control.shape[0]))
    A = X_control.dot(2 * V).dot(X_control.T) + w_pen_mat  # 5
    B = (
        X_treated.dot(2 * V).dot(X_control.T).T + 2 * w_pen / X_control.shape[0]
    )  # 6
    try:
        b = scipy.linalg.solve(A, B)
    except scipy.linalg.LinAlgError as exc:
        print("Unique weights not possible.")
        if w_pen == 0:
            print("Try specifying a very small w_pen rather than 0.")
        raise exc
    return b 
Example 4
Project: D-VAE   Author: muhanzhang   File: test_slinalg.py    License: MIT License 6 votes vote down vote up
def test_eigvalsh():
    if not imported_scipy:
        raise SkipTest("Scipy needed for the geigvalsh op.")
    import scipy.linalg

    A = theano.tensor.dmatrix('a')
    B = theano.tensor.dmatrix('b')
    f = function([A, B], eigvalsh(A, B))

    rng = numpy.random.RandomState(utt.fetch_seed())
    a = rng.randn(5, 5)
    a = a + a.T
    for b in [10 * numpy.eye(5, 5) + rng.randn(5, 5)]:
        w = f(a, b)
        refw = scipy.linalg.eigvalsh(a, b)
        numpy.testing.assert_array_almost_equal(w, refw)

    # We need to test None separatly, as otherwise DebugMode will
    # complain, as this isn't a valid ndarray.
    b = None
    B = theano.tensor.NoneConst
    f = function([A], eigvalsh(A, B))
    w = f(a)
    refw = scipy.linalg.eigvalsh(a, b)
    numpy.testing.assert_array_almost_equal(w, refw) 
Example 5
Project: D-VAE   Author: muhanzhang   File: test_slinalg.py    License: MIT License 6 votes vote down vote up
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) 
Example 6
Project: D-VAE   Author: muhanzhang   File: test_nlinalg.py    License: MIT License 6 votes vote down vote up
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 7
Project: vnpy_crypto   Author: birforce   File: linalg.py    License: MIT License 6 votes vote down vote up
def logdet_symm(m, check_symm=False):
    """
    Return log(det(m)) asserting positive definiteness of m.

    Parameters
    ----------
    m : array-like
        2d array that is positive-definite (and symmetric)

    Returns
    -------
    logdet : float
        The log-determinant of m.
    """
    from scipy import linalg
    if check_symm:
        if not np.all(m == m.T):  # would be nice to short-circuit check
            raise ValueError("m is not symmetric.")
    c, _ = linalg.cho_factor(m, lower=True)
    return 2*np.sum(np.log(c.diagonal())) 
Example 8
Project: OpenFermion   Author: quantumlib   File: _davidson_test.py    License: Apache License 2.0 6 votes vote down vote up
def setUp(self):
        """Sets up all variables needed for Davidson class."""
        dimension = 10
        matrix = generate_matrix(dimension)

        def mat_vec(vec):
            """Trivial matvec with a numpy matrix."""
            return numpy.dot(matrix, vec)

        self.linear_operator = scipy.sparse.linalg.LinearOperator(
            (dimension, dimension), matvec=mat_vec)
        self.diagonal = numpy.diag(matrix)

        self.davidson = Davidson(linear_operator=self.linear_operator,
                                 linear_operator_diagonal=self.diagonal)

        self.matrix = matrix
        self.initial_guess = numpy.eye(self.matrix.shape[0], 10)

        self.eigen_values = numpy.array([
            1.15675714, 1.59132505, 2.62268014, 4.44533793, 5.3722743,
            5.54393114, 7.73652405, 8.50089897, 9.4229309, 15.54405993,
        ]) 
Example 9
Project: OpenFermion   Author: quantumlib   File: _davidson_test.py    License: Apache License 2.0 6 votes vote down vote up
def setUp(self):
        """Sets up all variables needed for SparseDavidson class."""
        logging.basicConfig(level=logging.INFO)
        self.dimension = 1000
        self.sparse_matrix = generate_sparse_matrix(self.dimension)
        self.davidson_options = DavidsonOptions(max_subspace=100,
                                                max_iterations=50,
                                                real_only=True)

        # Checks for built-in eigh() function.
        self.eigen_values, self.eigen_vectors = numpy.linalg.eigh(
            self.sparse_matrix)
        self.assertAlmostEqual(get_difference(
            self.sparse_matrix, self.eigen_values, self.eigen_vectors), 0)

        # Makes sure eigenvalues are sorted.
        self.eigen_values = sorted(self.eigen_values) 
Example 10
Project: OpenFermion   Author: quantumlib   File: _sparse_tools.py    License: Apache License 2.0 6 votes vote down vote up
def variance(operator, state):
    """Compute variance of operator with a state.

    Args:
        operator(scipy.sparse.spmatrix or scipy.sparse.linalg.LinearOperator):
            The operator whose expectation value is desired.
        state(numpy.ndarray or scipy.sparse.spmatrix): A numpy array
            representing a pure state or a sparse matrix representing a density
            matrix.

    Returns:
        A complex number giving the variance.

    Raises:
        ValueError: Input state has invalid format.
    """
    return (expectation(operator ** 2, state) -
            expectation(operator, state) ** 2) 
Example 11
Project: OpenFermion   Author: quantumlib   File: _sparse_tools.py    License: Apache License 2.0 6 votes vote down vote up
def get_gap(sparse_operator, initial_guess=None):
    """Compute gap between lowest eigenvalue and first excited state.

    Args:
        sparse_operator (LinearOperator): Operator to find the ground state of.
        initial_guess (ndarray): Initial guess for eigenspace.  A good
            guess dramatically reduces the cost required to converge.
    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, k=2, v0=initial_guess, which='SA', maxiter=1e7)

    gap = abs(values[1] - values[0])
    return gap 
Example 12
Project: OpenFermion   Author: quantumlib   File: _linear_qubit_operator.py    License: Apache License 2.0 6 votes vote down vote up
def generate_linear_qubit_operator(qubit_operator, n_qubits=None, options=None):
    """ Generates a LinearOperator from a QubitOperator.

    Args:
        qubit_operator(QubitOperator): A qubit operator to be applied on
            vectors.
        n_qubits(int): The total number of qubits
        options(LinearQubitOperatorOptions): Options for the
            ParallelLinearQubitOperator.
    Returns:
        linear_operator(scipy.sparse.linalg.LinearOperator): A linear operator.
    """
    if options is None:
        linear_operator = LinearQubitOperator(qubit_operator, n_qubits)
    else:
        linear_operator = ParallelLinearQubitOperator(
            qubit_operator, n_qubits, options)
    return linear_operator 
Example 13
Project: OpenFermion   Author: quantumlib   File: _davidson.py    License: Apache License 2.0 6 votes vote down vote up
def generate_random_vectors(row, col, real_only=False):
    """Generates orthonormal random vectors with col columns.

    Args:
        row(int): Number of rows for the vectors.
        col(int): Number of columns for the vectors.
        real_only(bool): Real vectors or complex ones.

    Returns:
        random_vectors(numpy.ndarray(complex)): Orthonormal random vectors.
    """
    random_vectors = numpy.random.rand(row, col)
    if not real_only:
        random_vectors = random_vectors + numpy.random.rand(row, col) * 1.0j
    random_vectors = scipy.linalg.orth(random_vectors)
    return random_vectors 
Example 14
Project: attention-lvcsr   Author: rizar   File: test_slinalg.py    License: MIT License 6 votes vote down vote up
def test_eigvalsh():
    if not imported_scipy:
        raise SkipTest("Scipy needed for the geigvalsh op.")
    import scipy.linalg

    A = theano.tensor.dmatrix('a')
    B = theano.tensor.dmatrix('b')
    f = function([A, B], eigvalsh(A, B))

    rng = numpy.random.RandomState(utt.fetch_seed())
    a = rng.randn(5, 5)
    a = a + a.T
    for b in [10 * numpy.eye(5, 5) + rng.randn(5, 5)]:
        w = f(a, b)
        refw = scipy.linalg.eigvalsh(a, b)
        numpy.testing.assert_array_almost_equal(w, refw)

    # We need to test None separatly, as otherwise DebugMode will
    # complain, as this isn't a valid ndarray.
    b = None
    B = theano.tensor.NoneConst
    f = function([A], eigvalsh(A, B))
    w = f(a)
    refw = scipy.linalg.eigvalsh(a, b)
    numpy.testing.assert_array_almost_equal(w, refw) 
Example 15
Project: attention-lvcsr   Author: rizar   File: test_slinalg.py    License: MIT License 6 votes vote down vote up
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) 
Example 16
Project: attention-lvcsr   Author: rizar   File: test_nlinalg.py    License: MIT License 6 votes vote down vote up
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 17
Project: neural-belief-tracker   Author: nmrksic   File: nbt.py    License: Apache License 2.0 6 votes vote down vote up
def xavier_vector(word, D=300):
    """
    Returns a D-dimensional vector for the word. 

    We hash the word to always get the same vector for the given word. 
    """

    seed_value = hash_string(word)
    numpy.random.seed(seed_value)

    neg_value = - math.sqrt(6)/math.sqrt(D)
    pos_value = math.sqrt(6)/math.sqrt(D)

    rsample = numpy.random.uniform(low=neg_value, high=pos_value, size=(D,))
    norm = numpy.linalg.norm(rsample)
    rsample_normed = rsample/norm

    return rsample_normed 
Example 18
Project: Uranium   Author: Ultimaker   File: MeshData.py    License: GNU Lesser General Public License v3.0 6 votes vote down vote up
def calculateNormalsFromIndexedVertices(vertices: numpy.ndarray, indices: numpy.ndarray, face_count: int) -> numpy.ndarray:
    """Calculate the normals of this mesh of triagles using indexes.

    :param vertices: :type{narray} list of vertices as a 1D list of float triples
    :param indices: :type{narray} list of indices as a 1D list of integers
    :param face_count: :type{integer} the number of triangles defined by the indices array
    :return: :type{narray} list normals as a 1D array of floats, each group of 3 floats is a vector
    """

    start_time = time()
    # Numpy magic!
    # First, reset the normals
    normals = numpy.zeros((face_count*3, 3), dtype=numpy.float32)

    for face in indices[0:face_count]:
        normals[face[0]] = numpy.cross(vertices[face[0]] - vertices[face[1]], vertices[face[0]] - vertices[face[2]])
        length = numpy.linalg.norm(normals[face[0]])
        normals[face[0]] /= length
        normals[face[1]] = normals[face[0]]
        normals[face[2]] = normals[face[0]]
    end_time = time()
    Logger.log("d", "Calculating normals took %s seconds", end_time - start_time)
    return normals 
Example 19
Project: FreeCAD_assembly2   Author: hamish2014   File: __init__.py    License: GNU Lesser General Public License v2.1 6 votes vote down vote up
def rotation_matrix_axis_and_angle_2(R, debug=False, errorThreshold=10**-7, msg=None):
    w, v = numpy.linalg.eig(R) #this method is not used at the primary method as numpy.linalg.eig does not return answers in high enough precision
    angle, axis = None, None
    eigErrors = abs(w -1) #errors from 1+0j
    i = (eigErrors == min(eigErrors)).tolist().index(True)
    axis = numpy.real(v[:,i])
    if i != 1:
        angle = arccos2(  numpy.real( w[1] ) )
    else:
        angle = arccos2(  numpy.real( w[0] ) )
    error  = norm(axis_rotation_matrix(angle, *axis) - R)
    if debug: print('rotation_matrix_axis_and_angle error %1.1e' % error)
    if error > errorThreshold:
        angle = -angle
        error = norm(axis_rotation_matrix(angle, *axis) - R)
        if error > errorThreshold:
            R_pickle_str = pickle.dumps(R)
            #R_abs_minus_identity = abs(R) - numpy.eye(3)
            print(R*R.transpose())
            raise ValueError( 'rotation_matrix_axis_and_angle_2: no solution found! locals %s' % str(locals()))
    return axis, angle 
Example 20
Project: pyscf   Author: pyscf   File: 001-remove_linear_dep.py    License: Apache License 2.0 5 votes vote down vote up
def eig(h, s):
    d, t = numpy.linalg.eigh(s)
# Removing the eigenvectors assoicated to the smallest eigenvalue, the new
# basis defined by x matrix has 139 vectors.
    x = t[:,d>1e-8] / numpy.sqrt(d[d>1e-8])
    xhx = reduce(numpy.dot, (x.T, h, x))
    e, c = numpy.linalg.eigh(xhx)
    c = numpy.dot(x, c)
# Return 139 eigenvalues and 139 eigenvectors.
    return e, c
#
# Replacing the default eig function with the above one,  the HF solver
# generate only 139 canonical orbitals
# 
Example 21
Project: pyscf   Author: pyscf   File: 42-remove_linear_dep.py    License: Apache License 2.0 5 votes vote down vote up
def eig(h, s):
    d, t = numpy.linalg.eigh(s)
# Removing the eigenvectors assoicated to the smallest eigenvalue, the new
# basis defined by x matrix has 139 vectors.
    x = t[:,d>1e-8] / numpy.sqrt(d[d>1e-8])
    xhx = reduce(numpy.dot, (x.T, h, x))
    e, c = numpy.linalg.eigh(xhx)
    c = numpy.dot(x, c)
# Return 139 eigenvalues and 139 eigenvectors.
    return e, c
#
# Replacing the default eig function with the above one,  the HF solver
# generate only 139 canonical orbitals
# 
Example 22
Project: tangent   Author: google   File: grads.py    License: Apache License 2.0 5 votes vote down vote up
def adet(z, x):
  """d|A|/dA = adj(A).T

  See  Jacobi's formula: https://en.wikipedia.org/wiki/Jacobi%27s_formula
  """
  adjugate = numpy.linalg.det(x) * numpy.linalg.pinv(x)
  d[x] = d[z] * numpy.transpose(adjugate)


#
# Tangent adjoints
# 
Example 23
Project: recruit   Author: Frank-qlu   File: test_public_api.py    License: Apache License 2.0 5 votes vote down vote up
def test_numpy_linalg():
    bad_results = check_dir(np.linalg)
    assert bad_results == {} 
Example 24
Project: hoggorm   Author: olivertomic   File: statTools.py    License: BSD 2-Clause "Simplified" License 5 votes vote down vote up
def ortho(arr1, arr2):
    """
    This function orthogonalises arr1 with respect to arr2. The function then
    returns orthogonalised array arr1_orth.

    PARAMETERS
    ----------
    arr1 : numpy array
        A numpy array containing some data

    arr2 : numpy array
        A numpy array containing some data

    RETURNS
    -------
    numpy array
        A numpy array holding orthogonalised numpy array ``arr1``.

    Examples
    --------
    some examples

    """

    # Find number of rows, such that identity matrix I can be created
    numberRows = numpy.shape(arr1)[0]
    I = numpy.identity(numberRows, float)

    # Compute transpose of arr1
    arr2_T = numpy.transpose(arr2)

    term1 = numpy.linalg.inv(numpy.dot(arr2_T, arr2))
    term2 = numpy.dot(arr2, term1)
    term3 = numpy.dot(term2, arr2_T)
    arr1_orth = numpy.dot((I - term3), arr1)

    return arr1_orth 
Example 25
Project: hoggorm   Author: olivertomic   File: statTools.py    License: BSD 2-Clause "Simplified" License 5 votes vote down vote up
def matrixRank(arr, tol=1e-8):
    """
    Computes the rank of an array/matrix, i.e. number of linearly independent
    variables. This is not the same as numpy.rank() which only returns the
    number of ways (2-way, 3-way, etc) an array/matrix has.

    PARAMETERS
    ----------
    arrX : numpy array
        A numpy array containing the data

    RETURNS
    -------
    scalar
        Rank of matrix.

    Examples
    --------
    >>> import hoggorm as ho
    >>>
    >>> # Get the rank of the data
    >>> ho.matrixRank(myData)
    >>> 8

    """
    if len(arr.shape) != 2:
        raise ValueError('Input must be a 2-d array or Matrix object')

    s = numpy.linalg.svd(arr, compute_uv=0)
    return numpy.sum(numpy.where(s > tol, 1, 0)) 
Example 26
Project: D-VAE   Author: muhanzhang   File: test_slinalg.py    License: MIT License 5 votes vote down vote up
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 27
Project: D-VAE   Author: muhanzhang   File: test_slinalg.py    License: MIT License 5 votes vote down vote up
def test_solve_correctness(self):
        if not imported_scipy:
            raise SkipTest("Scipy needed for the Cholesky op.")
        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 28
Project: D-VAE   Author: muhanzhang   File: test_slinalg.py    License: MIT License 5 votes vote down vote up
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 29
Project: D-VAE   Author: muhanzhang   File: test_nlinalg.py    License: MIT License 5 votes vote down vote up
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", "full", "economic"]:
        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 30
Project: D-VAE   Author: muhanzhang   File: test_nlinalg.py    License: MIT License 5 votes vote down vote up
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)