[Scipy-svn] r3877 - in trunk/scipy: sparse splinalg splinalg/isolve splinalg/tests

scipy-svn@scip... scipy-svn@scip...
Wed Jan 30 20:44:36 CST 2008


Author: wnbell
Date: 2008-01-30 20:44:28 -0600 (Wed, 30 Jan 2008)
New Revision: 3877

Added:
   trunk/scipy/splinalg/interface.py
   trunk/scipy/splinalg/isolve/minres.py
   trunk/scipy/splinalg/tests/
   trunk/scipy/splinalg/tests/test_interface.py
Modified:
   trunk/scipy/sparse/compressed.py
   trunk/scipy/splinalg/__init__.py
   trunk/scipy/splinalg/setup.py
Log:
added LinearOperator
added draft MINRES implementation


Modified: trunk/scipy/sparse/compressed.py
===================================================================
--- trunk/scipy/sparse/compressed.py	2008-01-30 17:05:38 UTC (rev 3876)
+++ trunk/scipy/sparse/compressed.py	2008-01-31 02:44:28 UTC (rev 3877)
@@ -8,7 +8,7 @@
 import numpy
 from numpy import array, matrix, asarray, asmatrix, zeros, rank, intc, \
         empty, hstack, isscalar, ndarray, shape, searchsorted, empty_like, \
-        where, concatenate
+        where, concatenate, transpose
 
 from base import spmatrix, isspmatrix, SparseEfficiencyWarning
 from data import _data_matrix
@@ -364,9 +364,9 @@
 
     def rmatvec(self, other, conjugate=True):
         if conjugate:
-            return self.transpose().conj() * other
+            return transpose( self.transpose().conj().matvec(transpose(other)) )
         else:
-            return self.transpose() * other
+            return transpose( self.transpose().matvec(transpose(other)) )
 
     def getdata(self, ind):
         return self.data[ind]

Modified: trunk/scipy/splinalg/__init__.py
===================================================================
--- trunk/scipy/splinalg/__init__.py	2008-01-30 17:05:38 UTC (rev 3876)
+++ trunk/scipy/splinalg/__init__.py	2008-01-31 02:44:28 UTC (rev 3877)
@@ -4,6 +4,7 @@
 
 from isolve import *
 from dsolve import *
+from interface import *
 
 __all__ = filter(lambda s:not s.startswith('_'),dir())
 from scipy.testing.pkgtester import Tester

