Python scipy.sparse.linalg.norm() Examples

def frobeniusnorm(ar):
    """
    Compute the frobenius norm of an array (or matrix),

       sqrt( sum( each_element_of_a^2 ) )

    Parameters
    ----------
    ar : numpy array
        What to compute the frobenius norm of.  Note that ar can be any shape
        or number of dimenions.

    Returns
    -------
    float or complex
        depending on the element type of ar.
    """
    return _np.sqrt(_np.sum(ar**2)) 

def _similarity(self, q_vect: Union[csr_matrix, List]) -> List[float]:
        """Calculates cosine similarity between the user's query and product items.

        Parameters:
            q_cur: user's query

        Returns:
            cos_similarities: lits of similarity scores
        """

        norm = sparse_norm(q_vect) * sparse_norm(self.x_train_features, axis=1)
        cos_similarities = np.array(q_vect.dot(self.x_train_features.T).todense()) / norm

        cos_similarities = cos_similarities[0]
        cos_similarities = np.nan_to_num(cos_similarities)
        return cos_similarities 

def compute_rho_estimate(self):
        # Iterates
        x = self.work.x
        y = self.work.y
        z = self.work.z

        # Problem data
        P = self.work.data.P
        q = self.work.data.q
        A = self.work.data.A

        # Compute normalized residuals
        pri_res = la.norm(A.dot(x) - z, np.inf)
        pri_res /= (np.max([la.norm(A.dot(x), np.inf),
                            la.norm(z, np.inf)]) + 1e-10)
        dua_res = la.norm(P.dot(x) + q + A.T.dot(y), np.inf)
        dua_res /= (np.max([la.norm(A.T.dot(y), np.inf),
                           la.norm(P.dot(x), np.inf),
                           la.norm(q, np.inf)]) + 1e-10)

        # Compute new rho
        new_rho = self.work.settings.rho * np.sqrt(pri_res/(dua_res + 1e-10))
        return min(max(new_rho, RHO_MIN), RHO_MAX) 

def compute_dua_tol(self, eps_abs, eps_rel):
        """
        Compute dual tolerance
        """
        P = self.work.data.P
        q = self.work.data.q
        A = self.work.data.A
        if self.work.settings.scaling and not \
                self.work.settings.scaled_termination:
            cinv = self.work.scaling.cinv
            Dinv = self.work.scaling.Dinv
            max_rel_eps = cinv * np.max([
                la.norm(Dinv.dot(A.T.dot(self.work.y)), np.inf),
                la.norm(Dinv.dot(P.dot(self.work.x)), np.inf),
                la.norm(Dinv.dot(q), np.inf)])
        else:
            max_rel_eps = np.max([
                la.norm(A.T.dot(self.work.y), np.inf),
                la.norm(P.dot(self.work.x), np.inf),
                la.norm(q, np.inf)])

        eps_dua = eps_abs + eps_rel * max_rel_eps

        return eps_dua 

def compute_pri_tol(self, eps_abs, eps_rel):
        """
        Compute primal tolerance using problem data
        """
        A = self.work.data.A
        if self.work.settings.scaling and not \
                self.work.settings.scaled_termination:
            Einv = self.work.scaling.Einv
            max_rel_eps = np.max([
                la.norm(Einv.dot(A.dot(self.work.x)), np.inf),
                la.norm(Einv.dot(self.work.z), np.inf)])
        else:
            max_rel_eps = np.max([
                la.norm(A.dot(self.work.x), np.inf),
                la.norm(self.work.z, np.inf)])

        eps_pri = eps_abs + eps_rel * max_rel_eps

        return eps_pri 

def tracenorm(A):
    """
    Compute the trace norm of matrix A given by:

      Tr( sqrt{ A^dagger * A } )

    Parameters
    ----------
    A : numpy array
        The matrix to compute the trace norm of.
    """
    if _np.linalg.norm(A - _np.conjugate(A.T)) < 1e-8:
        #Hermitian, so just sum eigenvalue magnitudes
        return _np.sum(_np.abs(_np.linalg.eigvals(A)))
    else:
        #Sum of singular values (positive by construction)
        return _np.sum(_np.linalg.svd(A, compute_uv=False)) 

def safe_onenorm(A):
    """
    Computes the 1-norm of the dense or sparse matrix `A`.

    Parameters
    ----------
    A : ndarray or sparse matrix
        The matrix or vector to take the norm of.

    Returns
    -------
    float
    """
    if _sps.isspmatrix(A):
        return sparse_onenorm(A)
    else:
        return _np.linalg.norm(A, 1) 

