lobpcg.py

"""
Pure SciPy implementation of Locally Optimal Block Preconditioned Conjugate
Gradient Method (LOBPCG), see
http://www-math.cudenver.edu/~aknyazev/software/BLOPEX/

License: BSD

Authors: Robert Cimrman, Andrew Knyazev

Examples in tests directory contributed by Nils Wagner.
"""

import sys
import numpy as np
import scipy as sp

from scipy.sparse.linalg import aslinearoperator, LinearOperator

__all__ = ['lobpcg']

## try:
##     from symeig import symeig
## except:
##     raise ImportError('lobpcg requires symeig')

def symeig( mtxA, mtxB = None, eigenvectors = True, select = None ):
    import scipy.linalg as sla
    import scipy.lib.lapack as ll
    if select is None:
        if np.iscomplexobj( mtxA ):
            if mtxB is None:
                fun = ll.get_lapack_funcs( ['heev'], arrays = (mtxA,) )[0]
            else:
                fun = ll.get_lapack_funcs( ['hegv'], arrays = (mtxA,) )[0]
        else:
            if mtxB is None:
                fun = ll.get_lapack_funcs( ['syev'], arrays = (mtxA,) )[0]
            else:
                fun = ll.get_lapack_funcs( ['sygv'], arrays = (mtxA,) )[0]
##         print fun
        if mtxB is None:
            out = fun( mtxA )
        else:
            out = fun( mtxA, mtxB )
##         print w
##         print v
##         print info
##         from symeig import symeig
##         print symeig( mtxA, mtxB )
    else:
        out = sla.eig( mtxA, mtxB, right = eigenvectors )
        w = out[0]
        ii = np.argsort( w )
        w = w[slice( *select )]
        if eigenvectors:
            v = out[1][:,ii]
            v = v[:,slice( *select )]
            out = w, v, 0
        else:
            out = w, 0

    return out[:-1]

def pause():
    raw_input()

def save( ar, fileName ):
    from numpy import savetxt
    savetxt( fileName, ar, precision = 8 )

##
# 21.05.2007, c
def as2d( ar ):
    """
    If the input array is 2D return it, if it is 1D, append a dimension,
    making it a column vector.
    """
    if ar.ndim == 2:
        return ar
    else: # Assume 1!
        aux = np.array( ar, copy = False )
        aux.shape = (ar.shape[0], 1)
        return aux

class CallableLinearOperator(LinearOperator):
    def __call__(self, x):
        return self.matmat(x)

def makeOperator( operatorInput, expectedShape ):
    """Internal. Takes a dense numpy array or a sparse matrix or
    a function and makes an operator performing matrix * blockvector
    products.

    Examples
    --------
    >>> A = makeOperator( arrayA, (n, n) )
    >>> vectorB = A( vectorX )

    """
    if operatorInput is None:
        def ident(x):
            return x
        operator = LinearOperator(expectedShape, ident, matmat=ident)
    else:
        operator = aslinearoperator(operatorInput)

    if operator.shape != expectedShape:
        raise ValueError('operator has invalid shape')

    if sys.version_info[0] >= 3:
        # special methods are looked up on the class -- so make a new one
        operator.__class__ = CallableLinearOperator
    else:
        operator.__call__ = operator.matmat

    return operator



def applyConstraints( blockVectorV, factYBY, blockVectorBY, blockVectorY ):
    """Internal. Changes blockVectorV in place."""
    gramYBV = sp.dot( blockVectorBY.T, blockVectorV )
    import scipy.linalg as sla
    tmp = sla.cho_solve( factYBY, gramYBV )
    blockVectorV -= sp.dot( blockVectorY, tmp )


