Python scipy.sparse.linalg.lsqr() Examples

The following are 30 code examples of scipy.sparse.linalg.lsqr(). You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may also want to check out all available functions/classes of the module scipy.sparse.linalg , or try the search function .
Example #1
Source File: test_dwts.py    From pylops with GNU Lesser General Public License v3.0 6 votes vote down vote up 
def test_DWT_1dsignal(par):
    """Dot-test and inversion for DWT operator for 1d signal
    """
    DWTop = DWT(dims=[par['nt']], dir=0, wavelet='haar', level=3)
    x = np.random.normal(0., 1., par['nt']) + \
        par['imag'] * np.random.normal(0., 1., par['nt'])

    assert dottest(DWTop, DWTop.shape[0], DWTop.shape[1],
                   complexflag=0 if par['imag'] == 0 else 3)

    y = DWTop * x
    xadj = DWTop.H * y  # adjoint is same as inverse for dwt
    xinv = lsqr(DWTop, y, damp=1e-10, iter_lim=10, show=0)[0]

    assert_array_almost_equal(x, xadj, decimal=8)
    assert_array_almost_equal(x, xinv, decimal=8)
Example #2
Source File: ridge.py    From twitter-stock-recommendation with MIT License 6 votes vote down vote up 
def _solve_lsqr(X, y, alpha, max_iter=None, tol=1e-3):
    n_samples, n_features = X.shape
    coefs = np.empty((y.shape[1], n_features), dtype=X.dtype)
    n_iter = np.empty(y.shape[1], dtype=np.int32)

    # According to the lsqr documentation, alpha = damp^2.
    sqrt_alpha = np.sqrt(alpha)

    for i in range(y.shape[1]):
        y_column = y[:, i]
        info = sp_linalg.lsqr(X, y_column, damp=sqrt_alpha[i],
                              atol=tol, btol=tol, iter_lim=max_iter)
        coefs[i] = info[0]
        n_iter[i] = info[2]

    return coefs, n_iter
Example #3
Source File: optimization_internal.py    From chumpy with MIT License 6 votes vote down vote up 
def updateGN(self):
        tm = timer()
        if sp.issparse(self.A):
            self.A.eliminate_zeros()
            pif('sparse solve...sparsity infill is %.3f%% (hessian %dx%d)' % (100. * self.A.nnz / (self.A.shape[0] * self.A.shape[1]), self.A.shape[0], self.A.shape[1]))
            if self.g.size > 1:
                self.d_gn = self.solve(self.A, self.g).ravel()
                if np.any(np.isnan(self.d_gn)) or np.any(np.isinf(self.d_gn)):
                    from scipy.sparse.linalg import lsqr
                    warnings.warn("sparse solve failed, falling back to lsqr")
                    self.d_gn = lsqr(self.A, self.g)[0].ravel()
            else:
                self.d_gn = np.atleast_1d(self.g.ravel()[0]/self.A[0,0])
            pif('sparse solve...done in %.2fs' % tm())
        else:
            pif('dense solve...')
            try:
                self.d_gn = np.linalg.solve(self.A, self.g).ravel()
            except Exception:
                warnings.warn("dense solve failed, falling back to lsqr")
                self.d_gn = np.linalg.lstsq(self.A, self.g)[0].ravel()
            pif('dense solve...done in %.2fs' % tm())
Example #4
Source File: ridge.py    From Splunking-Crime with GNU Affero General Public License v3.0 6 votes vote down vote up 
def _solve_lsqr(X, y, alpha, max_iter=None, tol=1e-3):
    n_samples, n_features = X.shape
    coefs = np.empty((y.shape[1], n_features), dtype=X.dtype)
    n_iter = np.empty(y.shape[1], dtype=np.int32)

    # According to the lsqr documentation, alpha = damp^2.
    sqrt_alpha = np.sqrt(alpha)

    for i in range(y.shape[1]):
        y_column = y[:, i]
        info = sp_linalg.lsqr(X, y_column, damp=sqrt_alpha[i],
                              atol=tol, btol=tol, iter_lim=max_iter)
        coefs[i] = info[0]
        n_iter[i] = info[2]

    return coefs, n_iter
