Python scipy.sparse.linalg 模块，LinearOperator() 实例源码

def build_full(self, shape, fill_val):
        m, n = shape
        if fill_val == 0:
            return shape
        else:
            def matvec(v):
                return v.sum() * fill_val * np.ones(m)

            def rmatvec(v):
                return v.sum() * fill_val * np.ones(n)

            def matmat(M):
                return M.sum(axis=0) * fill_val * np.ones((m, M.shape[1]))

            return sla.LinearOperator(shape=shape,
                                      matvec=matvec,
                                      rmatvec=rmatvec,
                                      matmat=matmat,
                                      dtype=self.dtype)

def lagrangian_hessian(self, z, v):
        """Returns scaled Lagrangian Hessian"""
        # Compute Hessian in relation to x and s
        Hx = self.lagrangian_hessian_x(z, v)
        if self.n_ineq > 0:
            S_Hs_S = self.lagrangian_hessian_s(z, v)

        # The scaled Lagragian Hessian is:
        #     [ Hx    0    ]
        #     [ 0   S Hs S ]
        def matvec(vec):
            vec_x = self.get_variables(vec)
            vec_s = self.get_slack(vec)
            if self.n_ineq > 0:
                return np.hstack((Hx.dot(vec_x), S_Hs_S*vec_s))
            else:
                return Hx.dot(vec_x)
        return LinearOperator((self.n_vars+self.n_ineq,
                               self.n_vars+self.n_ineq),
                              matvec)

def J_q(x, c_lens):
    i = 0
    l = []
    for c_len in c_lens:
        x_ = x[i:i + c_len]

        x0 = x_[1:]
        x0norm = np.linalg.norm(x0)
        t0 = x_[0]

        if x0norm <= -t0:
            def matvec(v): return np.zeros(c_len)
        elif x0norm <= t0:
            def matvec(v): return v
        else:
            matvec = _soc(x0, x0norm, t0)
        l.append(sla.LinearOperator((c_len, c_len), matvec=matvec))
        i += c_len
    assert i == len(x), "{}, {}".format(i, len(x))
    return block_diag(l, arrtype=sla.LinearOperator)

def __call__(self, Y):
        """Return the Gram matrix associated with the data Y as a linear operator.

        .. math::
               K(X, Y)

        Parameters
        ----------
        Y : {array-like, sparse matrix}, shape = [n_samples1, n_features]
            Samples.

        Returns
        -------
        K(X, Y) : LinearOperator
            Returns K(X, Y).
        """
        return LinearOperator(
            (Y.shape[0] * self.p, self.n * self.p),
            dtype=self.X.dtype,
            matvec=lambda b: self._dot(self._Gram(Y), b),
            rmatvec=lambda b: self._dot(self._Gram(Y), b))

def __init__(self, X, A, scalar_kernel, scalar_kernel_params):
        """Initialize the Decomposable Operator-Valued Kernel.

        Parameters
        ----------

        X: {array-like, sparse matrix}, shape = [n_samples1, n_features]
            Support samples.

        A : {array, LinearOperator}, shape = [n_targets, n_targets]
            Linear operator acting on the outputs

        scalar_kernel : {callable}
            Callable which associate to the training points X the Gram matrix.

        scalar_kernel_params : {mapping of string to any}, optional
            Additional parameters (keyword arguments) for kernel function
            passed as callable object.
        """
        super(DecomposableKernelMap, self).__init__(A, scalar_kernel,
                                                    scalar_kernel_params)
        self.n = X.shape[0]
        self.d = X.shape[1]
        self.X = X
        self.Gs_train = None

def __call__(self, Y):
        """Return the Gram matrix associated with the data Y as a linear operator.

        .. math::
               K(X, Y)

        Parameters
        ----------
        Y : {array-like, sparse matrix}, shape = [n_samples1, n_features]
            Samples.

        Returns
        -------
        K(X, Y) : LinearOperator
            Returns K(X, Y).
        """
        return LinearOperator(
            (Y.shape[0] * self.p, self.n * self.p),
            dtype=self.X.dtype,
            matvec=lambda b: self._dot(self._Gram(Y), b),
            rmatvec=lambda b: self._dot(self._Gram(Y), b))

def __call__(self, Y):
        """Return the Gram matrix associated with the data Y as a linear operator.

        .. math::
               K(X, Y)

        Parameters
        ----------
        Y : {array-like, sparse matrix}, shape = [n_samples1, n_features]
            Samples.

        Returns
        -------
        K(X, Y) : LinearOperator
            Returns K(X, Y).
        """
        return LinearOperator(
            (Y.shape[0] * self.p, self.n * self.p),
            dtype=self.X.dtype,
            matvec=lambda b: self._dot(self._Gram(Y), b),
            rmatvec=lambda b: self._dot(self._Gram(Y), b))