def norm1to1(operator, n_samples=10000, mxBasis="gm", return_list=False):
    """
    Returns the Hermitian 1-to-1 norm of a superoperator represented in
    the standard basis, calculated via Monte-Carlo sampling. Definition
    of Hermitian 1-to-1 norm can be found in arxiv:1109.6887.
    """
    if mxBasis == 'gm':
        std_operator = change_basis(operator, 'gm', 'std')
    elif mxBasis == 'pp':
        std_operator = change_basis(operator, 'pp', 'std')
    elif mxBasis == 'std':
        std_operator = operator
    else:
        raise ValueError("mxBasis should be 'gm', 'pp' or 'std'!")

    rand_dim = int(_np.sqrt(float(len(std_operator))))
    vals = [norm1(unvec(_np.dot(std_operator, vec(random_hermitian(rand_dim)))))
            for n in range(n_samples)]
    if return_list:
        return vals
    else:
        return max(vals)


## ------------------------ General utility fns ----------------------------------- 

def frobeniusnorm2(ar):
    """
    Compute the squared frobenius norm of an array (or matrix),

       sum( each_element_of_a^2 ) )

    Parameters
    ----------
    ar : numpy array
        What to compute the squared frobenius norm of.  Note that ar can be any
        shape or number of dimenions.

    Returns
    -------
    float or complex
        depending on the element type of ar.
    """
    return _np.sum(ar**2) 

def test_matrix_norm(self):

        # Frobenius norm is the default
        assert_allclose(spnorm(self.b), 7.745966692414834)        
        assert_allclose(spnorm(self.b, 'fro'), 7.745966692414834)

        assert_allclose(spnorm(self.b, np.inf), 9)
        assert_allclose(spnorm(self.b, -np.inf), 2)
        assert_allclose(spnorm(self.b, 1), 7)
        assert_allclose(spnorm(self.b, -1), 6)

        # _multi_svd_norm is not implemented for sparse matrix
        assert_raises(NotImplementedError, spnorm, self.b, 2)
        assert_raises(NotImplementedError, spnorm, self.b, -2) 

def array_eq(a, b, tol=1e-8):
    """Test whether arrays `a` and `b` are equal, i.e. if `norm(a-b) < tol` """
    print(_np.linalg.norm(a - b))
    return _np.linalg.norm(a - b) < tol 

def rescale_matrix(matrix, scaling, axis, binary=True, return_scaling_values=False):
    '''Function to scale either rows or columns of the sparse rating matrix'''
    scaling_values = None
    if scaling == 1: # no scaling (standard SVD case)
        result = matrix

    if binary:
        norm = np.sqrt(matrix.getnnz(axis=axis)) # compute Euclidean norm as if values are binary
    else:
        norm = spnorm(matrix, axis=axis, ord=2) # compute Euclidean norm

    scaling_values = power(norm, scaling-1, where=norm != 0)
    scaling_matrix = diags(scaling_values)

    if axis == 0: # scale columns
        result = matrix.dot(scaling_matrix)
    if axis == 1: # scale rows
        result = scaling_matrix.dot(matrix)

    if return_scaling_values:
        result = (result, scaling_values)
    return result


# matvec implementation is based on
# http://stackoverflow.com/questions/18595981/improving-performance-of-multiplication-of-scipy-sparse-matrices 

def test_norm(self, count_matrix_dataset, vectorizer):
        """Normalization test."""
        _, x, y = count_matrix_dataset
        v = vectorizer(norm="l2")
        v.fit(x, y)
        matrix = v.transform(x)
        norm = sp_norm(matrix, axis=1)
        np.testing.assert_allclose(np.mean(norm), 1.) 

def test_row_norms(self):
    A = sp.rand(100, 100, 0.1)
    A = A.tocsc()

    n1 = spla.norm(A, axis=1, ord=2)
    n2 = rbf.sputils.row_norms(A, order=2)
    self.assertTrue(np.allclose(n1, n2))

    n1 = spla.norm(A, axis=1, ord=1)
    n2 = rbf.sputils.row_norms(A, order=1)
    self.assertTrue(np.allclose(n1, n2)) 

def test_uu_implicit():
    "Train and use user-user on an implicit data set."
    algo = knn.UserUser(20, center=False, aggregate='sum')
    data = ml_ratings.loc[:, ['user', 'item']]

    algo.fit(data)
    assert algo.user_means_ is None

    mat = algo.rating_matrix_.to_scipy()
    norms = spla.norm(mat, 2, 1)
    assert norms == approx(1.0)

    preds = algo.predict_for_user(50, [1, 2, 42])
    assert all(preds[preds.notna()] > 0) 

