Python scipy.linalg.schur() Examples

def matrix_sign(M):
    """ The "sign" matrix of `M` """
    #Notes: sign(M) defined s.t. eigvecs of sign(M) are evecs of M
    # and evals of sign(M) are +/-1 or 0 based on sign of eigenvalues of M

    #Using the extremely numerically stable (but expensive) Schur method
    # see http://www.maths.manchester.ac.uk/~higham/fm/OT104HighamChapter5.pdf
    N = M.shape[0]; assert(M.shape == (N, N)), "M must be square!"
    T, Z = _spl.schur(M, 'complex')  # M = Z T Z^H where Z is unitary and T is upper-triangular
    U = _np.zeros(T.shape, 'complex')  # will be sign(T), which is easy to compute
    # (U is also upper triangular), and then sign(M) = Z U Z^H

    # diagonals are easy
    U[_np.diag_indices_from(U)] = _np.sign(_np.diagonal(T))

    #Off diagonals: use U^2 = I or TU = UT
    # Note: Tij = Uij = 0 when i > j and i==j easy so just consider i<j case
    # 0 = sum_k Uik Ukj =  (i!=j b/c off-diag)
    # FUTURE: speed this up by using np.dot instead of sums below
    for j in range(1, N):
        for i in range(j - 1, -1, -1):
            S = U[i, i] + U[j, j]
            if _np.isclose(S, 0):  # then use TU = UT
                if _np.isclose(T[i, i] - T[j, j], 0):  # then just set to zero
                    U[i, j] = 0.0  # TODO: check correctness of this case
                else:
                    U[i, j] = T[i, j] * (U[i, i] - U[j, j]) / (T[i, i] - T[j, j]) + \
                        sum([U[i, k] * T[k, j] - T[i, k] * U[k, j] for k in range(i + 1, j)]) \
                        / (T[i, i] - T[j, j])
            else:  # use U^2 = I
                U[i, j] = - sum([U[i, k] * U[k, j] for k in range(i + 1, j)]) / S
    return _np.dot(Z, _np.dot(U, _np.conjugate(Z.T)))

    #Quick & dirty - not always stable:
    #U,_,Vt = _np.linalg.svd(M)
    #return _np.dot(U,Vt) 

def power(self, exponent: float):
        """Creates a unitary gate as `gate^exponent`.

        Args:
            exponent (float): Gate^exponent

        Returns:
            qiskit.extensions.UnitaryGate: To which `to_matrix` is self.to_matrix^exponent.

        Raises:
            CircuitError: If Gate is not unitary
        """
        from qiskit.quantum_info.operators import Operator  # pylint: disable=cyclic-import
        from qiskit.extensions.unitary import UnitaryGate  # pylint: disable=cyclic-import
        # Should be diagonalized because it's a unitary.
        decomposition, unitary = schur(Operator(self).data, output='complex')
        # Raise the diagonal entries to the specified power
        decomposition_power = list()

        decomposition_diagonal = decomposition.diagonal()
        # assert off-diagonal are 0
        if not np.allclose(np.diag(decomposition_diagonal), decomposition):
            raise CircuitError('The matrix is not diagonal')

        for element in decomposition_diagonal:
            decomposition_power.append(pow(element, exponent))
        # Then reconstruct the resulting gate.
        unitary_power = unitary @ np.diag(decomposition_power) @ unitary.conj().T
        return UnitaryGate(unitary_power, label='%s^%s' % (self.name, exponent)) 