Example #5
Source File: randomwalk.py    From pykernels with MIT License 6 votes vote down vote up 
def _compute(self, data_1, data_2):
        data_1 = basic.graphs_to_adjacency_lists(data_1)
        data_2 = basic.graphs_to_adjacency_lists(data_2)
        res = np.zeros((len(data_1), len(data_2)))
        N = len(data_1) * len(data_2)
        for i, graph1 in enumerate(data_1):
            for j, graph2 in enumerate(data_2):
                # norm1, norm2 - normalized adjacency matrixes
                norm1 = _norm(graph1)
                norm2 = _norm(graph2)
                # if graph is unweighted, W_prod = kron(a_norm(g1)*a_norm(g2))
                w_prod = kron(lil_matrix(norm1), lil_matrix(norm2))
                starting_prob = np.ones(w_prod.shape[0]) / (w_prod.shape[0])
                stop_prob = starting_prob
                # first solve (I - lambda * W_prod) * x = p_prod
                A = identity(w_prod.shape[0]) - (w_prod * self._lmb)
                x = lsqr(A, starting_prob)
                res[i, j] = stop_prob.T.dot(x[0])
                # print float(len(data_2)*i + j)/float(N), "%"
        return res
Example #6
Source File: inference.py    From ektelo with Apache License 2.0 6 votes vote down vote up 
def infer(self, Ms, ys, scale_factors=None):
        ''' Either:
            1) Ms is a single M and ys is a single y 
               (scale_factors ignored) or
            2) Ms and ys are lists of M matrices and y vectors
               and scale_factors is a list of the same length.
        '''
        A, y = _apply_scales(Ms, ys, scale_factors)

        if self.method == 'standard':
            assert self.l2_reg == 0, 'l2 reg not supported with method=standard'
            (x_est, _, rank, _) = linalg.lstsq(A.dense_matrix(), y, lapack_driver='gelsy')
        elif self.method == 'lsmr':
            res = lsmr(A, y, atol=0, btol=0, damp=self.l2_reg)
            x_est = res[0]
        elif self.method == 'lsqr':
            res = lsqr(A, y, atol=0, btol=0, damp=self.l2_reg)
            x_est = res[0]

        return x_est
Example #7
Source File: ridge.py    From Mastering-Elasticsearch-7.0 with MIT License 6 votes vote down vote up 
def _solve_lsqr(X, y, alpha, max_iter=None, tol=1e-3):
    n_samples, n_features = X.shape
    coefs = np.empty((y.shape[1], n_features), dtype=X.dtype)
    n_iter = np.empty(y.shape[1], dtype=np.int32)

    # According to the lsqr documentation, alpha = damp^2.
    sqrt_alpha = np.sqrt(alpha)

    for i in range(y.shape[1]):
        y_column = y[:, i]
        info = sp_linalg.lsqr(X, y_column, damp=sqrt_alpha[i],
                              atol=tol, btol=tol, iter_lim=max_iter)
        coefs[i] = info[0]
        n_iter[i] = info[2]

    return coefs, n_iter
Example #8
Source File: OLSUngnarStyle.py    From scattertext with Apache License 2.0 6 votes vote down vote up 
def get_scores_and_p_values(self, tdm, category):
		'''
		Parameters
		----------
		tdm: TermDocMatrix
		category: str, category name

		Returns
		-------
		pd.DataFrame(['coef', 'p-val'])
		'''
		X = tdm._X
		y = self._make_response_variable_1_or_negative_1(category, tdm)
		pX = X / X.sum(axis=1)
		ansX = self._anscombe_transform(pX.copy())
		B, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var\
			= lsqr(A=ansX, b=y, calc_var=True)
Example #9
Source File: test_dwts.py    From pylops with GNU Lesser General Public License v3.0 6 votes vote down vote up 
def test_DWT2D_3dsignal(par):
    """Dot-test and inversion for DWT operator for 3d signal
    """
    for dirs in [(0, 1), (0, 2), (1, 2)]:
        DWTop = DWT2D(dims=(par['nt'], par['nx'], par['ny']),
                      dirs=dirs, wavelet='haar', level=3)
        x = np.random.normal(0., 1., (par['nt'], par['nx'], par['ny'])) + \
            par['imag'] * np.random.normal(0., 1., (par['nt'], par['nx'], par['ny']))

        assert dottest(DWTop, DWTop.shape[0], DWTop.shape[1],
                       complexflag=0 if par['imag'] == 0 else 3)

        y = DWTop * x.ravel()
        xadj = DWTop.H * y  # adjoint is same as inverse for dwt
        xinv = lsqr(DWTop, y, damp=1e-10, iter_lim=10, show=0)[0]

        assert_array_almost_equal(x.ravel(), xadj, decimal=8)
        assert_array_almost_equal(x.ravel(), xinv, decimal=8)
