Python scipy.linalg.qr() Examples

def __init__(self,in_channel):
        super(InvConv,self).__init__()

        weight=np.random.randn(in_channel,in_channel)
        q,_=linalg.qr(weight)
        w_p,w_l,w_u=linalg.lu(q.astype(np.float32))
        w_s=np.diag(w_u)
        w_u=np.triu(w_u,1)
        u_mask=np.triu(np.ones_like(w_u),1)
        l_mask=u_mask.T

        self.register_buffer('w_p',torch.from_numpy(w_p))
        self.register_buffer('u_mask',torch.from_numpy(u_mask))
        self.register_buffer('l_mask',torch.from_numpy(l_mask))
        self.register_buffer('l_eye',torch.eye(l_mask.shape[0]))
        self.register_buffer('s_sign',torch.sign(torch.from_numpy(w_s)))
        self.w_l=torch.nn.Parameter(torch.from_numpy(w_l))
        self.w_s=torch.nn.Parameter(torch.log(1e-7+torch.abs(torch.from_numpy(w_s))))
        self.w_u=torch.nn.Parameter(torch.from_numpy(w_u))

        self.weight=None
        self.invweight=None

        return 

def simple_cnot_circuit(num_qubits, measure=True):
    """Creates a simple circuit composed by cnot gates, with measurements or
    not at the end of each qubit.

    Args:
        num_qubits (int): Number of qubits
        measure (bool): Add measurements at the end of each qubit

    Returns:
        QuantumCircuit: The simple quantum circuit
    """
    qr = QuantumRegister(num_qubits)
    circuit = QuantumCircuit(qr)
    for i in range(num_qubits):
        # for the last qubit, we exchange control and target qubits
        target_qubit = i + 1 if num_qubits - 1 > i else i - 1
        circuit.cx(qr[i], qr[target_qubit])

    if measure:
        circuit = _add_measurements(circuit, qr)
    return circuit


# pylint: disable=invalid-name 

def simple_u3_circuit(num_qubits, measure=True):
    """Creates a simple circuit composed by u3 gates, with measurements or not
    at the end of each qubit.

    Args:
        num_qubits (int): Number of qubits
        measure (bool): Add measurements at the end of each qubit

    Returns:
        QuantumCircuit: The simple quantum circuit
    """
    qr = QuantumRegister(num_qubits)
    circuit = QuantumCircuit(qr)
    for i in range(num_qubits):
        circuit.u3(1.1, 2.2, 3.3, qr[i])

    if measure:
        circuit = _add_measurements(circuit, qr)
    return circuit 

def qft_circuit(num_qubits, measure=True):
    """Create a qft circuit.

    Args:
        num_qubits (int): number of qubits
        measure (bool): include measurement in circuit.

    Returns:
        QftCircuit: A qft circuit.
    """
    # Create quantum/classical registers of size n
    qr = QuantumRegister(num_qubits)
    circuit = QuantumCircuit(qr)

    for i in range(num_qubits):
        for j in range(i):
            circuit.cu1(math.pi/float(2**(i-j)), qr[i], qr[j])
        circuit.h(qr[i])

    if measure is True:
        circuit = _add_measurements(circuit, qr)
    return circuit 

def starting_point(self, random=False):
        """Heuristic to find a starting point candidate

        Parameters
        ----------
        random : `bool`
            Use a random orthogonal matrix instead of identity

        Returns
        -------
        startint_point : `np.ndarray`, shape=(n_nodes, n_nodes)
            A starting point candidate
        """
        sqrt_C = sqrtm(self.covariance)
        sqrt_L = np.sqrt(self.mean_intensity)
        if random:
            random_matrix = np.random.rand(self.n_nodes, self.n_nodes)
            M, _ = qr(random_matrix)
        else:
            M = np.eye(self.n_nodes)
        initial = np.dot(np.dot(sqrt_C, M), np.diag(1. / sqrt_L))
        return initial 

