Python numpy.linalg.multi_dot() Examples

def test_boson_operator_sparse_multi_mode(self):
        op = BosonOperator('0^ 1 1^ 2')
        res = boson_operator_sparse(op, self.d).toarray()

        b0 = boson_ladder_sparse(3, 0, 0, self.d).toarray()
        b1 = boson_ladder_sparse(3, 1, 0, self.d).toarray()
        b2 = boson_ladder_sparse(3, 2, 0, self.d).toarray()

        expected = multi_dot([b0.T, b1, b1.T, b2])
        self.assertTrue(numpy.allclose(res, expected))

        op = QuadOperator('q0 p0 p1')
        res = boson_operator_sparse(op, self.d, self.hbar).toarray()

        expected = numpy.identity(self.d**2)
        for term in op.terms:
            for i, j in term:
                expected = expected.dot(single_quad_op_sparse(
                    2, i, j, self.hbar, self.d).toarray())
        self.assertTrue(numpy.allclose(res, expected)) 

def __init__(self):
        matrix_pairs = {}
        rng = np.random.RandomState(SEED)
        for s in SIZES:
            for n_terms in [2, 4]:
                for explicit in [True, False]:
                    arrays = [
                        rng.standard_normal((s if t % 2 == 0 else 2 * s,
                                             2 * s if t % 2 == 0 else s))
                        for t in range(n_terms)]
                    matrices_ = [
                        matrices.DenseRectangularMatrix(a) for a in arrays]
                    if explicit:
                        matrix = matrices.MatrixProduct(matrices_)
                    else:
                        matrix = reduce(lambda a, b: a @ b, matrices_)
                    matrix_pairs[(s, n_terms, explicit)] = (
                        matrix, nla.multi_dot(arrays))
        super().__init__(matrix_pairs, rng) 

def test_hadamard_decomposition(self, tol):
        """Tests that the decomposition of the Hadamard gate is correct"""
        op = qml.Hadamard(wires=0)
        res = op.decomposition(0)

        assert len(res) == 3

        assert res[0].name == "PhaseShift"
        assert res[0].wires == [0]  #qml.wires.Wires([0])
        assert res[0].params[0] == np.pi / 2

        assert res[1].name == "RX"
        assert res[1].wires == [0]  #qml.wires.Wires([0])
        assert res[0].params[0] == np.pi / 2
        
        assert res[2].name == "PhaseShift"
        assert res[2].wires == [0]  #qml.wires.Wires([0])
        assert res[0].params[0] == np.pi / 2
        
        decomposed_matrix = np.linalg.multi_dot([i.matrix for i in reversed(res)])
        assert np.allclose(decomposed_matrix, op.matrix, atol=tol, rtol=0) 

def test_y_decomposition(self, tol):
        """Tests that the decomposition of the PauliY is correct"""
        op = qml.PauliY(wires=0)
        res = op.decomposition(0)

        assert len(res) == 3

        assert res[0].name == "PhaseShift"
        assert res[0].wires == [0]  #qml.wires.Wires([0])
        assert res[0].params[0] == np.pi / 2
        
        assert res[1].name == "RY"
        assert res[1].wires == [0]  #qml.wires.Wires([0])
        assert res[1].params[0] == np.pi
        
        assert res[2].name == "PhaseShift"
        assert res[2].wires == [0]  #qml.wires.Wires([0])
        assert res[2].params[0] == np.pi / 2
        
        decomposed_matrix = np.linalg.multi_dot([i.matrix for i in reversed(res)])
        assert np.allclose(decomposed_matrix, op.matrix, atol=tol, rtol=0) 

def test_x_decomposition(self, tol):
        """Tests that the decomposition of the PauliX is correct"""
        op = qml.PauliX(wires=0)
        res = op.decomposition(0)

        assert len(res) == 3

        assert res[0].name == "PhaseShift"
        assert res[0].wires == [0]  #qml.wires.Wires([0])
        assert res[0].params[0] == np.pi / 2
        
        assert res[1].name == "RX"
        assert res[1].wires == [0]  #qml.wires.Wires([0])
        assert res[1].params[0] == np.pi
        
        assert res[2].name == "PhaseShift"
        assert res[2].wires == [0]  #qml.wires.Wires([0])
        assert res[2].params[0] == np.pi / 2

        decomposed_matrix = np.linalg.multi_dot([i.matrix for i in reversed(res)])
        assert np.allclose(decomposed_matrix, op.matrix, atol=tol, rtol=0) 