Example #10
Source File: test_dwts.py    From pylops with GNU Lesser General Public License v3.0 6 votes vote down vote up 
def test_DWT2D_2dsignal(par):
    """Dot-test and inversion for DWT2D operator for 2d signal
    """
    DWTop = DWT2D(dims=(par['nt'], par['nx']),
                dirs=(0, 1), wavelet='haar', level=3)
    x = np.random.normal(0., 1., (par['nt'], par['nx'])) + \
        par['imag'] * np.random.normal(0., 1., (par['nt'], par['nx']))

    assert dottest(DWTop, DWTop.shape[0], DWTop.shape[1],
                   complexflag=0 if par['imag'] == 0 else 3)

    y = DWTop * x.ravel()
    xadj = DWTop.H * y  # adjoint is same as inverse for dwt
    xinv = lsqr(DWTop, y, damp=1e-10, iter_lim=10, show=0)[0]

    assert_array_almost_equal(x.ravel(), xadj, decimal=8)
    assert_array_almost_equal(x.ravel(), xinv, decimal=8)
Example #11
Source File: test_dwts.py    From pylops with GNU Lesser General Public License v3.0 6 votes vote down vote up 
def test_DWT_3dsignal(par):
    """Dot-test and inversion for DWT operator for 3d signal
    """
    for dir in [0, 1, 2]:
        DWTop = DWT(dims=(par['nt'], par['nx'], par['ny']),
                    dir=dir, wavelet='haar', level=3)
        x = np.random.normal(0., 1., (par['nt'], par['nx'], par['ny'])) + \
            par['imag'] * np.random.normal(0., 1., (par['nt'], par['nx'], par['ny']))

        assert dottest(DWTop, DWTop.shape[0], DWTop.shape[1],
                       complexflag=0 if par['imag'] == 0 else 3)

        y = DWTop * x.ravel()
        xadj = DWTop.H * y  # adjoint is same as inverse for dwt
        xinv = lsqr(DWTop, y, damp=1e-10, iter_lim=10, show=0)[0]

        assert_array_almost_equal(x.ravel(), xadj, decimal=8)
        assert_array_almost_equal(x.ravel(), xinv, decimal=8)
Example #12
Source File: test_dwts.py    From pylops with GNU Lesser General Public License v3.0 6 votes vote down vote up 
def test_DWT_2dsignal(par):
    """Dot-test and inversion for DWT operator for 2d signal
    """
    for dir in [0, 1]:
        DWTop = DWT(dims=(par['nt'], par['nx']),
                    dir=dir, wavelet='haar', level=3)
        x = np.random.normal(0., 1., (par['nt'], par['nx'])) + \
            par['imag'] * np.random.normal(0., 1., (par['nt'], par['nx']))

        assert dottest(DWTop, DWTop.shape[0], DWTop.shape[1],
                       complexflag=0 if par['imag'] == 0 else 3)

        y = DWTop * x.ravel()
        xadj = DWTop.H * y  # adjoint is same as inverse for dwt
        xinv = lsqr(DWTop, y, damp=1e-10, iter_lim=10, show=0)[0]

        assert_array_almost_equal(x.ravel(), xadj, decimal=8)
        assert_array_almost_equal(x.ravel(), xinv, decimal=8)
Example #13
Source File: LaplacianMesh.py    From laplacian-meshes with GNU General Public License v3.0 6 votes vote down vote up 
def solveLaplacianMesh(mesh, anchors, anchorsIdx, cotangent=True):
    n = mesh.VPos.shape[0] # N x 3
    k = anchorsIdx.shape[0]

    operator = (getLaplacianMatrixUmbrella, getLaplacianMatrixCotangent)

    L = operator[1](mesh, anchorsIdx) if cotangent else operator[0](mesh, anchorsIdx)
    delta = np.array(L.dot(mesh.VPos))
    
    # augment delta solution matrix with weighted anchors
    for i in range(k):
        delta[n + i, :] = WEIGHT * anchors[i, :]

    # update mesh vertices with least-squares solution
    for i in range(3):
        mesh.VPos[:, i] = lsqr(L, delta[:, i])[0]
    
    return mesh