def _absorb_constraints(design_matrix, constraints):
    """Absorb model parameters constraints into the design matrix.

    :param design_matrix: The (2-d array) initial design matrix.
    :param constraints: The 2-d array defining initial model parameters
     (``betas``) constraints (``np.dot(constraints, betas) = 0``).
    :return: The new design matrix with absorbed parameters constraints.

    :raise ImportError: if scipy is not found, used for ``scipy.linalg.qr()``
      which is cleaner than numpy's version requiring a call like
      ``qr(..., mode='complete')`` to get a full QR decomposition.
    """
    try:
        from scipy import linalg
    except ImportError: # pragma: no cover
        raise ImportError("Cubic spline functionality requires scipy.")

    m = constraints.shape[0]
    q, r = linalg.qr(np.transpose(constraints))

    return np.dot(design_matrix, q[:, m:]) 

def predict(self, u=0):
        """
        Predict next state (prior) using the Kalman filter state propagation
        equations.

        Parameters
        ----------

        u : np.array, optional
            Optional control vector. If non-zero, it is multiplied by B
            to create the control input into the system.
        """

        # x = Fx + Bu
        self.x = dot(self.F, self.x) + dot(self.B, u)

        # P = FPF' + Q
        _, P2 = qr(np.hstack([dot(self.F, self._P1_2), self._Q1_2]).T)
        self._P1_2 = P2[:self.dim_x, :self.dim_x].T

        # copy prior
        self.x_prior = np.copy(self.x)
        self._P1_2_prior = np.copy(self._P1_2) 

def haar_measure(n):
    """A Random matrix distributed with the Haar measure.

    For more details, see :cite:`mezzadri2006`.

    Args:
        n (int): matrix size
    Returns:
        array: an nxn random matrix
    """
    z = (sp.randn(n, n) + 1j * sp.randn(n, n)) / np.sqrt(2.0)
    q, r = qr(z)
    d = sp.diagonal(r)
    ph = d / np.abs(d)
    q = np.multiply(q, ph, q)
    return q 

def haar_measure(n):
    """A Random matrix distributed with the Haar measure.

    For more details, see :cite:`mezzadri2006`.

    Args:
        n (int): matrix size
    Returns:
        array: an nxn random matrix
    """
    z = (sp.randn(n, n) + 1j * sp.randn(n, n)) / np.sqrt(2.0)
    q, r = qr(z)
    d = sp.diagonal(r)
    ph = d / np.abs(d)
    q = np.multiply(q, ph, q)
    return q 

def test_delete_1x1_row_col(self):
        a, q, r = self.generate('1x1')
        q1, r1 = qr_delete(q, r, 0, 1, 'row')
        assert_equal(q1, np.ndarray(shape=(0, 0), dtype=q.dtype))
        assert_equal(r1, np.ndarray(shape=(0, r.shape[1]), dtype=r.dtype))

        a, q, r = self.generate('1x1')
        q1, r1 = qr_delete(q, r, 0, 1, 'col')
        assert_unitary(q1)
        assert_(q1.dtype == q.dtype)
        assert_(q1.shape == q.shape)
        assert_equal(r1, np.ndarray(shape=(r.shape[0], 0), dtype=r.dtype))

    # all full qr, row deletes and single column deletes should be able to
    # handle any non negative strides. (only row and column vector
    # operations are used.) p column delete require fortran ordered
    # Q and R and will make a copy as necessary.  Economic qr row deletes 
    # requre a contigous q. 