def __init__(self, A, scalar_kernel=rbf_kernel, scalar_kernel_params=None):
        """Initialize the Decomposable Operator-Valued Kernel.

        Parameters
        ----------

        A : {array, LinearOperator}, shape = [n_targets, n_targets]
            Linear operator acting on the outputs

        scalar_kernel : {callable}
            Callable which associate to the training points X the Gram matrix.

        scalar_kernel_params : {mapping of string to any}, optional
            Additional parameters (keyword arguments) for kernel function
            passed as callable object.
        """
        self.A = A
        self.scalar_kernel = scalar_kernel
        self.scalar_kernel_params = scalar_kernel_params
        self.p = A.shape[0]

def get_subset_lin_op(lin_op, sub_idx):
    """ Subset a linear operator to the indices in `sub_idx`. Equivalent to A' = A[sub_idx, :]
    :param LinearOperator lin_op: input linear operator
    :param np.ndarray[int] sub_idx: subset index
    :return: the subset linear operator
    :rtype: LinearOperator
    """
    if lin_op is None:
        return None
    if type(lin_op) is IndexOperator:
        # subsetting IndexOperator yields a new IndexOperator
        return IndexOperator(lin_op.index_map[sub_idx], dim_x=lin_op.dim_x)
    elif isinstance(lin_op, MatrixLinearOperator):
        # subsetting a matrix multiplication operation yields a new matrix
        return MatrixLinearOperator(lin_op.A[sub_idx, :])
    # in the general case, append a sub-indexing operator
    return IndexOperator(sub_idx, dim_x=lin_op.shape[0]) * lin_op

def rmatvec_nd(lin_op, x):
    """
    Project a 1D or 2D numpy or sparse array using rmatvec. This is different from rmatvec
    because it applies rmatvec to each row and column. If x is n x n and lin_op is n x k,
    the result will be k x k.

    :param LinearOperator lin_op: The linear operator to apply to x
    :param np.ndarray|sp.spmatrix x: array/matrix to be projected
    :return: the projected array
    :rtype: np.ndarray|sp.spmatrix
    """
    if x is None or lin_op is None:
        return x
    if isinstance(x, sp.spmatrix):
        y = x.toarray()
    elif np.isscalar(x):
        y = np.array(x, ndmin=1)
    else:
        y = np.copy(x)
    proj_func = lambda z: lin_op.rmatvec(z)
    for j in range(y.ndim):
        if y.shape[j] == lin_op.shape[0]:
            y = np.apply_along_axis(proj_func, j, y)
    return y

def __init__(self,shape,matvec,dtype=None):
        '''
        Constructor.

        Parameters
        ----------
        shape : 2-tuple
            The shape of the linear operator.
        matvec : callable
            The matrix-vector multiplication function.
        dtype : np.float64, np.complex128, etc, optional
            The data type of the linear operator.
        '''
        super(LinearOperator,self).__init__(dtype=dtype,shape=shape)
        self.count=0
        self._matvec_=matvec

def _matvec_mull(self, v, sv):
        """
        helper function for linalg.LinearOperator class, acts as multiplication function for big matrix and vector
        without explicit forming a often big square matrix 
        """

        X = self._X
        l = self.l2reg

        y = l * v
        y[-1] = 0

        # Check which method is faster, with slicing support vectors before computation or after

        # 1. Method one
        # Xsv = X[sv]
        # compute dot products on support vectors only
        # z = Xsv.dot(v[0:-1]) + v[-1]
        # y = y + np.append(z.dot(Xsv), z.sum())

        #2. Method two
        # compute dot products on whole dataset
        z = X.dot(v[0:-1]) + v[-1]
        zz = np.zeros(z.shape[0])
        # choose support vectors values only
        zz[sv] = z[sv]
        y = y + np.append(zz.dot(X), zz.sum())

        return y