def b_orthonormalize( B, blockVectorV,
                      blockVectorBV = None, retInvR = False ):
    """Internal."""
    import scipy.linalg as sla
    if blockVectorBV is None:
        if B is not None:
            blockVectorBV = B( blockVectorV )
        else:
            blockVectorBV = blockVectorV # Shared data!!!
    gramVBV = sp.dot( blockVectorV.T, blockVectorBV )
    gramVBV = sla.cholesky( gramVBV )
    gramVBV = sla.inv( gramVBV, overwrite_a = True )
    # gramVBV is now R^{-1}.
    blockVectorV = sp.dot( blockVectorV, gramVBV )
    if B is not None:
        blockVectorBV = sp.dot( blockVectorBV, gramVBV )

    if retInvR:
        return blockVectorV, blockVectorBV, gramVBV
    else:
        return blockVectorV, blockVectorBV

def lobpcg( A, X,
            B=None, M=None, Y=None,
            tol= None, maxiter=20,
            largest = True, verbosityLevel = 0,
            retLambdaHistory = False, retResidualNormsHistory = False ):
    """Solve symmetric partial eigenproblems with optional preconditioning

    This function implements the Locally Optimal Block Preconditioned
    Conjugate Gradient Method (LOBPCG).

    Parameters
    ----------
    A : {sparse matrix, dense matrix, LinearOperator}
        The symmetric linear operator of the problem, usually a
        sparse matrix.  Often called the "stiffness matrix".
    X : array_like
        Initial approximation to the k eigenvectors. If A has
        shape=(n,n) then X should have shape shape=(n,k).
    B : {dense matrix, sparse matrix, LinearOperator}, optional
        the right hand side operator in a generalized eigenproblem.
        by default, B = Identity
        often called the "mass matrix"
    M : {dense matrix, sparse matrix, LinearOperator}, optional
        preconditioner to A; by default M = Identity
        M should approximate the inverse of A
    Y : array_like, optional
        n-by-sizeY matrix of constraints, sizeY < n
        The iterations will be performed in the B-orthogonal complement
        of the column-space of Y. Y must be full rank.

    Returns
    -------
    w : array
        Array of k eigenvalues
    v : array
        An array of k eigenvectors.  V has the same shape as X.

    Other Parameters
    ----------------
    tol : scalar, optional
        Solver tolerance (stopping criterion)
        by default: tol=n*sqrt(eps)
    maxiter: integer, optional
        maximum number of iterations
        by default: maxiter=min(n,20)
    largest : boolean, optional
        when True, solve for the largest eigenvalues, otherwise the smallest
    verbosityLevel : integer, optional
        controls solver output.  default: verbosityLevel = 0.
    retLambdaHistory : boolean, optional
        whether to return eigenvalue history
    retResidualNormsHistory : boolean, optional
        whether to return history of residual norms


    Notes
    -----
    If both retLambdaHistory and retResidualNormsHistory are True, the
    return tuple has the following format
    (lambda, V, lambda history, residual norms history)

    """
    failureFlag = True
    import scipy.linalg as sla

    blockVectorX = X
    blockVectorY = Y
    residualTolerance = tol
    maxIterations = maxiter

    if blockVectorY is not None:
        sizeY = blockVectorY.shape[1]
    else:
        sizeY = 0

    # Block size.
    if len(blockVectorX.shape) != 2:
        raise ValueError('expected rank-2 array for argument X')

    n, sizeX = blockVectorX.shape
    if sizeX > n:
        raise ValueError('X column dimension exceeds the row dimension')

    A = makeOperator(A, (n,n))
    B = makeOperator(B, (n,n))
    M = makeOperator(M, (n,n))

    if (n - sizeY) < (5 * sizeX):
        #warn('The problem size is small compared to the block size.' \
        #        ' Using dense eigensolver instead of LOBPCG.')

        if blockVectorY is not None:
            raise NotImplementedError('symeig does not support constraints')

        if largest:
            lohi = (n - sizeX, n)
        else:
            lohi = (1, sizeX)

        A_dense = A(np.eye(n))

        if B is not None:
            B_dense = B(np.eye(n))
            _lambda, eigBlockVector = symeig(A_dense, B_dense, select=lohi )
        else:
            _lambda, eigBlockVector = symeig(A_dense, select=lohi )

        return _lambda, eigBlockVector


    if residualTolerance is None:
        residualTolerance = np.sqrt( 1e-15 ) * n

    maxIterations = min( n, maxIterations )

    if verbosityLevel:
        aux = "Solving "
        if B is None:
            aux += "standard"
        else:
            aux += "generalized"
        aux += " eigenvalue problem with"
        if M is None:
            aux += "out"
        aux += " preconditioning\n\n"
        aux += "matrix size %d\n" % n
        aux += "block size %d\n\n" % sizeX
        if blockVectorY is None:
            aux += "No constraints\n\n"
        else:
            if sizeY > 1:
                aux += "%d constraints\n\n" % sizeY
            else:
                aux += "%d constraint\n\n" % sizeY
        print aux

    ##
    # Apply constraints to X.
    if blockVectorY is not None:

        if B is not None:
            blockVectorBY = B( blockVectorY )
        else:
            blockVectorBY = blockVectorY

        # gramYBY is a dense array.
        gramYBY = sp.dot( blockVectorY.T, blockVectorBY )
        try:
            # gramYBY is a Cholesky factor from now on...
            gramYBY = sla.cho_factor( gramYBY )
        except:
            raise ValueError('cannot handle linearly dependent constraints')

        applyConstraints( blockVectorX, gramYBY, blockVectorBY, blockVectorY )

    ##
    # B-orthonormalize X.
    blockVectorX, blockVectorBX = b_orthonormalize( B, blockVectorX )

    ##
    # Compute the initial Ritz vectors: solve the eigenproblem.
    blockVectorAX = A( blockVectorX )
    gramXAX = sp.dot( blockVectorX.T, blockVectorAX )
    # gramXBX is X^T * X.
    gramXBX = sp.dot( blockVectorX.T, blockVectorX )

    _lambda, eigBlockVector = symeig( gramXAX )
    ii = np.argsort( _lambda )[:sizeX]
    if largest:
        ii = ii[::-1]
    _lambda = _lambda[ii]

    eigBlockVector = np.asarray( eigBlockVector[:,ii] )
    blockVectorX  = sp.dot( blockVectorX,  eigBlockVector )
    blockVectorAX = sp.dot( blockVectorAX, eigBlockVector )
    if B is not None:
        blockVectorBX = sp.dot( blockVectorBX, eigBlockVector )

    ##
    # Active index set.
    activeMask = np.ones( (sizeX,), dtype = np.bool )

    lambdaHistory = [_lambda]
    residualNormsHistory = []

    previousBlockSize = sizeX
    ident  = np.eye( sizeX, dtype = A.dtype )
    ident0 = np.eye( sizeX, dtype = A.dtype )

    ##
    # Main iteration loop.
    for iterationNumber in xrange( maxIterations ):
        if verbosityLevel > 0:
            print 'iteration %d' %  iterationNumber

        aux = blockVectorBX * _lambda[np.newaxis,:]
        blockVectorR = blockVectorAX - aux

        aux = np.sum( blockVectorR.conjugate() * blockVectorR, 0 )
        residualNorms = np.sqrt( aux )

        residualNormsHistory.append( residualNorms )

        ii = np.where( residualNorms > residualTolerance, True, False )
        activeMask = activeMask & ii
        if verbosityLevel > 2:
            print activeMask

        currentBlockSize = activeMask.sum()
        if currentBlockSize != previousBlockSize:
            previousBlockSize = currentBlockSize
            ident = np.eye( currentBlockSize, dtype = A.dtype )

        if currentBlockSize == 0:
            failureFlag = False # All eigenpairs converged.
            break

        if verbosityLevel > 0:
            print 'current block size:', currentBlockSize
            print 'eigenvalue:', _lambda
            print 'residual norms:', residualNorms
        if verbosityLevel > 10:
            print eigBlockVector

        activeBlockVectorR = as2d( blockVectorR[:,activeMask] )

        if iterationNumber > 0:
            activeBlockVectorP  = as2d( blockVectorP [:,activeMask] )
            activeBlockVectorAP = as2d( blockVectorAP[:,activeMask] )
            activeBlockVectorBP = as2d( blockVectorBP[:,activeMask] )

        if M is not None:
            # Apply preconditioner T to the active residuals.
            activeBlockVectorR = M( activeBlockVectorR )

        ##
        # Apply constraints to the preconditioned residuals.
        if blockVectorY is not None:
            applyConstraints( activeBlockVectorR,
                              gramYBY, blockVectorBY, blockVectorY )

        ##
        # B-orthonormalize the preconditioned residuals.

        aux = b_orthonormalize( B, activeBlockVectorR )
        activeBlockVectorR, activeBlockVectorBR = aux

        activeBlockVectorAR = A( activeBlockVectorR )

        if iterationNumber > 0:
            aux = b_orthonormalize( B, activeBlockVectorP,
                                    activeBlockVectorBP, retInvR = True )
            activeBlockVectorP, activeBlockVectorBP, invR = aux
            activeBlockVectorAP = sp.dot( activeBlockVectorAP, invR )

        ##
        # Perform the Rayleigh Ritz Procedure:
        # Compute symmetric Gram matrices:

        xaw = sp.dot( blockVectorX.T,       activeBlockVectorAR )
        waw = sp.dot( activeBlockVectorR.T, activeBlockVectorAR )
        xbw = sp.dot( blockVectorX.T,       activeBlockVectorBR )

        if iterationNumber > 0:
            xap = sp.dot( blockVectorX.T,       activeBlockVectorAP )
            wap = sp.dot( activeBlockVectorR.T, activeBlockVectorAP )
            pap = sp.dot( activeBlockVectorP.T, activeBlockVectorAP )
            xbp = sp.dot( blockVectorX.T,       activeBlockVectorBP )
            wbp = sp.dot( activeBlockVectorR.T, activeBlockVectorBP )

            gramA = np.bmat( [[np.diag( _lambda ),   xaw,  xap],
                              [             xaw.T,   waw,  wap],
                              [             xap.T, wap.T,  pap]] )

            gramB = np.bmat( [[ident0,    xbw,    xbp],
                              [ xbw.T,  ident,    wbp],
                              [ xbp.T,  wbp.T,  ident]] )
        else:
            gramA = np.bmat( [[np.diag( _lambda ),  xaw],
                              [             xaw.T,  waw]] )
            gramB = np.bmat( [[ident0,    xbw],
                              [ xbw.T,  ident]] )

        try:
            assert np.allclose( gramA.T, gramA )
        except:
            print gramA.T - gramA
            raise

        try:
            assert np.allclose( gramB.T, gramB )
        except:
            print gramB.T - gramB
            raise

        if verbosityLevel > 10:
            save( gramA, 'gramA' )
            save( gramB, 'gramB' )

        ##
        # Solve the generalized eigenvalue problem.