def __init__(self, in_channel):
        super().__init__()

        weight = np.random.randn(in_channel, in_channel)
        q, _ = la.qr(weight)
        w_p, w_l, w_u = la.lu(q.astype(np.float32))
        w_s = np.diag(w_u)
        w_u = np.triu(w_u, 1)
        u_mask = np.triu(np.ones_like(w_u), 1)
        l_mask = u_mask.T

        w_p = torch.from_numpy(w_p)
        w_l = torch.from_numpy(w_l)
        w_s = torch.from_numpy(w_s)
        w_u = torch.from_numpy(w_u)

        self.register_buffer('w_p', w_p)
        self.register_buffer('u_mask', torch.from_numpy(u_mask))
        self.register_buffer('l_mask', torch.from_numpy(l_mask))
        self.register_buffer('s_sign', torch.sign(w_s))
        self.register_buffer('l_eye', torch.eye(l_mask.shape[0]))
        self.w_l = nn.Parameter(w_l)
        self.w_s = nn.Parameter(logabs(w_s))
        self.w_u = nn.Parameter(w_u) 

def _absorb_constraints(design_matrix, constraints):
    """Absorb model parameters constraints into the design matrix.

    :param design_matrix: The (2-d array) initial design matrix.
    :param constraints: The 2-d array defining initial model parameters
     (``betas``) constraints (``np.dot(constraints, betas) = 0``).
    :return: The new design matrix with absorbed parameters constraints.

    :raise ImportError: if scipy is not found, used for ``scipy.linalg.qr()``
      which is cleaner than numpy's version requiring a call like
      ``qr(..., mode='complete')`` to get a full QR decomposition.
    """
    try:
        from scipy import linalg
    except ImportError: # pragma: no cover
        raise ImportError("Cubic spline functionality requires scipy.")

    m = constraints.shape[0]
    q, r = linalg.qr(np.transpose(constraints))

    return np.dot(design_matrix, q[:, m:]) 

def null_space_method(Ao, bo, Co, do, verbose):
    A = deepcopy(Ao)
    C = deepcopy(Co)
    b = deepcopy(bo)
    d = deepcopy(do)
    m, n = A.shape
    p, n = C.shape
    Q, R = np.linalg.qr(C.T, 'complete')
    Q1 = Q[0:n, 0:p]
    Q2 = Q[0:n, p:n]
    # Lower triangular matrix!
    L = R.T
    L = L[0:p, 0:p]
    y1 = least_squares(L, d, verbose)
    c = b - np.dot( np.dot(A , Q1) , y1)
    AQ2 = np.dot(A , Q2)
    y2 = least_squares(AQ2 , c, verbose)
    x = np.dot(Q1 , y1) + np.dot(Q2 , y2)
    cond = np.linalg.cond(AQ2)
    if verbose is True:
        print('The condition number of the matrix is '+str(cond)+'.')
    return x 

def minimum_norm(A, b):
    Q, R, pvec = qr(A, pivoting=True)
    m, n = A.shape
    r = np.linalg.matrix_rank(A)
    Q1 = Q[0:r, 0:r]
    R1 = R[0:r, 0:r]
    indices = np.argsort(pvec)
    P = np.eye(n)
    temp = P[indices,:]
    P1 = temp[0:n, 0:r]
    x = np.dot(P1 ,  np.dot( np.linalg.inv(R1)  , np.dot( Q1.T , b ) ) )
    x = x.reshape(n, 1)
    return x 

def orthonormalize(A, overwrite_a=True, check_finite=False):
    """orthonormalize the columns of A via QR decomposition"""
    # NOTE: for A(m, n) 'economic' returns Q(m, k), R(k, n) where k is min(m, n)
    # TODO: when does overwrite_a even work? (fortran?)
    Q, _ = linalg.qr(A, overwrite_a=overwrite_a, check_finite=check_finite,
                     mode='economic', pivoting=False)
    return Q 

def _remove_rank_def_rows(b_mat, nvelpar, ifgv, sel):
    """
    Remove rank deficient rows of design matrix
    """
    _, _, e_var = qr(b_mat, mode='economic', pivoting=True)
    licols = e_var[matrix_rank(b_mat):nvelpar]
    [rmrow, _] = where(b_mat[:, licols] != 0)
    b_mat = delete(b_mat, rmrow, axis=0)
    ifgv = delete(ifgv, rmrow)
    sel = delete(sel, rmrow)
    return b_mat, ifgv, sel, rmrow 