def ngstep(x0, obj0, objgrad0, obj_and_kl_func, hvpx0_func, max_kl, damping, max_cg_iter,
           enable_bt):
    '''
    Natural gradient step using hessian-vector products

    Args:
        x0: current point
        obj0: objective value at x0
        objgrad0: grad of objective value at x0
        obj_and_kl_func: function mapping a point x to the objective and kl values
        hvpx0_func: function mapping a vector v to the KL Hessian-vector product H(x0)v
        max_kl: max kl divergence limit. Triggers a line search.
        damping: multiple of I to mix with Hessians for Hessian-vector products
        max_cg_iter: max conjugate gradient iterations for solving for natural gradient step
    '''

    assert x0.ndim == 1 and x0.shape == objgrad0.shape

    # Solve for step direction
    damped_hvp_func = lambda v: hvpx0_func(v) + damping * v
    hvpop = ssl.LinearOperator(shape=(x0.shape[0], x0.shape[0]), matvec=damped_hvp_func)
    step, _ = ssl.cg(hvpop, -objgrad0, maxiter=max_cg_iter)
    fullstep = step / np.sqrt(.5 * step.dot(damped_hvp_func(step)) / max_kl + 1e-8)

    # Line search on objective with a hard KL wall
    if not enable_bt:
        return x0 + fullstep, 0

    def barrierobj(p):
        obj, kl = obj_and_kl_func(p)
        return np.inf if kl > 2 * max_kl else obj

    xnew, num_bt_steps = btlinesearch(f=barrierobj, x0=x0, fx0=obj0, g=objgrad0, dx=fullstep,
                                      accept_ratio=.1, shrink_factor=.5, max_steps=10)
    return xnew, num_bt_steps

def _get_backend(rows, dtype, arrtype):
    if arrtype == np.ndarray and dtype is not None:
        return NumpyBackend(arrtype, dtype)
    elif arrtype == sla.LinearOperator:
        return LinearOperatorBackend(dtype)
    elif arrtype is not None and re.search('torch\..*Tensor', repr(arrtype)):
        return TorchBackend(dtype)
    elif arrtype is not None and re.search('torch\..*(Variable|Parameter)', repr(arrtype)):
        return TorchVariableBackend(dtype)
    else:
        npb = NumpyBackend()
        tb = TorchBackend()
        lob = LinearOperatorBackend()
        tvb = TorchVariableBackend()
        for row in rows:
            for elem in row:
                if npb.is_complete(elem) and elem.size > 0:
                    if dtype is None:
                        dtype = type(elem[0, 0])
                    if arrtype is None:
                        arrtype = type(elem)
                    return NumpyBackend(dtype, arrtype)
                elif tb.is_complete(elem):
                    return TorchBackend(type(elem))
                elif lob.is_complete(elem):
                    return LinearOperatorBackend(elem.dtype)
                elif tvb.is_complete(elem):
                    return TorchVariableBackend(type(elem.data))

    assert(False)

def build(self, rows):
        col_sizes = [lo.shape[1] if self.is_complete(lo) else lo[1]
                     for lo in rows[0]]
        col_idxs = np.cumsum([0] + col_sizes)
        row_sizes = [row[0].shape[0] if self.is_complete(row[0]) else row[0][0]
                     for row in rows]
        row_idxs = np.cumsum([0] + row_sizes)
        m, n = sum(row_sizes), sum(col_sizes)

        def matvec(v):
            out = np.zeros(m)
            for row, i, j in zip(rows, row_idxs[:-1], row_idxs[1:]):
                out[i:j] = sum(lo.matvec(v[k:l]) for lo, k, l in
                               zip(row, col_idxs[:-1], col_idxs[1:])
                               if self.is_complete(lo))
            return out

        # The transposed list
        cols = zip(*rows)

        def rmatvec(v):
            out = np.zeros(n)
            for col, i, j in zip(cols, col_idxs[:-1], col_idxs[1:]):
                out[i:j] = sum(lo.rmatvec(v[k:l]) for lo, k, l in
                               zip(col, row_idxs[:-1], row_idxs[1:])
                               if self.is_complete(lo))
            return out

        def matmat(M):
            out = np.zeros((m, M.shape[1]))
            for row, i, j in zip(rows, row_idxs[:-1], row_idxs[1:]):
                out[i:j] = sum(lo.matmat(M[k:l]) for lo, k, l in
                               zip(row, col_idxs[:-1], col_idxs[1:])
                               if self.is_complete(lo))
            return out

        return sla.LinearOperator(shape=(m, n),
                                  matvec=matvec,
                                  rmatvec=rmatvec,
                                  matmat=matmat,
                                  dtype=self.dtype)

def is_complete(self, x):
        return isinstance(x, sla.LinearOperator)