#Purpose: Given a few RGB colors on a mesh, smoothly interpolate those colors
#by using their values as anchors and 
#Inputs: mesh (polygon mesh object), anchors (a K x 3 numpy array of anchor
#coordinates), colorsIdx (a parallel array of the indices of the RGB anchor indices)
Example #14
Source File: test_diagonal.py    From pylops with GNU Lesser General Public License v3.0 6 votes vote down vote up 
def test_Diagonal_3dsignal(par):
    """Dot-test and inversion for Diagonal operator for 3d signal
    """
    for idim, ddim in enumerate((par['ny'], par['nx'], par['nt'])):
        d = np.arange(ddim) + 1. +\
            par['imag'] * (np.arange(ddim) + 1.)

        Dop = Diagonal(d, dims=(par['ny'], par['nx'], par['nt']),
                       dir=idim, dtype=par['dtype'])
        assert dottest(Dop, par['ny']*par['nx']*par['nt'],
                       par['ny']*par['nx']*par['nt'],
                       complexflag=0 if par['imag'] == 0 else 3)

        x = np.ones((par['ny'], par['nx'], par['nt'])) + \
            par['imag']*np.ones((par['ny'], par['nx'], par['nt']))
        xlsqr = lsqr(Dop, Dop * x.ravel(), damp=1e-20, iter_lim=300, show=0)[0]

        assert_array_almost_equal(x.ravel(), xlsqr.ravel(), decimal=4)
Example #15
Source File: test_diagonal.py    From pylops with GNU Lesser General Public License v3.0 6 votes vote down vote up 
def test_Diagonal_2dsignal(par):
    """Dot-test and inversion for Diagonal operator for 2d signal
    """
    for idim, ddim in enumerate((par['nx'], par['nt'])):
        d = np.arange(ddim) + 1. +\
            par['imag'] * (np.arange(ddim) + 1.)

        Dop = Diagonal(d, dims=(par['nx'], par['nt']),
                       dir=idim, dtype=par['dtype'])
        assert dottest(Dop, par['nx']*par['nt'], par['nx']*par['nt'],
                       complexflag=0 if par['imag'] == 0 else 3)

        x = np.ones((par['nx'], par['nt'])) + \
            par['imag']*np.ones((par['nx'], par['nt']))
        xlsqr = lsqr(Dop, Dop * x.ravel(), damp=1e-20, iter_lim=300, show=0)[0]

        assert_array_almost_equal(x.ravel(), xlsqr.ravel(), decimal=4)
Example #16
Source File: test_fredholm.py    From pylops with GNU Lesser General Public License v3.0 6 votes vote down vote up 
def test_Fredholm1(par):
    """Dot-test and inversion for Fredholm1 operator
    """
    np.random.seed(10)

    _F = np.arange(par['nsl'] * par['nx'] * par['ny']).reshape(par['nsl'],
                                                               par['nx'],
                                                               par['ny'])
    F = _F - par['imag'] * _F

    x = np.ones((par['nsl'], par['ny'], par['nz'])) + \
        par['imag'] * np.ones((par['nsl'], par['ny'], par['nz']))

    Fop = Fredholm1(F, nz=par['nz'], saveGt=par['saveGt'],
                    usematmul=par['usematmul'], dtype=par['dtype'])
    assert dottest(Fop, par['nsl']*par['nx']*par['nz'],
                   par['nsl']*par['ny']*par['nz'],
                   complexflag=0 if par['imag'] == 0 else 3)
    xlsqr = lsqr(Fop, Fop * x.flatten(), damp=1e-20,
                 iter_lim=30, show=0)[0]
    xlsqr = xlsqr.reshape(par['nsl'], par['ny'], par['nz'])
    assert_array_almost_equal(x, xlsqr, decimal=3)