Added: trunk/scipy/splinalg/interface.py
===================================================================
--- trunk/scipy/splinalg/interface.py	2008-01-30 17:05:38 UTC (rev 3876)
+++ trunk/scipy/splinalg/interface.py	2008-01-31 02:44:28 UTC (rev 3877)
@@ -0,0 +1,87 @@
+import numpy as np
+from scipy.sparse.sputils import isshape
+from scipy.sparse import isspmatrix
+
+__all__ = ['LinearOperator', 'aslinearoperator']
+
+class LinearOperator:
+    def __init__( self, shape, matvec, rmatvec=None, dtype=None ):
+        """Common interface for performing matrix vector products
+        """
+        shape = tuple(shape)
+
+        if not isshape(shape):
+            raise ValueError('invalid shape')
+
+        self.shape  = shape
+        self.matvec = matvec
+
+        if rmatvec is None:
+            def rmatvec(x):
+                raise NotImplementedError('rmatvec is not defined')
+            self.rmatvec = rmatvec
+        else:
+            self.rmatvec = rmatvec
+
+        if dtype is not None:
+            self.dtype = np.dtype(dtype)
+
+    def __repr__(self):
+        M,N = self.shape
+        if hasattr(self,'dtype'):
+            dt = 'dtype=' + str(self.dtype)
+        else:
+            dt = 'unspecified dtype'
+
+        return '<%dx%d LinearOperator with %s>' % (M,N,dt)
+
+def aslinearoperator(A):
+    """Return A as a LinearOperator
+
+    'A' may be any of the following types
+        - ndarray
+        - matrix
+        - sparse matrix (e.g. csr_matrix, lil_matrix, etc.)
+        - LinearOperator
+        - An object with .shape and .matvec attributes
+
+    See the LinearOperator documentation for additonal information.
+
+    Examples
+    ========
+    
+    >>> from scipy import matrix
+    >>> M = matrix( [[1,2,3],[4,5,6]], dtype='int32' )
+    >>> aslinearoperator( M )
+    <2x3 LinearOperator with dtype=int32>
+
+    """
+
+    if isinstance(A, LinearOperator):
+        return A
+
+    elif isinstance(A, np.ndarray) or isinstance(A,np.matrix):
+        def matvec(x):
+            return np.dot(np.asarray(A),x)
+        def rmatvec(x):
+            return np.dot(x,np.asarray(A))
+        return LinearOperator( A.shape, matvec, rmatvec=rmatvec, dtype=A.dtype )
+
+    elif isspmatrix(A):
+        return LinearOperator( A.shape, A.matvec, rmatvec=A.rmatvec, dtype=A.dtype ) 
+
+    else:
+        if hasattr(A,'shape') and hasattr(A,'matvec'):
+            rmatvec = None
+            dtype = None
+
+            if hasattr(A,'rmatvec'):
+                rmatvec = A.rmatvec
+            if hasattr(A,'dtype'):
+                dtype = A.dtype
+            return LinearOperator(A.shape, A.matvec, rmatvec=rmatvec, dtype=dtype)
+
+        else:
+            raise TypeError('type not understood')
+
+