def _choi_to_kraus(data, input_dim, output_dim, atol=ATOL_DEFAULT):
    """Transform Choi representation to Kraus representation."""
    # Check if hermitian matrix
    if is_hermitian_matrix(data, atol=atol):
        # Get eigen-decomposition of Choi-matrix
        # This should be a call to la.eigh, but there is an OpenBlas
        # threading issue that is causing segfaults.
        # Need schur here since la.eig does not
        # guarentee orthogonality in degenerate subspaces
        w, v = la.schur(data, output='complex')
        w = w.diagonal().real
        # Check eigenvalues are non-negative
        if len(w[w < -atol]) == 0:
            # CP-map Kraus representation
            kraus = []
            for val, vec in zip(w, v.T):
                if abs(val) > atol:
                    k = np.sqrt(val) * vec.reshape(
                        (output_dim, input_dim), order='F')
                    kraus.append(k)
            # If we are converting a zero matrix, we need to return a Kraus set
            # with a single zero-element Kraus matrix
            if not kraus:
                kraus.append(np.zeros((output_dim, input_dim), dtype=complex))
            return kraus, None
    # Non-CP-map generalized Kraus representation
    mat_u, svals, mat_vh = la.svd(data)
    kraus_l = []
    kraus_r = []
    for val, vec_l, vec_r in zip(svals, mat_u.T, mat_vh.conj()):
        kraus_l.append(
            np.sqrt(val) * vec_l.reshape((output_dim, input_dim), order='F'))
        kraus_r.append(
            np.sqrt(val) * vec_r.reshape((output_dim, input_dim), order='F'))
    return kraus_l, kraus_r 

def get_cholesky_like_decomposition(mat: np.array) -> np.array:
    """Given a PSD matrix A, finds a matrix T such that TT^{dagger}
    is an approximation of A
    Args:
        mat: A nxn matrix, assumed to be positive semidefinite.
    Returns:
        A matrix T such that TT^{dagger} approximates A
    """
    decomposition, unitary = schur(mat, output='complex')
    eigenvals = np.array(decomposition.diagonal())
    # if a 0 eigenvalue is represented by infinitisimal negative float
    eigenvals[eigenvals < 0] = 0
    DD = np.diag(np.sqrt(eigenvals))
    return unitary @ DD 

def schur_ordered(A, ct=False):
    r"""Returns block ordered complex Schur form of matrix :math:`\mathbf{A}`

    .. math:: \mathbf{TAT}^H = \mathbf{A}_s = \begin{bmatrix} A_{11} & A_{12} \\ 0 & A_{22} \end{bmatrix}

    where :math:`A_{11}\in\mathbb{C}^{s\times s}` contains the :math:`s` stable
    eigenvalues of :math:`\mathbf{A}\in\mathbb{R}^{m\times m}`.

    Args:
        A (np.ndarray): Matrix to decompose.
        ct (bool): Continuous time system.

    Returns:
        tuple: Tuple containing the Schur decomposition of :math:`\mathbf{A}`, :math:`\mathbf{A}_s`; the transformation
        :math:`\mathbf{T}\in\mathbb{C}^{m\times m}`; and the number of stable eigenvalues of :math:`\mathbf{A}`.

    Notes:
        This function is a wrapper of ``scipy.linalg.schur`` imposing the settings required for this application.

    """
    if ct:
        sort_eigvals = 'lhp'
    else:
        sort_eigvals = 'iuc'

    # if A.dtype == complex:
    #     output_form = 'complex'
    # else:
    #     output_form = 'real'
    # issues when not using the complex form of the Schur decomposition

    output_form = 'complex'
    As, Tt, n_stable1 = sclalg.schur(A, output=output_form, sort=sort_eigvals)

    if sort_eigvals == 'lhp':
        n_stable = np.sum(np.linalg.eigvals(A).real <= 0)
    elif sort_eigvals == 'iuc':
        n_stable = np.sum(np.abs(np.linalg.eigvals(A)) <= 1.)
    else:
        raise NameError('Unknown sorting of eigenvalues. Either iuc or lhp')

    assert n_stable == n_stable1, 'Number of stable eigenvalues not equal in Schur output and manual calculation'

    assert (np.abs(As-np.conj(Tt.T).dot(A.dot(Tt))) < 1e-4).all(), 'Schur breakdown - A_schur != T^H A T'
    return As, Tt.T, n_stable