def get_svd_subset_selection(Ao, number_of_subsamples):
    """
    Singular value decomposition and pivoted QR factorization, where the pivots
    are used as a heuristic for subsampling.
    """
    A = deepcopy(Ao)
    _, _, V = svd(A.T)
    _, _, pvec = qr(V[:, 0:number_of_subsamples].T , pivoting=True )
    z = pvec[0:number_of_subsamples]
    return z 

def test_crs_with_specific_constraint():
    from patsy.highlevel import incr_dbuilder, build_design_matrices, dmatrix
    x = (-1.5)**np.arange(20)
    # Hard coded R values for smooth: s(x, bs="cr", k=5)
    # R> knots <- smooth$xp
    knots_R = np.array([-2216.837820053100585937,
                        -50.456909179687500000,
                        -0.250000000000000000,
                        33.637939453125000000,
                        1477.891880035400390625])
    # R> centering.constraint <- t(qr.X(attr(smooth, "qrc")))
    centering_constraint_R = np.array([[0.064910676323168478574,
                                        1.4519875239407085132,
                                        -2.1947446912471946234,
                                        1.6129783104357671153,
                                        0.064868180547550072235]])
    # values for which we want a prediction
    new_x = np.array([-3000., -200., 300., 2000.])
    result1 = dmatrix("cr(new_x, knots=knots_R[1:-1], "
                      "lower_bound=knots_R[0], upper_bound=knots_R[-1], "
                      "constraints=centering_constraint_R)")

    data_chunked = [{"x": x[:10]}, {"x": x[10:]}]
    new_data = {"x": new_x}
    builder = incr_dbuilder("cr(x, df=4, constraints='center')",
                            lambda: iter(data_chunked))
    result2 = build_design_matrices([builder], new_data)[0]

    assert np.allclose(result1, result2, rtol=1e-12, atol=0.) 

def potential(g, div_neg=None):
    """potential is related to the intrinsic time. Note that the returned value from this function is the negative of
    potential. Thus small potential is related to smaller intrinsic time and vice versa."""

    div_neg = -div(g) if div_neg is None else div_neg
    g_undirected = g.copy()
    g_undirected.to_undirected()
    L = np.array(g_undirected.laplacian())
    Q, R = qr(L)
    p = np.linalg.pinv(R).dot(Q.T).dot(div_neg)

    res = p - p.min()
    return res 

def __hocur_find_li_cols(matrix, tol=1e-14):
    """Find linearly independent columns of a matrix.

    Parameters
    ----------
    matrix: np.ndarray (m,n)
        rectangular matrix

    Returns
    -------
    cols: list[int]
        indices of linearly independent columns
    """

    # define column list
    cols = []

    # apply QR decomposition with pivoting
    _, r, p = splin.qr(matrix, pivoting=True, mode='economic')

    if tol == 0:
        cols = [p[i] for i in range(matrix.shape[0])]
    else:
        for i in range(r.shape[0]):
            if np.abs(r[i, i]) > tol:
                cols.append(p[i])

    return cols 

def pivot_columns(a, rank=None, threshold=None):
    """Computes the QR decomposition of a matrix with column pivoting, i.e. solves the equation AP=QR such that Q is
    orthogonal, R is upper triangular, and P is a permutation matrix.

    Args:
        a (np.ndarray):    Matrix for which to compute QR decomposition.
        threshold (float): Threshold specifying approximate rank of a. All singular values less than threshold * (largest singular value) will be set to 0
        rank (int):        The approximate rank.
    Returns:
        np.array: The permutation p.
    """
    assert (threshold is None) != (rank is None), "Exactly one of threshold and rank should be specified."
    if threshold is not None:
        rank = approx_matrix_rank(a, threshold)
    return qr(a, pivoting=True)[2][:rank] 