#        _lambda, eigBlockVector = la.eig( gramA, gramB )
        _lambda, eigBlockVector = symeig( gramA, gramB )
        ii = np.argsort( _lambda )[:sizeX]
        if largest:
            ii = ii[::-1]
        if verbosityLevel > 10:
            print ii

        _lambda = _lambda[ii].astype( np.float64 )
        eigBlockVector = np.asarray( eigBlockVector[:,ii].astype( np.float64 ) )

        lambdaHistory.append( _lambda )

        if verbosityLevel > 10:
            print 'lambda:', _lambda
##         # Normalize eigenvectors!
##         aux = np.sum( eigBlockVector.conjugate() * eigBlockVector, 0 )
##         eigVecNorms = np.sqrt( aux )
##         eigBlockVector = eigBlockVector / eigVecNorms[np.newaxis,:]
#        eigBlockVector, aux = b_orthonormalize( B, eigBlockVector )

        if verbosityLevel > 10:
            print eigBlockVector
            pause()

        ##
        # Compute Ritz vectors.
        if iterationNumber > 0:
            eigBlockVectorX = eigBlockVector[:sizeX]
            eigBlockVectorR = eigBlockVector[sizeX:sizeX+currentBlockSize]
            eigBlockVectorP = eigBlockVector[sizeX+currentBlockSize:]

            pp  = sp.dot( activeBlockVectorR, eigBlockVectorR )
            pp += sp.dot( activeBlockVectorP, eigBlockVectorP )

            app  = sp.dot( activeBlockVectorAR, eigBlockVectorR )
            app += sp.dot( activeBlockVectorAP, eigBlockVectorP )

            bpp  = sp.dot( activeBlockVectorBR, eigBlockVectorR )
            bpp += sp.dot( activeBlockVectorBP, eigBlockVectorP )
        else:
            eigBlockVectorX = eigBlockVector[:sizeX]
            eigBlockVectorR = eigBlockVector[sizeX:]

            pp  = sp.dot( activeBlockVectorR,  eigBlockVectorR )
            app = sp.dot( activeBlockVectorAR, eigBlockVectorR )
            bpp = sp.dot( activeBlockVectorBR, eigBlockVectorR )

        if verbosityLevel > 10:
            print pp
            print app
            print bpp
            pause()

        blockVectorX  = sp.dot( blockVectorX, eigBlockVectorX )  + pp
        blockVectorAX = sp.dot( blockVectorAX, eigBlockVectorX ) + app
        blockVectorBX = sp.dot( blockVectorBX, eigBlockVectorX ) + bpp

        blockVectorP, blockVectorAP, blockVectorBP = pp, app, bpp

    aux = blockVectorBX * _lambda[np.newaxis,:]
    blockVectorR = blockVectorAX - aux

    aux = np.sum( blockVectorR.conjugate() * blockVectorR, 0 )
    residualNorms = np.sqrt( aux )


    if verbosityLevel > 0:
        print 'final eigenvalue:', _lambda
        print 'final residual norms:', residualNorms

    if retLambdaHistory:
        if retResidualNormsHistory:
            return _lambda, blockVectorX, lambdaHistory, residualNormsHistory
        else:
            return _lambda, blockVectorX, lambdaHistory
    else:
        if retResidualNormsHistory:
            return _lambda, blockVectorX, residualNormsHistory
        else:
            return _lambda, blockVectorX