def ngstep(x0, obj0, objgrad0, obj_and_kl_func, hvpx0_func, max_kl, damping, max_cg_iter, enable_bt):
    '''
    Natural gradient step using hessian-vector products

    Args:
        x0: current point
        obj0: objective value at x0
        objgrad0: grad of objective value at x0
        obj_and_kl_func: function mapping a point x to the objective and kl values
        hvpx0_func: function mapping a vector v to the KL Hessian-vector product H(x0)v
        max_kl: max kl divergence limit. Triggers a line search.
        damping: multiple of I to mix with Hessians for Hessian-vector products
        max_cg_iter: max conjugate gradient iterations for solving for natural gradient step
    '''

    assert x0.ndim == 1 and x0.shape == objgrad0.shape

    # Solve for step direction
    damped_hvp_func = lambda v: hvpx0_func(v) + damping*v
    hvpop = ssl.LinearOperator(shape=(x0.shape[0], x0.shape[0]), matvec=damped_hvp_func)
    step, _ = ssl.cg(hvpop, -objgrad0, maxiter=max_cg_iter)
    fullstep = step / np.sqrt(.5 * step.dot(damped_hvp_func(step)) / max_kl + 1e-8)

    # Line search on objective with a hard KL wall
    if not enable_bt:
        return x0+fullstep, 0

    def barrierobj(p):
        obj, kl = obj_and_kl_func(p)
        return np.inf if kl > 2*max_kl else obj
    xnew, num_bt_steps = btlinesearch(
        f=barrierobj,
        x0=x0,
        fx0=obj0,
        g=objgrad0,
        dx=fullstep,
        accept_ratio=.1, shrink_factor=.5, max_steps=10)
    return xnew, num_bt_steps

def _linear_operator_difference(fun, x0, f0, h, method):
    m = f0.size
    n = x0.size

    if method == '2-point':
        def matvec(p):
            if np.array_equal(p, np.zeros_like(p)):
                return np.zeros(m)
            dx = h / norm(p)
            x = x0 + dx*p
            df = fun(x) - f0
            return df / dx

    elif method == '3-point':
        def matvec(p):
            if np.array_equal(p, np.zeros_like(p)):
                return np.zeros(m)
            dx = 2*h / norm(p)
            x1 = x0 - (dx/2)*p
            x2 = x0 + (dx/2)*p
            f1 = fun(x1)
            f2 = fun(x2)
            df = f2 - f1
            return df / dx

    elif method == 'cs':
        def matvec(p):
            if np.array_equal(p, np.zeros_like(p)):
                return np.zeros(m)
            dx = h / norm(p)
            x = x0 + dx*p*1.j
            f1 = fun(x)
            df = f1.imag
            return df / dx

    else:
        raise RuntimeError("Never be here.")

    return LinearOperator((m, n), matvec)

def scaling(self, z):
        """Returns scaling vector.
        Given by:
            scaling = [ones(n_vars), s]
        """
        s = self.get_slack(z)
        diag_elements = np.hstack((np.ones(self.n_vars), s))

        # Diagonal Matrix
        def matvec(vec):
            return diag_elements*vec
        return LinearOperator((self.n_vars+self.n_ineq,
                               self.n_vars+self.n_ineq),
                              matvec)

def J_f(x):
    n = len(x)
    return sla.LinearOperator((n, n), matvec=lambda v: v)

def J_l(x):
    n = len(x)
    d = (x >= 0).astype(np.float64)
    return sla.LinearOperator((n, n), matvec=lambda v: d * v)

def J_ep(x, n_cones):
    Js = _exp_Js(x, n_cones)
    return block_diag(out, arrtype=sla.LinearOperator)

def J_ed(x, n_cones):
    Js = _exp_Js(-x, n_cones)
    I = np.eye(3)
    Js0 = [I - J for J in Js]
    return block_diag(Js0, arrtype=sla.LinearOperator)

def gen(self):
        shape_U = ((self.ns + self.ls) * self.p, (self.ns + self.ls) * self.p)
        shape_JT = ((self.ns + self.ls) * self.p, (self.ns + self.ls) * self.p)

        # return U, JT
        return (LinearOperator(shape_U,
                               dtype=self.L.dtype,
                               matvec=lambda b: self._dot_U(b),
                               rmatvec=lambda b: self._dot_U(b)),
                LinearOperator(shape_JT,
                               dtype=self.L.dtype,
                               matvec=lambda b: self._dot_JT(b),
                               rmatvec=lambda b: self._dot_JT(b)))