def ortho_matrix(n, m=None):
    if m is None:
        m = n
    max_nm = max(n, m)
    from scipy.linalg import qr
    H = np.random.randn(max_nm, max_nm) + 1j * np.random.randn(max_nm, max_nm)
    Q, _ = qr(H)
    M = Q.dot(np.conjugate(Q.T))
    return M[:n, :] 

def __init__(self, in_channel):
        super().__init__()

        weight = torch.randn(in_channel, in_channel)
        q, _ = torch.qr(weight)
        weight = q.unsqueeze(2).unsqueeze(3)
        self.weight = nn.Parameter(weight) 

def randomized_range_finder(A, size, n_iter):
    """
    Computes an orthonormal matrix whose range
    approximates the range of A.

    Parameters
    ----------
    A: 2D array
        The input data matrix

    size: integer
        Size of the return array

    n_iter: integer
        Number of power iterations used to stabilize the result

    Returns
    -------
    Q: 2D array
        A (size x size) projection matrix, the range of which
        approximates well the range of the input matrix A.

    Notes
    -----
    scikit-learn implementation
    """
    # Generating normal random vectors with shape: (A.shape[1], size)
    Q = np.random.normal(size=(A.shape[1], size))

    # Perform power iterations with Q to further 'imprint' the top
    # singular vectors of A in Q
    for i in range(n_iter):
        Q, _ = linalg.lu(np.dot(A, Q), permute_l=True)
        Q, _ = linalg.lu(np.dot(A.T, Q), permute_l=True)

    # Sample the range of A using by linear projection of Q
    # Extract an orthonormal basis
    Q, _ = linalg.qr(np.dot(A, Q), mode='economic')
    return Q 

def block_arnoldi_krylov(r, F, G, approx_type='Pade', side='controllability'):

    n = G.shape[0]
    m = G.shape[1]

    Q0, R0, P0 = sclalg.qr(G, pivoting=True)
    Q0 = Q0[:, :m]

    Q = np.zeros((n,m*r), dtype=complex)
    V = np.zeros((n,m*r), dtype=complex)

    for k in range(r):

        if k == 0:
            Q[:, 0:m] = F.dot(Q0)
        else:
            Q[:, k*m: k*m + m] = F.dot(Q[:, (k-1)*m:(k-1)*m + m])

        Q[:, :k*m + m] = mgs_ortho(Q[:, :k*m+m])

        Qf, R, P = sclalg.qr(Q[:, k*m: k*m + m], pivoting=True)
        Q[:, k*m: k*m + m ] = Qf[:, :m]
        if R[0,0] >= 1e-6:
            V[:, k*m:k*m + m] = Q[:, k*m: k*m + m ]
        else:
            print('Deflating')
            k -= 1

    V = mgs_ortho(V)

    return V 

def mgs_ortho(X):
    r"""
    Modified Gram-Schmidt Orthogonalisation

    Orthogonalises input matrix :math:`\mathbf{X}` column by column.

    Args:
        X (np.ndarray): Input matrix of dimensions :math:`n` by :math:`m`.

    Returns:
        np.ndarray: Orthogonalised matrix of dimensions :math:`n` by :math:`m`.

    Notes:
        This method is faster than scipy's :func:`scipy.linalg.qr` method that returns an orthogonal matrix as part of
        the QR decomposition, albeit at a higher number of function calls.
    """

    # Q, R = sclalg.qr(X)
    n = X.shape[1]
    m = X.shape[0]

    Q = np.zeros((m, n), dtype=float)

    for i in range(n):
        w = X[:, i]
        for j in range(i):
            h = Q[:, j].T.dot(w)
            w = w - h * Q[:, j]
        Q[:, i] = w / sclalg.norm(w)

    return Q 