Added: trunk/scipy/splinalg/isolve/minres.py
===================================================================
--- trunk/scipy/splinalg/isolve/minres.py	2008-01-30 17:05:38 UTC (rev 3876)
+++ trunk/scipy/splinalg/isolve/minres.py	2008-01-31 02:44:28 UTC (rev 3877)
@@ -0,0 +1,297 @@
+from numpy import sqrt, inner, finfo
+from numpy.linalg import norm
+
+def psolve(x): return x
+def check_sizes(A,x,b): pass
+
+def minres(A, b, x0=None, shift=0.0, tol=1e-5, maxiter=None, xtype=None,
+        precond=None, callback=None, show=False, check=False):
+    """Use the Minimum Residual Method (MINRES) to solve Ax=b 
+    
+    MINRES minimizes norm(A*x - b) for the symmetric matrix A.  Unlike
+    the Conjugate Gradient method, A can be indefinite or singular.
+
+    If shift != 0 then the method solves (A - shift*I)x = b
+
+
+    Parameters
+    ==========
+        TODO
+
+    References
+    ==========
+        
+        Solution of sparse indefinite systems of linear equations,
+            C. C. Paige and M. A. Saunders (1975),
+            SIAM J. Numer. Anal. 12(4), pp. 617-629.
+            http://www.stanford.edu/group/SOL/software/minres.html
+
+        This file is a translation of the following MATLAB implementation:
+            http://www.stanford.edu/group/SOL/software/minres/matlab/
+
+    """ 
+
+    show  = True  #TODO remove
+    check = True  #TODO remove
+
+    first = 'Enter minres.   '
+    last  = 'Exit  minres.   '
+
+    assert(A.shape[0] == A.shape[1])
+    assert(A.shape[1] == len(b))
+
+    b = asarray(b).ravel()
+    n = A.shape[0]
+
+    if maxiter is None:
+        maxiter = 5 * n
+
+    matvec = A.matvec
+
+    msg   =[' beta2 = 0.  If M = I, b and x are eigenvectors    ',   # -1
+            ' beta1 = 0.  The exact solution is  x = 0          ',   #  0
+            ' A solution to Ax = b was found, given rtol        ',   #  1
+            ' A least-squares solution was found, given rtol    ',   #  2
+            ' Reasonable accuracy achieved, given eps           ',   #  3
+            ' x has converged to an eigenvector                 ',   #  4
+            ' acond has exceeded 0.1/eps                        ',   #  5
+            ' The iteration limit was reached                   ',   #  6
+            ' Aname  does not define a symmetric matrix         ',   #  7
+            ' Mname  does not define a symmetric matrix         ',   #  8
+            ' Mname  does not define a pos-def preconditioner   ']   #  9
+
+     
+    if show:
+        print first + 'Solution of symmetric Ax = b'
+        print first + 'n      =  %3g     shift  =  %23.14e'  % (n,shift)
+        print first + 'itnlim =  %3g     rtol   =  %11.2e'   % (maxiter,tol)
+        print
+        
+    istop = 0;   itn   = 0;   Anorm = 0;    Acond = 0;
+    rnorm = 0;   ynorm = 0; 
+
+    xtype = A.dtype #TODO update
+
+    eps = finfo(xtype).eps 
+
+    x = zeros( n, dtype=xtype )
+
+    # Set up y and v for the first Lanczos vector v1.
+    # y  =  beta1 P' v1,  where  P = C**(-1).
+    # v is really P' v1.
+    
+    y  = b
+    r1 = b
+
+    y = psolve(b)
+    
+    beta1 = inner(b,y)
+
+    if beta1 < 0:
+        raise ValueError('indefinite preconditioner')
+    elif beta1 == 0:
+        return x
+    
+    beta1 = sqrt( beta1 )
+
+    if check:
+        # see if M is symmetric
+        r2   = psolve(y)
+        s    = inner(y,y)
+        t    = inner(r1,r2)
+        z    = abs( s - t )
+        epsa = (s + eps) * eps**(1.0/3.0)
+        if z > epsa:
+            raise ValueError('non-symmetric preconditioner')
+        
+        # see if A is symmetric
+        w    = matvec(y)
+        r2   = matvec(w)
+        s    = inner(w,w)
+        t    = inner(y,r2)
+        epsa = (s + eps) * eps**(1.0/3.0)
+        if z > epsa:
+            raise ValueError('non-symmetric matrix')
+
+
+    # Initialize other quantities
+    oldb   = 0;          beta   = beta1;   dbar   = 0;       epsln  = 0;
+    qrnorm = beta1;      phibar = beta1;   rhs1   = beta1;
+    rhs2   = 0;          tnorm2 = 0;       ynorm2 = 0;
+    cs     = -1;         sn     = 0;
+    w      = zeros(n, dtype=xtype)
+    w2     = zeros(n, dtype=xtype)
+    r2     = r1
+
+    if show:
+      print
+      print
+      print '   Itn     x(1)     Compatible    LS       norm(A)  cond(A) gbar/|A|'
+
+    while itn < maxiter:
+        itn += 1
+
+        s = 1.0/beta 
+        v = s*y
+        
+        y  = matvec(v)
+        y  = y - shift * v
+
+        if itn >= 2:
+            y  = y - (beta/oldb)*r1
+        
+        alfa   = inner(v,y)
+        y      = y - (alfa/beta)*r2
+        r1     = r2
+        r2     = y
+        y      = psolve(r2)
+        oldb   = beta
+        beta   = inner(r2,y)
+        if beta < 0:
+            raise ValueError('non-symmetric matrix')
+        beta    = sqrt(beta) 
+        tnorm2 += alfa**2 + oldb**2 + beta**2
+
+        if itn == 1:
+            if beta/beta1 <= 10*eps: 
+                istop = -1  # Terminate later
+            #tnorm2 = alfa**2 ??
+            gmax = abs(alfa)
+            gmin = gmax
+
+        # Apply previous rotation Qk-1 to get
+        #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
+        #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
+
+        oldeps = epsln
+        delta  = cs * dbar  +  sn * alfa   # delta1 = 0         deltak
+        gbar   = sn * dbar  -  cs * alfa   # gbar 1 = alfa1     gbar k
+        epsln  =               sn * beta   # epsln2 = 0         epslnk+1
+        dbar   =            -  cs * beta   # dbar 2 = beta2     dbar k+1
+        root   = norm([gbar, dbar])
+        Arnorm = phibar * root
+
+        # Compute the next plane rotation Qk
+        
+        gamma  = norm([gbar, beta])       # gammak
+        gamma  = max(gamma, eps) 
+        cs     = gbar / gamma             # ck
+        sn     = beta / gamma             # sk
+        phi    = cs * phibar              # phik
+        phibar = sn * phibar              # phibark+1
+        
+        # Update  x.
+        
+        denom = 1.0/gamma
+        w1    = w2
+        w2    = w
+        w     = (v - oldeps*w1 - delta*w2) * denom
+        x     = x + phi*w
+        
+        # Go round again.
+        
+        gmax   = max(gmax, gamma)
+        gmin   = min(gmin, gamma)
+        z      = rhs1 / gamma
+        ynorm2 = z**2  +  ynorm2
+        rhs1   = rhs2 -  delta*z
+        rhs2   =      -  epsln*z
+        
+        # Estimate various norms and test for convergence.
+        
+        Anorm  = sqrt( tnorm2 )
+        ynorm  = sqrt( ynorm2 )
+        epsa   = Anorm * eps
+        epsx   = Anorm * ynorm * eps
+        epsr   = Anorm * ynorm * tol
+        diag   = gbar
+
+        if diag == 0: diag = epsa
+        
+        qrnorm = phibar
+        rnorm  = qrnorm
+        test1  = rnorm / (Anorm*ynorm)    #  ||r|| / (||A|| ||x||)
+        test2  = root  /  Anorm           # ||Ar|| / (||A|| ||r||)
+
+        # Estimate  cond(A).
+        # In this version we look at the diagonals of  R  in the
+        # factorization of the lower Hessenberg matrix,  Q * H = R,
+        # where H is the tridiagonal matrix from Lanczos with one
+        # extra row, beta(k+1) e_k^T.
+        
+        Acond  = gmax/gmin
+        
+        # See if any of the stopping criteria are satisfied.
+        # In rare cases, istop is already -1 from above (Abar = const*I).
+        
+        if istop == 0:
+            t1 = 1 + test1      # These tests work if tol < eps
+            t2 = 1 + test2
+            if t2    <= 1       : istop = 2
+            if t1    <= 1       : istop = 1
+            
+            if itn   >= maxiter : istop = 6
+            if Acond >= 0.1/eps : istop = 4
+            if epsx  >= beta1   : istop = 3
+            #if rnorm <= epsx   : istop = 2
+            #if rnorm <= epsr   : istop = 1
+            if test2 <= tol     : istop = 2
+            if test1 <= tol     : istop = 1
+    
+        # See if it is time to print something.
+        
+        prnt = False
+        if n        <= 40         : prnt = True
+        if itn      <= 10         : prnt = True
+        if itn      >= maxiter-10 : prnt = True
+        if itn % 10 == 0          : prnt = True
+        if qrnorm   <= 10*epsx    : prnt = True
+        if qrnorm   <= 10*epsr    : prnt = True
+        if Acond    <= 1e-2/eps   : prnt = True
+        if istop  !=  0           : prnt = True
+        
+        if show and prnt:
+            str1 = '%6g %12.5e %10.3e'  % (itn, x[0], test1)
+            str2 = ' %10.3e'            % (test2,)
+            str3 = ' %8.1e %8.1e %8.1e' % (Anorm, Acond, gbar/Anorm)
+
+            print str1 + str2 + str3
+        
+            if itn % 10 == 0: print
+
+        if callback is not None:
+            callback(x)
+
+        if istop > 0: break
+        
+
+    if show:
+        print
+        print last + ' istop   =  %3g               itn   =%5g' % (istop,itn) 
+        print last + ' Anorm   =  %12.4e      Acond =  %12.4e'  % (Anorm,Acond)
+        print last + ' rnorm   =  %12.4e      ynorm =  %12.4e'  % (rnorm,ynorm)
+        print last + ' Arnorm  =  %12.4e'                       %  (Arnorm,)
+        print last + msg[istop+1]
+
+    return x
+
+
+if __name__ == '__main__':
+    from scipy import *
+    from scipy.sparse import *
+    from scipy.splinalg import *
+    from scipy.sandbox.multigrid import *
+
+    n = 100
+
+    residuals = []
+
+    def cb(x):
+        residuals.append(norm(b - A*x))
+
+    #A = poisson((10,),format='csr')
+    A = spdiags( [arange(1,n+1,dtype=float)], [0], n, n, format='csr')
+    b = ones( A.shape[0] )
+    #x = minres(A,b,tol=1e-12,maxiter=None,callback=cb)
+    x = cg(A,b,x0=b,tol=1e-12,maxiter=None,callback=cb)[0]
+