Example #17
Source File: LaplacianMesh.py    From laplacian-meshes with GNU General Public License v3.0 6 votes vote down vote up 
def smoothColors(mesh, anchors, colorsIdx):
    colorsIdx = np.array(colorsIdx)
    n = mesh.VPos.shape[0] 
    k = anchors.shape[0]
    colors = np.zeros((n, 3))
    delta = np.zeros((n + k, 3))

    L = getLaplacianMatrixUmbrella(mesh, colorsIdx);

     # augment delta solution matrix with weighted anchors
    for i in range(k):
        delta[n + i, :] = WEIGHT * anchors[i, :]
    
    # update RGB values with least-squares solution
    for i in range(3):
        colors[:, i] = lsqr(L, delta[:, i])[0]

    return colors

#Purpose: Given a mesh, to smooth it by subtracting off the delta coordinates
#from each vertex, normalized by the degree of that vertex
#Inputs: mesh (polygon mesh object)
#Returns: Nothing (should update mesh.VPos)
Example #18
Source File: test_basicoperators.py    From pylops with GNU Lesser General Public License v3.0 6 votes vote down vote up 
def test_Regression(par):
    """Dot-test, inversion and apply for Regression operator
    """
    np.random.seed(10)
    order = 4
    t = np.arange(par['ny'], dtype=np.float32)
    LRop = Regression(t, order=order, dtype=par['dtype'])
    assert dottest(LRop, par['ny'], order+1)

    x = np.array([1., 2., 0., 3., -1.], dtype=np.float32)
    xlsqr = lsqr(LRop, LRop*x, damp=1e-10, iter_lim=300, show=0)[0]
    assert_array_almost_equal(x, xlsqr, decimal=3)

    y = LRop * x
    y1 = LRop.apply(t, x)
    assert_array_almost_equal(y, y1, decimal=3)
Example #19
Source File: test_kronecker.py    From pylops with GNU Lesser General Public License v3.0 6 votes vote down vote up 
def test_Kroneker(par):
    """Dot-test, inversion and comparison with np.kron for Kronecker operator
    """
    np.random.seed(10)
    G1 = np.random.normal(0, 10, (par['ny'], par['nx'])).astype(par['dtype'])
    G2 = np.random.normal(0, 10, (par['ny'], par['nx'])).astype(par['dtype'])
    x = np.ones(par['nx']**2) + par['imag']*np.ones(par['nx']**2)

    Kop = Kronecker(MatrixMult(G1, dtype=par['dtype']),
                    MatrixMult(G2, dtype=par['dtype']),
                    dtype=par['dtype'])
    assert dottest(Kop, par['ny']**2, par['nx']**2,
                   complexflag=0 if par['imag'] == 0 else 3)

    xlsqr = lsqr(Kop, Kop * x, damp=1e-20, iter_lim=300, show=0)[0]
    assert_array_almost_equal(x, xlsqr, decimal=2)

    # Comparison with numpy
    assert_array_almost_equal(np.kron(G1, G2), Kop * np.eye(par['nx']**2),
                              decimal=3)
Example #20
Source File: test_basicoperators.py    From pylops with GNU Lesser General Public License v3.0 6 votes vote down vote up 
def test_MatrixMult_repeated(par):
    """Dot-test and inversion for test_MatrixMult operator repeated
    along another dimension
    """
    np.random.seed(10)
    G = np.random.normal(0, 10, (par['ny'], par['nx'])).astype('float32') + \
        par['imag'] * np.random.normal(0, 10, (par['ny'],
                                               par['nx'])).astype('float32')
    Gop = MatrixMult(G, dims=5, dtype=par['dtype'])
    assert dottest(Gop, par['ny']*5, par['nx']*5,
                   complexflag=0 if par['imag'] == 1 else 3)

    x = (np.ones((par['nx'], 5)) +
         par['imag'] * np.ones((par['nx'], 5))).flatten()
    xlsqr = lsqr(Gop, Gop*x, damp=1e-20, iter_lim=300, show=0)[0]
    assert_array_almost_equal(x, xlsqr, decimal=4)
Example #21
Source File: test_lsqr.py    From Computable with MIT License 5 votes vote down vote up 
def test_gh_2466():
    row = np.array([0, 0])
    col = np.array([0, 1])
    val = np.array([1, -1])
    A = scipy.sparse.coo_matrix((val, (row, col)), shape=(1, 2))
    b = np.asarray([4])
    lsqr(A, b)