def _normalize(self, rmat):
        rmat = rmat.to_scipy()
        # compute column norms
        norms = spla.norm(rmat, 2, axis=0)
        # and multiply by a diagonal to normalize columns
        recip_norms = norms.copy()
        is_nz = recip_norms > 0
        recip_norms[is_nz] = np.reciprocal(recip_norms[is_nz])
        norm_mat = rmat @ sps.diags(recip_norms)
        assert norm_mat.shape[1] == rmat.shape[1]
        # and reset NaN
        norm_mat.data[np.isnan(norm_mat.data)] = 0
        _logger.info('[%s] normalized rating matrix columns', self._timer)
        return matrix.CSR.from_scipy(norm_mat, False) 

def is_hermitian(operator):
    """
    Check if matrix or operator is hermitian.

    :param (numpy.ndarray|qutip.Qobj) operator: The operator or matrix to be tested.
    :return: True if the operator is hermitian.
    :rtype: bool
    """
    if isinstance(operator, qt.Qobj):
        return (operator.dag() - operator).norm(FROBENIUS) / operator.norm(FROBENIUS) < EPS
    if isinstance(operator, np.ndarray):
        return np.linalg.norm(operator.T.conj() - operator) / np.linalg.norm(operator) < EPS
    return spnorm(operator.H - operator) / spnorm(operator) < EPS 

def test_sparse_vector_norms(self):
        for sparse_type in self._sparse_types:
            for M in self._test_matrices:
                S = sparse_type(M)
                for axis in (0, 1, -1, -2, (0, ), (1, ), (-1, ), (-2, )):
                    assert_allclose(spnorm(S, axis=axis), npnorm(M, axis=axis))
                    for ord in None, 2, np.inf, -np.inf, 1, 0.5, 0.42:
                        assert_allclose(spnorm(S, ord, axis=axis),
                                        npnorm(M, ord, axis=axis)) 

def test_sparse_matrix_norms(self):
        for sparse_type in self._sparse_types:
            for M in self._test_matrices:
                S = sparse_type(M)
                assert_allclose(spnorm(S), npnorm(M))
                assert_allclose(spnorm(S, 'fro'), npnorm(M, 'fro'))
                assert_allclose(spnorm(S, np.inf), npnorm(M, np.inf))
                assert_allclose(spnorm(S, -np.inf), npnorm(M, -np.inf))
                assert_allclose(spnorm(S, 1), npnorm(M, 1))
                assert_allclose(spnorm(S, -1), npnorm(M, -1)) 

def test_norm_exceptions(self):
        m = self.b
        assert_raises(TypeError, spnorm, m, None, 1.5)
        assert_raises(TypeError, spnorm, m, None, [2])
        assert_raises(ValueError, spnorm, m, None, ())
        assert_raises(ValueError, spnorm, m, None, (0, 1, 2))
        assert_raises(ValueError, spnorm, m, None, (0, 0))
        assert_raises(ValueError, spnorm, m, None, (0, 2))
        assert_raises(ValueError, spnorm, m, None, (-3, 0))
        assert_raises(ValueError, spnorm, m, None, 2)
        assert_raises(ValueError, spnorm, m, None, -3)
        assert_raises(ValueError, spnorm, m, 'plate_of_shrimp', 0)
        assert_raises(ValueError, spnorm, m, 'plate_of_shrimp', (0, 1)) 

def test_vector_norm(self):
        v = [4.5825756949558398, 4.2426406871192848, 4.5825756949558398]
        for m, a in (self.b, 0), (self.b.T, 1):
            for axis in a, (a, ), a-2, (a-2, ):
                assert_allclose(spnorm(m, 1, axis=axis), [7, 6, 7])
                assert_allclose(spnorm(m, np.inf, axis=axis), [4, 3, 4])
                assert_allclose(spnorm(m, axis=axis), v)
                assert_allclose(spnorm(m, ord=2, axis=axis), v)
                assert_allclose(spnorm(m, ord=None, axis=axis), v) 

def test_matrix_norm_axis(self):
        for m, axis in ((self.b, None), (self.b, (0, 1)), (self.b.T, (1, 0))):
            assert_allclose(spnorm(m, axis=axis), 7.745966692414834)        
            assert_allclose(spnorm(m, 'fro', axis=axis), 7.745966692414834)
            assert_allclose(spnorm(m, np.inf, axis=axis), 9)
            assert_allclose(spnorm(m, -np.inf, axis=axis), 2)
            assert_allclose(spnorm(m, 1, axis=axis), 7)
            assert_allclose(spnorm(m, -1, axis=axis), 6) 