def __mul__(self, Ky):
        """Syntaxic sugar.

           If Kx is a compatible DotProduct kernel, returns

           .. math::
                K(X, Y) = K_X^T K_Y

        Parameters
        ----------
        Ky : {DotProductKernelMap}
            Compatible kernel Map (e.g. same kernel but different support data
            X).

        Returns
        -------
        K(X, Y) : LinearOperator
            Returns K(X, Y).

        Examples
        --------
        >>> import operalib as ovk
        >>> import numpy as np
        >>> X = np.random.randn(100, 10)
        >>> K = ovk.DotProductKernel(mu=0.2, p=5)
        >>> Gram = K(X, X)
        >>> Gram
        # doctest: +NORMALIZE_WHITESPACE +ELLIPSIS
        <500x500 _CustomLinearOperator with dtype=float64>
        >>> C = np.random.randn(Gram.shape[0])
        >>> Kx = K(X)  # The kernel map.
        >>> Ky = K(X)
        >>> np.allclose(Gram * C, (Kx.T * Ky) * C)
        # doctest: +NORMALIZE_WHITESPACE +ELLIPSIS
        True
        """
        # TODO: Check that Kx is compatible
        return self.__call__(Ky.X)

def __mul__(self, Ky):
        """Syntaxic sugar.

           If Kx is a compatible decomposable kernel, returns

           .. math::
                K(X, Y) = K_X^T K_Y

        Parameters
        ----------
        Ky : {DecomposableKernelMap}
            Compatible kernel Map (e.g. same kernel but different support data
            X).

        Returns
        -------
        K(X, Y) : LinearOperator
            Returns K(X, Y).

        Examples
        --------
        >>> import operalib as ovk
        >>> import numpy as np
        >>> X = np.random.randn(100, 10)
        >>> K = ovk.DecomposableKernel(np.eye(2))
        >>> Gram = K(X, X)
        >>> Gram
        # doctest: +NORMALIZE_WHITESPACE +ELLIPSIS
        <200x200 _CustomLinearOperator with dtype=float64>
        >>> C = np.random.randn(Gram.shape[0])
        >>> Kx = K(X)  # The kernel map.
        >>> Ky = K(X)
        >>> np.allclose(Gram * C, (Kx.T * Ky) * C)
        # doctest: +NORMALIZE_WHITESPACE +ELLIPSIS
        True
        """
        # TODO: Check that Kx is compatible
        return self.__call__(Ky.X)

def __call__(self, X, Y=None):
        r"""Return the kernel map associated with the data X.

        .. math::
               K_x: \begin{cases}
               Y \mapsto K(X, Y) \enskip\text{if } Y \text{is None,} \\
               K(X, Y) \enskip\text{otherwise.}
               \end{cases}

        Parameters
        ----------
        X : {array-like, sparse matrix}, shape = [n_samples1, n_features]
            Samples.

        Y : {array-like, sparse matrix}, shape = [n_samples2, n_features],
                                          default = None
            Samples.

        Returns
        -------
        K_x : DotProductKernelMap, callable or LinearOperator

            .. math::
               K_x: \begin{cases}
               Y \mapsto K(X, Y) \enskip\text{if } Y \text{is None,} \\
               K(X, Y) \enskip\text{otherwise}
               \end{cases}
        """
        Kmap = self.get_kernel_map(X)
        if Y is None:
            return Kmap
        else:
            return Kmap(Y)

def __call__(self, X, Y=None):
        r"""Return the kernel map associated with the data X.

        .. math::
               K_x: \begin{cases}
               Y \mapsto K(X, Y) \enskip\text{if } Y \text{is None,} \\
               K(X, Y) \enskip\text{otherwise.}
               \end{cases}

        Parameters
        ----------
        X : {array-like, sparse matrix}, shape = [n_samples1, n_features]
            Samples.

        Y : {array-like, sparse matrix}, shape = [n_samples2, n_features],
                                          default = None
            Samples.

        Returns
        -------
        K_x : DecomposableKernelMap, callable or LinearOperator

            .. math::
               K_x: \begin{cases}
               Y \mapsto K(X, Y) \enskip\text{if } Y \text{is None,} \\
               K(X, Y) \enskip\text{otherwise}
               \end{cases}
        """
        Kmap = self.get_kernel_map(X)
        if Y is None:
            return Kmap
        else:
            return Kmap(Y)

def __call__(self, X, Y=None):
        r"""Return the kernel map associated with the data X.

        .. math::
               K_x: \begin{cases}
               Y \mapsto K(X, Y) \enskip\text{if } Y \text{is None,} \\
               K(X, Y) \enskip\text{otherwise.}
               \end{cases}

        Parameters
        ----------
        X : {array-like, sparse matrix}, shape = [n_samples1, n_features]
            Samples.

        Y : {array-like, sparse matrix}, shape = [n_samples2, n_features],
                                          default = None
            Samples.

        Returns
        -------
        K_x : DecomposableKernelMap, callable or LinearOperator

        .. math::
            K_x: \begin{cases}
            Y \mapsto K(X, Y) \enskip\text{if } Y \text{is None,} \\
            K(X, Y) \enskip\text{otherwise}
            \end{cases}
        """
        Kmap = self.get_kernel_map(X)
        if Y is None:
            return Kmap
        else:
            return Kmap(Y)