def test_diagonalization(self, obs, mat, eigs, tol):
        """Test the method transforms standard observables into the Z-gate."""
        ob = obs(wires=0)
        A = ob.matrix

        diag_gates = ob.diagonalizing_gates()
        U = np.eye(2)

        if diag_gates:
            mats = [i.matrix for i in diag_gates]
            # Need to revert the order in which the matrices are applied such that they adhere to the order
            # of matrix multiplication
            # E.g. for PauliY: [PauliZ(wires=self.wires), S(wires=self.wires), Hadamard(wires=self.wires)]
            # becomes Hadamard @ S @ PauliZ, where @ stands for matrix multiplication
            mats = mats[::-1]
            U = multi_dot([np.eye(2)] + mats)

        res = U @ A @ U.conj().T
        expected = np.diag(eigs)
        assert np.allclose(res, expected, atol=tol, rtol=0) 

def __init__(self):
        matrix_pairs = {}
        rng = np.random.RandomState(SEED)
        for s in SIZES:
            for n_terms in [2, 5]:
                for explicit in [True, False]:
                    arrays = [
                        rng.standard_normal((s, s)) for _ in range(n_terms)]
                    matrices_ = [
                        matrices.DenseSquareMatrix(a) for a in arrays]
                    if explicit:
                        matrix = matrices.InvertibleMatrixProduct(matrices_)
                    else:
                        matrix = reduce(lambda a, b: a @ b, matrices_)
                    matrix_pairs[(s, n_terms, explicit)] = (
                        matrix, nla.multi_dot(arrays))
        super().__init__(matrix_pairs, rng) 

def assoc_test(weights, gwas, ldmat, heterogeneity=False):
    """
    TWAS association test.

    :param weights: numpy.ndarray of eQTL weights
    :param gwas: pyfocus.GWAS object
    :param ldmat: numpy.ndarray LD matrix
    :param heterogeneity:  bool estimate variance from multiplicative random effect

    :return: tuple (beta, se)
    """

    p = ldmat.shape[0]
    assoc = np.dot(weights, gwas.Z)
    if heterogeneity:
        resid = assoc - gwas.Z
        resid_var = mdot([resid, lin.pinvh(ldmat), resid]) / p
    else:
        resid_var = 1

    se = np.sqrt(resid_var * mdot([weights, ldmat, weights]))

    return assoc, se 

def estimate_cor(wmat, ldmat, intercept=False):
    """
    Estimate the sample correlation structure for predicted expression.

    :param wmat: numpy.ndarray eQTL weight matrix for a risk region
    :param ldmat: numpy.ndarray LD matrix for a risk region
    :param intercept: bool to return the intercept variable or not

    :return: tuple (pred_expr correlation, intercept variable; None if intercept=False)
    """
    wcov = mdot([wmat.T, ldmat, wmat])
    scale = np.diag(1 / np.sqrt(np.diag(wcov)))
    wcor = mdot([scale, wcov, scale])

    if intercept:
        inter = mdot([scale, wmat.T, ldmat])
        return wcor, inter
    else:
        return wcor, None 

def unOrthoDen(self):
        """Routine to unorthogonalize the orthonormal Density matrix to AO basis"""
        self.P = np.dot(self.X,np.dot(self.PO,self.X.T)) 

def test_vector_as_last_argument(self):
        # The last argument can be 1-D
        A = np.random.random((6, 2))
        B = np.random.random((2, 6))
        C = np.random.random((6, 2))
        D1d = np.random.random(2)  # 1-D

        # the result should be 1-D
        assert_equal(multi_dot([A, B, C, D1d]).shape, (6,)) 

def test_too_few_input_arrays(self):
        assert_raises(ValueError, multi_dot, [])
        assert_raises(ValueError, multi_dot, [np.random.random((3, 3))]) 