def safenorm(A, part=None):
    """
    Returns the frobenius norm of a matrix or vector, `A` when it is either
    a dense array or a sparse matrix.

    Parameters
    ----------
    A : ndarray or sparse matrix
        The matrix or vector to take the norm of.

    part : {None,'real','imag'}
        If not None, return the norm of the real or imaginary
        part of `A`.

    Returns
    -------
    float
    """
    if part == 'real': takepart = _np.real
    elif part == 'imag': takepart = _np.imag
    else: takepart = lambda x: x
    if _sps.issparse(A):
        assert(_sps.isspmatrix_csr(A)), "Non-CSR sparse formats not implemented"
        return _np.linalg.norm(takepart(A.data))
    else:
        return _np.linalg.norm(takepart(A))
    # could also use _spsl.norm(A) 

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 compute_dua_res(self, x, y):
        """
        Compute dual residual ||Px + q + A'y||
        """

        dua_res = self.work.data.P.dot(x) +\
            self.work.data.q + self.work.data.A.T.dot(y)

        if self.work.settings.scaling and not \
                self.work.settings.scaled_termination:
            # Use unscaled residual
            dua_res = self.work.scaling.cinv * \
                self.work.scaling.Dinv.dot(dua_res)

        return la.norm(dua_res, np.inf) 

def unitary_superoperator_matrix_log(M, mxBasis):
    """
    Construct the logarithm of superoperator matrix `M`
    that acts as a unitary on density-matrix space,
    (`M: rho -> U rho Udagger`) so that log(M) can be
    written as the action by Hamiltonian `H`:
    `log(M): rho -> -i[H,rho]`.


    Parameters
    ----------
    M : numpy array
        The superoperator matrix whose logarithm is taken

    mxBasis : {'std', 'gm', 'pp', 'qt'} or Basis object
        The source and destination basis, respectively.  Allowed
        values are Matrix-unit (std), Gell-Mann (gm), Pauli-product (pp),
        and Qutrit (qt) (or a custom basis object).

    Returns
    -------
    numpy array
        A matrix `logM`, of the same shape as `M`, such that `M = exp(logM)`
        and `logM` can be written as the action `rho -> -i[H,rho]`.
    """
    from . import lindbladtools as _lt  # (would create circular imports if at top)
    from . import optools as _gt  # (would create circular imports if at top)

    M_std = change_basis(M, mxBasis, "std")
    evals = _np.linalg.eigvals(M_std)
    assert(_np.allclose(_np.abs(evals), 1.0))  # simple but technically incomplete check for a unitary superop
    # (e.g. could be anti-unitary: diag(1, -1, -1, -1))
    U = _gt.process_mx_to_unitary(M_std)
    H = _spl.logm(U) / -1j  # U = exp(-iH)
    logM_std = _lt.hamiltonian_to_lindbladian(H)  # rho --> -i[H, rho] * sqrt(d)/2
    logM = change_basis(logM_std * (2.0 / _np.sqrt(H.shape[0])), "std", mxBasis)
    assert(_np.linalg.norm(_spl.expm(logM) - M) < 1e-8)  # expensive b/c of expm - could comment for performance
    return logM 

def compute_pri_res(self, x, z):
        """
        Compute primal residual ||Ax - z||
        """

        # Primal residual
        Ax = self.work.data.A.dot(x)
        pri_res = Ax - z

        if self.work.settings.scaling and not \
                self.work.settings.scaled_termination:
            pri_res = self.work.scaling.Einv.dot(pri_res)

        return la.norm(pri_res, np.inf) 

def _norm_KKT_cols(self, P, A):
        """
        Compute the norm of the KKT matrix from P and A
        """

        # First half
        norm_P_cols = spspa.linalg.norm(P, np.inf, axis=0)
        norm_A_cols = spspa.linalg.norm(A, np.inf, axis=0)
        norm_first_half = np.maximum(norm_P_cols, norm_A_cols)

        # Second half (norm cols of A')
        norm_second_half = spspa.linalg.norm(A, np.inf, axis=1)

        return np.hstack((norm_first_half, norm_second_half)) 

def norm1(matr):
    """
    Returns the 1 norm of a matrix
    """
    return float(_np.real(_np.trace(_sqrtm(_np.dot(matr.conj().T, matr))))) 

def sparse_onenorm(A):
    """
    Computes the 1-norm of the scipy sparse matrix `A`.

    Parameters
    ----------
    A : scipy sparse matrix
        The matrix or vector to take the norm of.

    Returns
    -------
    float
    """
    return max(abs(A).sum(axis=0).flat)