def __call__(self, X, Y=None):
        r"""Return the kernel map associated with the data X.

        .. math::
               K_x: \begin{cases}
               Y \mapsto K(X, Y) \enskip\text{if } Y \text{is None,} \\
               K(X, Y) \enskip\text{otherwise.}
               \end{cases}

        Parameters
        ----------
        X : {array-like, sparse matrix}, shape = [n_samples1, n_features]
            Samples.

        Y : {array-like, sparse matrix}, shape = [n_samples2, n_features],
                                          default = None
            Samples.

        Returns
        -------
        K_x : DecomposableKernelMap, callable or LinearOperator

        .. math::
            K_x: \begin{cases}
            Y \mapsto K(X, Y) \enskip\text{if } Y \text{is None,} \\
            K(X, Y) \enskip\text{otherwise}
            \end{cases}
        """
        Kmap = self.get_kernel_map(X)
        if Y is None:
            return Kmap
        else:
            return Kmap(Y)

def eigsh(A,max_try=6,return_eigenvectors=True,**karg):
    '''
    Find the eigenvalues and eigenvectors of the real symmetric square matrix or complex hermitian matrix A.
    This is a wrapper for scipy.sparse.linalg.eigsh to handle the exceptions it raises.

    Parameters
    ----------
    A : An NxN matrix, array, sparse matrix, or LinearOperator
        The matrix whose eigenvalues and eigenvectors is to be computed.
    max_try : integer, optional
        The maximum number of tries to do the computation when the computed eigenvalues/eigenvectors do not converge.
    return_eigenvectors : logical, optional
        True for returning the eigenvectors and False for not.
    karg : dict
        Please refer to https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.linalg.eigsh.html for details.
    '''
    if A.shape==(1,1):
        assert 'M' not in karg and 'sigma' not in karg and karg.get('k',1)==1
        if return_eigenvectors:
            return A.dot(np.ones(1)).reshape(-1),np.ones((1,1),dtype=A.dtype)
        else:
            return A.dot(np.ones(1)).reshape(-1)
    else:
        ntry=1
        while True:
            try:
                result=pl.eigsh(A,return_eigenvectors=return_eigenvectors,**karg)
                break
            except pl.ArpackNoConvergence as err:
                if ntry<max_try:
                    ntry+=1
                else:
                    raise err
        return result

def _computePreconditioner(self, system_matrix):
        '''
        Create a preconditioner base on
        an incomplete LU decomposition
        '''
        # convert to compressed column storage format (efficiency)
        sm_csc = system_matrix.tocsc()
        lu = sparse_solvers.spilu(sm_csc)
        m_operator = lambda x: lu.solve(x)
        shape = (self.nx*self.ny, self.nx*self.ny)
        self.preconditioner = sparse_solvers.LinearOperator(shape, m_operator)

def solvemdbi_cg(ah, rho, b, axisM, axisK, tol=1e-5, mit=1000, isn=None):
    r"""
    Solve a multiple diagonal block linear system with a scaled
    identity term using Conjugate Gradient (CG) via
    :func:`scipy.sparse.linalg.cg`.

    The solution is obtained by independently solving a set of linear
    systems of the form (see :cite:`wohlberg-2016-efficient`)

     .. math::
      (\rho I + \mathbf{a}_0 \mathbf{a}_0^H + \mathbf{a}_1 \mathbf{a}_1^H +
       \; \ldots \; + \mathbf{a}_{K-1} \mathbf{a}_{K-1}^H) \; \mathbf{x} =
       \mathbf{b}

    where each :math:`\mathbf{a}_k` is an :math:`M`-vector.
    The inner products and matrix products in this equation are taken
    along the M and K axes of the corresponding multi-dimensional arrays;
    the solutions are independent over the other axes.

    Parameters
    ----------
    ah : array_like
      Linear system component :math:`\mathbf{a}^H`
    rho : float
      Parameter rho
    b : array_like
      Linear system component :math:`\mathbf{b}`
    axisM : int
      Axis in input corresponding to index m in linear system
    axisK : int
      Axis in input corresponding to index k in linear system
    tol : float
      CG tolerance
    mit : int
      CG maximum iterations
    isn : array_like
      CG initial solution

    Returns
    -------
    x : ndarray
      Linear system solution :math:`\mathbf{x}`
    cgit : int
      Number of CG iterations
    """

    a = np.conj(ah)
    if isn is not None:
        isn = isn.ravel()
    Aop = lambda x: inner(ah, x, axis=axisM)
    AHop = lambda x: inner(a, x, axis=axisK)
    AHAop = lambda x: AHop(Aop(x))
    vAHAoprI = lambda x: AHAop(x.reshape(b.shape)).ravel() + rho*x.ravel()
    lop = LinearOperator((b.size, b.size), matvec=vAHAoprI, dtype=b.dtype)
    vx, cgit = cg(lop, b.ravel(), isn, tol, mit)
    return vx.reshape(b.shape), cgit