Example #22
Source File: test_lsqr.py    From Computable with MIT License 5 votes vote down vote up 
def test_basic():
    svx = np.linalg.solve(G, b)
    X = lsqr(G, b, show=show, atol=tol, btol=tol, iter_lim=maxit)
    xo = X[0]
    assert_(norm(svx - xo) < 1e-5)
Example #23
Source File: test_smoothing.py    From pylops with GNU Lesser General Public License v3.0 5 votes vote down vote up 
def test_Smoothing2D(par):
    """Dot-test for smoothing
    """
    # 2d kernel on 2d signal
    if par['dir'] < 2:
        D2op = Smoothing2D(nsmooth=(5, 5), dims=(par['ny'], par['nx']), dtype='float32')
        assert dottest(D2op, par['ny']*par['nx'], par['ny']*par['nx'])

        x = np.zeros((par['ny'], par['nx']))
        x[par['ny']//2, par['nx']//2] = 1.
        x = x.flatten()
        y = D2op*x
        xlsqr = lsqr(D2op, y, damp=1e-10, iter_lim=400, show=0)[0]
        assert_array_almost_equal(x, xlsqr, decimal=1)

    # 2d kernel on 3d signal
    D2op = Smoothing2D(nsmooth=(5, 5), dims=(par['nz'], par['ny'], par['nx']),
                       nodir=par['dir'], dtype='float32')
    assert dottest(D2op, par['nz'] * par['ny'] * par['nx'], par['nz'] * par['ny'] * par['nx'])

    x = np.zeros((par['nz'], par['ny'], par['nx']))
    x[par['nz']//2, par['ny']//2, par['nx']//2] = 1.
    x = x.flatten()
    y = D2op*x
    xlsqr = lsqr(D2op, y, damp=1e-10, iter_lim=400, show=0)[0]
    assert_array_almost_equal(x, xlsqr, decimal=1)
Example #24
Source File: test_smoothing.py    From pylops with GNU Lesser General Public License v3.0 5 votes vote down vote up 
def test_Smoothing1D(par):
    """Dot-test and inversion for smoothing
    """
    # 1d kernel on 1d signal
    D1op = Smoothing1D(nsmooth=5, dims=par['nx'], dtype='float32')
    assert dottest(D1op, par['nx'], par['nx'])

    x = np.random.normal(0, 1, par['nx'])
    y = D1op * x
    xlsqr = lsqr(D1op, y, damp=1e-10, iter_lim=100, show=0)[0]
    assert_array_almost_equal(x, xlsqr, decimal=3)

    # 1d kernel on 2d signal
    D1op = Smoothing1D(nsmooth=5, dims=(par['ny'], par['nx']), dir=par['dir'], dtype='float32')
    assert dottest(D1op, par['ny']*par['nx'], par['ny']*par['nx'])

    x = np.random.normal(0, 1, (par['ny'], par['nx'])).flatten()
    y = D1op*x
    xlsqr = lsqr(D1op, y, damp=1e-10, iter_lim=100, show=0)[0]
    assert_array_almost_equal(x, xlsqr, decimal=3)

    # 1d kernel on 3d signal
    D1op = Smoothing1D(nsmooth=5, dims=(par['nz'], par['ny'], par['nx']),
                       dir=par['dir'], dtype='float32')
    assert dottest(D1op, par['nz'] * par['ny'] * par['nx'], par['nz'] * par['ny'] * par['nx'])

    x = np.random.normal(0, 1, (par['nz'], par['ny'], par['nx'])).flatten()
    y = D1op*x
    xlsqr = lsqr(D1op, y, damp=1e-10, iter_lim=100, show=0)[0]
    assert_array_almost_equal(x, xlsqr, decimal=3)
Example #25
Source File: test_lsqr.py    From GraphicDesignPatternByPython with MIT License 5 votes vote down vote up 
def test_basic():
    b_copy = b.copy()
    X = lsqr(G, b, show=show, atol=tol, btol=tol, iter_lim=maxit)
    assert_(np.all(b_copy == b))

    svx = np.linalg.solve(G, b)
    xo = X[0]
    assert_(norm(svx - xo) < 1e-5)
Example #26
Source File: test_lsqr.py    From GraphicDesignPatternByPython with MIT License 5 votes vote down vote up 
def test_gh_2466():
    row = np.array([0, 0])
    col = np.array([0, 1])
    val = np.array([1, -1])
    A = scipy.sparse.coo_matrix((val, (row, col)), shape=(1, 2))
    b = np.asarray([4])
    lsqr(A, b)
Example #27
Source File: test_lsqr.py    From GraphicDesignPatternByPython with MIT License 5 votes vote down vote up 
def test_well_conditioned_problems():
    # Test that sparse the lsqr solver returns the right solution
    # on various problems with different random seeds.
    # This is a non-regression test for a potential ZeroDivisionError
    # raised when computing the `test2` & `test3` convergence conditions.
    n = 10
    A_sparse = scipy.sparse.eye(n, n)
    A_dense = A_sparse.toarray()

    with np.errstate(invalid='raise'):
        for seed in range(30):
            rng = np.random.RandomState(seed + 10)
            beta = rng.rand(n)
            beta[beta == 0] = 0.00001  # ensure that all the betas are not null
            b = A_sparse * beta[:, np.newaxis]
            output = lsqr(A_sparse, b, show=show)

            # Check that the termination condition corresponds to an approximate
            # solution to Ax = b
            assert_equal(output[1], 1)
            solution = output[0]

            # Check that we recover the ground truth solution
            assert_array_almost_equal(solution, beta)

            # Sanity check: compare to the dense array solver
            reference_solution = np.linalg.solve(A_dense, b).ravel()
            assert_array_almost_equal(solution, reference_solution)
Example #28
Source File: test_lsqr.py    From GraphicDesignPatternByPython with MIT License 5 votes vote down vote up 
def test_b_shapes():
    # Test b being a scalar.
    A = np.array([[1.0, 2.0]])
    b = 3.0
    x = lsqr(A, b)[0]
    assert_almost_equal(norm(A.dot(x) - b), 0)

    # Test b being a column vector.
    A = np.eye(10)
    b = np.ones((10, 1))
    x = lsqr(A, b)[0]
    assert_almost_equal(norm(A.dot(x) - b.ravel()), 0)
Example #29
Source File: test_lsqr.py    From GraphicDesignPatternByPython with MIT License 5 votes vote down vote up 
def test_initialization():
    # Test the default setting is the same as zeros
    b_copy = b.copy()
    x_ref = lsqr(G, b, show=show, atol=tol, btol=tol, iter_lim=maxit)
    x0 = np.zeros(x_ref[0].shape)
    x = lsqr(G, b, show=show, atol=tol, btol=tol, iter_lim=maxit, x0=x0)
    assert_(np.all(b_copy == b))
    assert_array_almost_equal(x_ref[0], x[0])

    # Test warm-start with single iteration
    x0 = lsqr(G, b, show=show, atol=tol, btol=tol, iter_lim=1)[0]
    x = lsqr(G, b, show=show, atol=tol, btol=tol, iter_lim=maxit, x0=x0)
    assert_array_almost_equal(x_ref[0], x[0])
    assert_(np.all(b_copy == b))
Example #30
Source File: mesh_edit.py    From hmd with MIT License 5 votes vote down vote up 
def solveLaplacianMesh(mesh, anchors, anchorsIdx, cotangent=True):
    n = mesh.n_vertices()
    k = anchorsIdx.shape[0]

    operator = (getLaplacianMatrixUmbrella, getLaplacianMatrixCotangent)

    L = operator[1](mesh, anchorsIdx) if cotangent else operator[0](mesh, anchorsIdx)
    delta = np.array(L.dot(mesh.points()))

    # augment delta solution matrix with weighted anchors
    for i in range(k):
        delta[n + i, :] = WEIGHT * anchors[i, :]

    # update mesh vertices with least-squares solution
    for i in range(3):
        #mesh.points()[:, i] = lsqr(L, delta[:, i])[0]
        mesh.points()[:, i] = sparseqr.solve(L, delta[:, i], tolerance = 1e-8)
    
    return mesh



##############################################################
##           High Speed Laplacian Mesh Editing              ##
##############################################################
# using umbrella weights for higher speed