def orthoFock(self):
        """Routine to orthogonalize the AO Fock matrix to orthonormal basis"""
        self.FO = np.dot(self.X.T,np.dot(self.F,self.X)) 

def unOrthoFock(self):
        """Routine to unorthogonalize the orthonormal Fock matrix to AO basis"""
        self.F = np.dot(self.U.T,np.dot(self.FO,self.U)) 

def orthoDen(self):
        """Routine to orthogonalize the AO Density matrix to orthonormal basis"""
        self.PO = np.dot(self.U,np.dot(self.P,self.U.T)) 

def updateDIIS(self,F,P):
        FPS =   dot([F,P,self.S])
        SPF =   self.adj(FPS) 
        # error must be in orthonormal basis
        error = dot([self.X,FPS-SPF,self.X]) 
        self.fockSet.append(self.F)
        self.errorSet.append(error) 
        numFock = len(self.fockSet)
        # limit subspace, hardcoded for now
        if numFock > 8:
            del self.fockSet[0] 
            del self.errorSet[0] 
            numFock -= 1
        B = np.zeros((numFock + 1,numFock + 1)) 
        B[-1,:] = B[:,-1] = -1.0
        B[-1,-1] = 0.0
        # B is symmetric
        for i in range(numFock):
            for j in range(i+1):
                B[i,j] = B[j,i] = \
                    np.real(np.trace(np.dot(self.adj(self.errorSet[i]),
                                                     self.errorSet[j])))
        residual = np.zeros(numFock + 1)
        residual[-1] = -1.0
        weights = np.linalg.solve(B,residual)

        # weights is 1 x numFock + 1, but first numFock values
        # should sum to one if we are doing DIIS correctly
        assert np.isclose(sum(weights[:-1]),1.0)

        F = np.zeros((self.nbasis,self.nbasis),dtype='complex')
        for i, Fock in enumerate(self.fockSet):
            F += weights[i] * Fock

        return F 

def test_basic_function_with_dynamic_programing_optimization(self):
        # multi_dot with four or more arguments uses the dynamic programing
        # optimization and therefore deserve a separate
        A = np.random.random((6, 2))
        B = np.random.random((2, 6))
        C = np.random.random((6, 2))
        D = np.random.random((2, 1))
        assert_almost_equal(multi_dot([A, B, C, D]), A.dot(B).dot(C).dot(D)) 

def test_basic_function_with_three_arguments(self):
        # multi_dot with three arguments uses a fast hand coded algorithm to
        # determine the optimal order. Therefore test it separately.
        A = np.random.random((6, 2))
        B = np.random.random((2, 6))
        C = np.random.random((6, 2))

        assert_almost_equal(multi_dot([A, B, C]), A.dot(B).dot(C))
        assert_almost_equal(multi_dot([A, B, C]), np.dot(A, np.dot(B, C))) 

def test_basic_function_with_three_arguments(self):
        # multi_dot with three arguments uses a fast hand coded algorithm to
        # determine the optimal order. Therefore test it separately.
        A = np.random.random((6, 2))
        B = np.random.random((2, 6))
        C = np.random.random((6, 2))

        assert_almost_equal(multi_dot([A, B, C]), A.dot(B).dot(C))
        assert_almost_equal(multi_dot([A, B, C]), np.dot(A, np.dot(B, C))) 

def __init__(self):
        matrix_pairs = {}
        rng = np.random.RandomState(SEED)
        for s in SIZES:
            for n_terms in [2, 5]:
                arrays = [
                    rng.standard_normal((s, s)) for _ in range(n_terms)]
                matrix = matrices.SquareMatrixProduct([
                    matrices.DenseSquareMatrix(a) for a in arrays])
                matrix_pairs[(s, n_terms)] = (matrix, nla.multi_dot(arrays))
        super().__init__(matrix_pairs, rng) 