def _solve_CG(self, X, Y):
        """
        Solve the primal SVM problem with Newton method without computing hessina matrix explicit,
        good for big sparse matrix
        :param X: 
        :param Y: 
        """
        [n, d] = X.shape

        # we add one last component, which is b (bias)
        self.w = np.zeros(d + 1)

        # helper variable for storing 1-Y*(np.dot(X,w))
        self.out = np.ones(n)

        l = self.l2reg
        # the number of alg. iteration
        iter = 0

        sv = np.where(self.out > 0)[0]

        # create linear operator, acts as matrix vector multiplication, without storing full matrix(hessian)
        #hess_vec = linalg.LinearOperator((d + 1, d + 1), matvec=self._matvec_mull)

        # This is a hack in order to pass additional parameters to linear matvec function
        mv2 = lambda v: self._matvec_mull(v, sv)
        # create linear operator, acts as matrix vector multiplication, without storing full matrix(hessian)
        hess_vec = linalg.LinearOperator((d + 1, d + 1), matvec=mv2)

        while True:
            iter = iter + 1
            if iter > self.newton_iter:
                print("Maximum {0} of Newton steps reached, change newton_iter parameter or try larger lambda".format(
                    iter))
                break

            obj, grad = self._obj_func(self.w, X, Y, self.out)

            # np.where returns a tuple, we take the first dim
            sv = np.where(self.out > 0)[0]

            step, info = linalg.minres(hess_vec, -grad)

            t, self.out = self._line_search(self.w, step, self.out)

            self.w += t * step

            if -step.dot(grad) < self._prec * obj:
                break

def gram_block_circulant(ds, n_times_valid, method='full',
                         sample_weights=None):
    """Returns ...

    Parameters
    ----------
    ds : array, shape (n_atoms, n_times_atom)
        The atoms
    n_times_valid : int
        n_times - n_times_atom + 1
    method : string
        If 'full', returns full circulant matrix.
        If 'scipy', returns scipy linear operator.
        If 'custom', returns custom linear operator.
    sample_weights : array, shape (n_times, )
        The sample weights for one trial
    """
    from scipy.sparse.linalg import LinearOperator
    from functools import partial

    n_atoms, n_times_atom = ds.shape
    n_times = n_times_valid + n_times_atom - 1

    if method == 'full':
        D = np.zeros((n_times, n_atoms * n_times_valid))
        for k_idx in range(n_atoms):
            d_padded = np.zeros((n_times, ))
            d_padded[:n_times_atom] = ds[k_idx]
            start = k_idx * n_times_valid
            stop = start + n_times_valid
            D[:, start:stop] = linalg.circulant((d_padded))[:, :n_times_valid]
        if sample_weights is not None:
            wD = sample_weights[:, None] * D
            return np.dot(D.T, wD)
        else:
            return np.dot(D.T, D)

    elif method == 'scipy':

        def matvec(v, ds):
            assert v.shape[0] % ds.shape[0] == 0
            return _fprime(ds, v, Xi=None, reg=None,
                           sample_weights=sample_weights)

        D = LinearOperator((n_atoms * n_times_valid, n_atoms * n_times_valid),
                           matvec=partial(matvec, ds=ds))

    elif method == 'custom':
        return CustomLinearOperator(ds, n_times_valid, sample_weights)
    else:
        raise ValueError('Unkown method %s.' % method)
    return D

def build_eye(self, n):
        def identity(v): return v
        return sla.LinearOperator(shape=(n, n),
                                  matvec=identity,
                                  rmatvec=identity,
                                  matmat=identity,
                                  dtype=self.dtype)