def nullspace_qr(m, tol=1e-7):
    """
    Compute the nullspace of a matrix using the QR decomposition.

    The QR decomposition is faster but less accurate than the SVD
    used by :func:`nullspace`.

    Parameters
    ----------
    m : numpy array
       An matrix of shape (M,N) whose nullspace to compute.

    tol : float (optional)
       Nullspace tolerance, used when comparing diagonal values of R with zero.

    Returns
    -------
    An matrix of shape (M,K) whose columns contain nullspace basis vectors.
    """
    #if M,N = m.shape, and q,r,p = _spl.qr(...)
    # q.shape == (N,N), r.shape = (N,M), p.shape = (M,)
    q, r, _ = _spl.qr(m.T, mode='full', pivoting=True)
    rank = (_np.abs(_np.diagonal(r)) > tol).sum()

    #DEBUG: requires q,r,p = _sql.qr(...) above
    #assert( _np.linalg.norm(_np.dot(q,r) - m.T[:,p]) < 1e-8) #check QR decomp
    #print("Rank QR = ",rank)
    #print('\n'.join(map(str,_np.abs(_np.diagonal(r)))))
    #print("Ret = ", q[:,rank:].shape, " Q = ",q.shape, " R = ",r.shape)

    return q[:, rank:] 

def __init__(self, subsampling_algorithm):
        self.subsampling_algorithm = subsampling_algorithm
        if self.subsampling_algorithm == 'qr':
            self.algorithm = lambda A, k : get_qr_column_pivoting(A, k)
        elif self.subsampling_algorithm == 'svd':
            self.algorithm = lambda A, k : get_svd_subset_selection(A, k)
        elif self.subsampling_algorithm == 'newton':
            self.algorithm = lambda A, k : get_newton_determinant_maximization(A, k)
        elif self.subsampling_algorithm == 'random':
            self.algorithm = 0 #np.random.choice(int(m), m_refined, replace=False)
        elif self.subsampling_algorithm == None:
            self.algorithm = lambda A, k: _get_all_pivots(A, k) 

def _weightsamples(self):
        reps = int(np.ceil(self.n / self.d))
        Q = np.empty((self.d, self.d*reps))

        for r in range(reps):
            W = self._random.randn(self.d, self.d)
            Q[:, (r * self.d):((r + 1) * self.d)] = qr(W)[0]

        S = np.sqrt(self._random.chisquare(df=self.d, size=self.d))
        weights = np.diag(S).dot(Q[:, :self.n])
        return weights 

def test_crs_with_specific_constraint():
    from patsy.highlevel import incr_dbuilder, build_design_matrices, dmatrix
    x = (-1.5)**np.arange(20)
    # Hard coded R values for smooth: s(x, bs="cr", k=5)
    # R> knots <- smooth$xp
    knots_R = np.array([-2216.837820053100585937,
                        -50.456909179687500000,
                        -0.250000000000000000,
                        33.637939453125000000,
                        1477.891880035400390625])
    # R> centering.constraint <- t(qr.X(attr(smooth, "qrc")))
    centering_constraint_R = np.array([[0.064910676323168478574,
                                        1.4519875239407085132,
                                        -2.1947446912471946234,
                                        1.6129783104357671153,
                                        0.064868180547550072235]])
    # values for which we want a prediction
    new_x = np.array([-3000., -200., 300., 2000.])
    result1 = dmatrix("cr(new_x, knots=knots_R[1:-1], "
                      "lower_bound=knots_R[0], upper_bound=knots_R[-1], "
                      "constraints=centering_constraint_R)")

    data_chunked = [{"x": x[:10]}, {"x": x[10:]}]
    new_data = {"x": new_x}
    builder = incr_dbuilder("cr(x, df=4, constraints='center')",
                            lambda: iter(data_chunked))
    result2 = build_design_matrices([builder], new_data)[0]

    assert np.allclose(result1, result2, rtol=1e-12, atol=0.)