###########################################################################
if __name__ == '__main__':
    from scipy.sparse import spdiags, speye, issparse
    import time

##     def B( vec ):
##         return vec

    n = 100
    vals = [np.arange( n, dtype = np.float64 ) + 1]
    A = spdiags( vals, 0, n, n )
    B = speye( n, n )
#    B[0,0] = 0
    B = np.eye( n, n )
    Y = np.eye( n, 3 )


#    X = sp.rand( n, 3 )
    xfile = {100 : 'X.txt', 1000 : 'X2.txt', 10000 : 'X3.txt'}
    X = np.fromfile( xfile[n], dtype = np.float64, sep = ' ' )
    X.shape = (n, 3)

    ivals = [1./vals[0]]
    def precond( x ):
        invA = spdiags( ivals, 0, n, n )
        y = invA  * x
        if issparse( y ):
            y = y.toarray()

        return as2d( y )

    precond = spdiags( ivals, 0, n, n )
#    precond = None
    tt = time.clock()
#    B = None
    eigs, vecs = lobpcg( X, A, B, blockVectorY = Y,
                         M = precond,
                         residualTolerance = 1e-4, maxIterations = 40,
                         largest = False, verbosityLevel = 1 )
    print 'solution time:', time.clock() - tt

    print vecs
    print eigs