def get_orff_map(self, X, D=100, random_state=0):
        r"""Return the Random Fourier Feature map associated with the data X.

        .. math::
               K_x: Y \mapsto \tilde{\Phi}(X)

        Parameters
        ----------
        X : {array-like, sparse matrix}, shape = [n_samples, n_features]
            Samples.

        Returns
        -------
        \tilde{\Phi}(X) : Linear Operator, callable
        """
        self.r = 1
        if not hasattr(self, 'Xb_'):
            self.phi_ = RBFSampler(gamma=self.gamma,
                                   n_components=D, random_state=random_state)
            self.phi_.fit(X)
            self.Xb_ = self.phi_.transform(X)
            self.Xb_ = (self.Xb_.reshape((self.Xb_.shape[0],
                                          1, self.Xb_.shape[1])) *
                        self.phi_.random_weights_.reshape((1, -1,
                                                           self.Xb_.shape[1])))
            self.Xb_ = self.Xb_.reshape((-1, self.Xb_.shape[2]))

        D = self.phi_.n_components
        if X is self.Xb_:
            return LinearOperator(self.Xb_.shape,
                                  matvec=lambda b: dot(self.Xb_ * b),
                                  rmatvec=lambda r: dot(self.Xb_.T * r))
        else:
            Xb = self.phi_.transform(X)
            # TODO:
            # w = self.phi_.random_weights_.reshape((1, -1, Xb.shape[1]))
            # wn = np.linalg.norm(w)
            # Xb = (Xb.reshape((Xb.shape[0], 1, Xb.shape[1])) *
            #       wn * np.eye()w np.dot(w.T, w) / wn)
            Xb = Xb.reshape((-1, Xb.shape[2]))
            return LinearOperator(Xb.shape,
                                  matvec=lambda b: dot(Xb, b),
                                  rmatvec=lambda r: dot(Xb.T, r))

def factor_aug(z, DPhival, G, A):
    r"""Set up augmented system and return.

    Parameters
    ----------
    z : (N+P+M+M,) ndarray
        Current iterate, z = (x, nu, l, s)
    DPhival : LinearOperator
        Jacobian of the variational inequality mapping
    G : (M, N) ndarray or sparse matrix
        Inequality constraints
    A : (P, N) ndarray or sparse matrix
        Equality constraints  

    Returns
    -------
    J : LinearOperator
        Augmented system

    """
    M, N = G.shape
    P, N = A.shape
    """Multiplier for inequality constraints"""
    l = z[N+P:N+P+M]

    """Slacks"""
    s = z[N+P+M:]

    """Sigma matrix"""
    SIG = diags(l/s, 0)
    # SIG = diags(l*s, 0)

    """Convert A"""
    if not issparse(A):
        A = csr_matrix(A)

    """Convert G"""
    if not issparse(G):
        G = csr_matrix(G)

    """Since we expect symmetric DPhival, we need to change A"""
    sign = np.zeros(N)
    sign[0:N/2] = 1.0
    sign[N/2:] = -1.0
    T = diags(sign, 0)

    A_new = A.dot(T)

    W = AugmentedSystem(DPhival, G, SIG, A_new)
    return W

def _solve_sparse_cg(X, y, alpha, max_iter=None, tol=1e-3, verbose=0):
    n_samples, n_features = X.shape
    X1 = sp_linalg.aslinearoperator(X)
    coefs = np.empty((y.shape[1], n_features))

    if n_features > n_samples:
        def create_mv(curr_alpha):
            def _mv(x):
                return X1.matvec(X1.rmatvec(x)) + curr_alpha * x
            return _mv
    else:
        def create_mv(curr_alpha):
            def _mv(x):
                return X1.rmatvec(X1.matvec(x)) + curr_alpha * x
            return _mv

    for i in range(y.shape[1]):
        y_column = y[:, i]

        mv = create_mv(alpha[i])
        if n_features > n_samples:
            # kernel ridge
            # w = X.T * inv(X X^t + alpha*Id) y
            C = sp_linalg.LinearOperator(
                (n_samples, n_samples), matvec=mv, dtype=X.dtype)
            coef, info = sp_linalg.cg(C, y_column, tol=tol)
            coefs[i] = X1.rmatvec(coef)
        else:
            # linear ridge
            # w = inv(X^t X + alpha*Id) * X.T y
            y_column = X1.rmatvec(y_column)
            C = sp_linalg.LinearOperator(
                (n_features, n_features), matvec=mv, dtype=X.dtype)
            coefs[i], info = sp_linalg.cg(C, y_column, maxiter=max_iter,
                                          tol=tol)
        if info < 0:
            raise ValueError("Failed with error code %d" % info)

        if max_iter is None and info > 0 and verbose:
            warnings.warn("sparse_cg did not converge after %d iterations." %
                          info)

    return coefs