Modified: trunk/scipy/splinalg/setup.py
===================================================================
--- trunk/scipy/splinalg/setup.py	2008-01-30 17:05:38 UTC (rev 3876)
+++ trunk/scipy/splinalg/setup.py	2008-01-31 02:44:28 UTC (rev 3877)
@@ -8,6 +8,8 @@
     config.add_subpackage(('isolve'))
     config.add_subpackage(('dsolve'))
     
+    config.add_data_dir('tests')
+    
     return config
 
 if __name__ == '__main__':

Added: trunk/scipy/splinalg/tests/test_interface.py
===================================================================
--- trunk/scipy/splinalg/tests/test_interface.py	2008-01-30 17:05:38 UTC (rev 3876)
+++ trunk/scipy/splinalg/tests/test_interface.py	2008-01-31 02:44:28 UTC (rev 3877)
@@ -0,0 +1,64 @@
+#!/usr/bin/env python
+""" Test functions for the splinalg.interface module
+"""
+
+from scipy.testing import *
+
+import numpy
+from numpy import array, matrix, ones, ravel
+from scipy.sparse import csr_matrix
+
+from scipy.splinalg.interface import *
+
+
+class TestInterface(TestCase):
+    def test_aslinearoperator(self):
+        cases = []
+
+        cases.append( matrix([[1,2,3],[4,5,6]]) )
+        cases.append( array([[1,2,3],[4,5,6]]) )
+        cases.append( csr_matrix([[1,2,3],[4,5,6]]) )
+
+        class matlike:
+            def __init__(self):
+                self.dtype = numpy.dtype('int')
+                self.shape = (2,3)
+            def matvec(self,x):
+                y = array([ 1*x[0] + 2*x[1] + 3*x[2],
+                            4*x[0] + 5*x[1] + 6*x[2]])
+                if len(x.shape) == 2:
+                    y = y.reshape(-1,1)
+                return y
+
+            def rmatvec(self,x):
+                if len(x.shape) == 1:
+                    y = array([ 1*x[0] + 4*x[1],
+                                2*x[0] + 5*x[1],
+                                3*x[0] + 6*x[1]])
+                    return y
+                else:
+                    y = array([ 1*x[0,0] + 4*x[0,1],
+                                2*x[0,0] + 5*x[0,1],
+                                3*x[0,0] + 6*x[0,1]])
+                    return y.reshape(1,-1)
+
+                return y
+               
+        cases.append( matlike() )
+
+
+        for M in cases:
+            A = aslinearoperator(M)
+            M,N = A.shape
+
+            assert_equal(A.matvec(array([1,2,3])),      [14,32])
+            assert_equal(A.matvec(array([[1],[2],[3]])),[[14],[32]])
+
+            assert_equal(A.rmatvec(array([1,2])),  [9,12,15])
+            assert_equal(A.rmatvec(array([[1,2]])),[[9,12,15]])
+
+            if hasattr(M,'dtype'):
+                assert_equal(A.dtype, M.dtype)
+
+if __name__ == "__main__":
+    nose.run(argv=['', __file__])



More information about the Scipy-svn mailing list