def get_resid(zscores, swld, wcor):
    """
    Regress out the average pleiotropic signal tagged by TWAS at the region

    :param zscores: numpy.ndarray TWAS zscores
    :param swld: numpy.ndarray intercept variable
    :param wcor: numpy.ndarray predicted expression correlation

    :return: tuple (residual TWAS zscores, intercept z-score)
    """
    m, m = wcor.shape
    m, p = swld.shape

    # create mean factor
    intercept = swld.dot(np.ones(p))

    # estimate under the null for variance components, i.e. V = SW LD SW
    wcor_inv, rank = lin.pinvh(wcor, return_rank=True)

    numer = mdot([intercept.T, wcor_inv, zscores])
    denom = mdot([intercept.T, wcor_inv, intercept])
    alpha = numer / denom
    resid = zscores - intercept * alpha

    s2 = mdot([resid, wcor_inv, resid]) / (rank - 1)
    inter_se = np.sqrt(s2 / denom)
    inter_z = alpha / inter_se

    return resid, inter_z 

def test_vector_as_first_and_last_argument(self):
        # The first and last arguments can be 1-D
        A1d = np.random.random(2)  # 1-D
        B = np.random.random((2, 6))
        C = np.random.random((6, 2))
        D1d = np.random.random(2)  # 1-D

        # the result should be a scalar
        assert_equal(multi_dot([A1d, B, C, D1d]).shape, ()) 

def test_vector_as_last_argument(self):
        # The last argument can be 1-D
        A = np.random.random((6, 2))
        B = np.random.random((2, 6))
        C = np.random.random((6, 2))
        D1d = np.random.random(2)  # 1-D

        # the result should be 1-D
        assert_equal(multi_dot([A, B, C, D1d]).shape, (6,)) 

def test_vector_as_first_argument(self):
        # The first argument can be 1-D
        A1d = np.random.random(2)  # 1-D
        B = np.random.random((2, 6))
        C = np.random.random((6, 2))
        D = np.random.random((2, 2))

        # the result should be 1-D
        assert_equal(multi_dot([A1d, B, C, D]).shape, (2,)) 

def test_basic_function_with_dynamic_programing_optimization(self):
        # multi_dot with four or more arguments uses the dynamic programing
        # optimization and therefore deserve a separate
        A = np.random.random((6, 2))
        B = np.random.random((2, 6))
        C = np.random.random((6, 2))
        D = np.random.random((2, 1))
        assert_almost_equal(multi_dot([A, B, C, D]), A.dot(B).dot(C).dot(D)) 

def test_basic_function_with_three_arguments(self):
        # multi_dot with three arguments uses a fast hand coded algorithm to
        # determine the optimal order. Therefore test it separately.
        A = np.random.random((6, 2))
        B = np.random.random((2, 6))
        C = np.random.random((6, 2))

        assert_almost_equal(multi_dot([A, B, C]), A.dot(B).dot(C))
        assert_almost_equal(multi_dot([A, B, C]), np.dot(A, np.dot(B, C))) 

def test_vector_as_first_and_last_argument(self):
        # The first and last arguments can be 1-D
        A1d = np.random.random(2)  # 1-D
        B = np.random.random((2, 6))
        C = np.random.random((6, 2))
        D1d = np.random.random(2)  # 1-D

        # the result should be a scalar
        assert_equal(multi_dot([A1d, B, C, D1d]).shape, ()) 

def test_vector_as_last_argument(self):
        # The last argument can be 1-D
        A = np.random.random((6, 2))
        B = np.random.random((2, 6))
        C = np.random.random((6, 2))
        D1d = np.random.random(2)  # 1-D

        # the result should be 1-D
        assert_equal(multi_dot([A, B, C, D1d]).shape, (6,)) 

def test_vector_as_first_argument(self):
        # The first argument can be 1-D
        A1d = np.random.random(2)  # 1-D
        B = np.random.random((2, 6))
        C = np.random.random((6, 2))
        D = np.random.random((2, 2))

        # the result should be 1-D
        assert_equal(multi_dot([A1d, B, C, D]).shape, (2,)) 

def test_basic_function_with_dynamic_programing_optimization(self):
        # multi_dot with four or more arguments uses the dynamic programing
        # optimization and therefore deserve a separate
        A = np.random.random((6, 2))
        B = np.random.random((2, 6))
        C = np.random.random((6, 2))
        D = np.random.random((2, 1))
        assert_almost_equal(multi_dot([A, B, C, D]), A.dot(B).dot(C).dot(D))