[Scipy-svn] r3039 - in trunk/Lib/sandbox/lobpcg: . tests
scipy-svn@scip...
scipy-svn@scip...
Thu May 24 07:30:14 CDT 2007
Author: rc
Date: 2007-05-24 07:30:05 -0500 (Thu, 24 May 2007)
New Revision: 3039
Added:
trunk/Lib/sandbox/lobpcg/tests/
trunk/Lib/sandbox/lobpcg/tests/benchmark.py
trunk/Lib/sandbox/lobpcg/tests/test_lobpcg.py
Modified:
trunk/Lib/sandbox/lobpcg/info.py
trunk/Lib/sandbox/lobpcg/lobpcg.py
trunk/Lib/sandbox/lobpcg/setup.py
Log:
more docs, examples and return values
Modified: trunk/Lib/sandbox/lobpcg/info.py
===================================================================
--- trunk/Lib/sandbox/lobpcg/info.py 2007-05-24 08:08:11 UTC (rev 3038)
+++ trunk/Lib/sandbox/lobpcg/info.py 2007-05-24 12:30:05 UTC (rev 3039)
@@ -1,16 +1,88 @@
"""
The algorithm of LOBPCG is described in detail in:
-A. V. Knyazev, Toward the Optimal Preconditioned Eigensolver: Locally Optimal Block Preconditioned Conjugate Gradient Method. SIAM Journal on Scientific Computing 23 (2001), no. 2, pp. 517-541. http://dx.doi.org/10.1137/S1064827500366124
+A. V. Knyazev, Toward the Optimal Preconditioned Eigensolver: Locally Optimal
+Block Preconditioned Conjugate Gradient Method. SIAM Journal on Scientific
+Computing 23 (2001), no. 2,
+pp. 517-541. http://dx.doi.org/10.1137/S1064827500366124
-A. V. Knyazev, I. Lashuk, M. E. Argentati, and E. Ovchinnikov, Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) in hypre and PETSc (2007). http://arxiv.org/abs/0705.2626
+A. V. Knyazev, I. Lashuk, M. E. Argentati, and E. Ovchinnikov, Block Locally
+Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) in hypre and PETSc
+(2007). http://arxiv.org/abs/0705.2626
+Call the function lobpcg - see help for lobpcg.lobpcg. See also lobpcg.as2d,
+which can be used in the preconditioner (example below)
-Depends upon symeig (http://mdp-toolkit.sourceforge.net/symeig.html) for the
-moment, as the symmetric eigenvalue solvers were not available in scipy.
+Example:
-Usage: XXXXX
+ # Solve A x = lambda B x with constraints and preconditioning.
+ n = 100
+ vals = [nm.arange( n, dtype = nm.float64 ) + 1]
+
+ # Matrix A.
+ operatorA = spdiags( vals, 0, n, n )
+
+ # Matrix B
+ operatorB = nm.eye( n, n )
+
+ # Constraints.
+ Y = nm.eye( n, 3 )
+
+ # Initial guess for eigenvectors, should have linearly independent
+ # columns. Column dimension = number of requested eigenvalues.
+ X = sc.rand( n, 3 )
+
+ # Preconditioner - inverse of A.
+ ivals = [1./vals[0]]
+ def precond( x ):
+ invA = spdiags( ivals, 0, n, n )
+ y = invA * x
+ if sp.issparse( y ):
+ y = y.toarray()
+
+ return as2d( y )
+
+ # Alternative way of providing the same preconditioner.
+ #precond = spdiags( ivals, 0, n, n )
+
+ tt = time.clock()
+ eigs, vecs = lobpcg( X, operatorA, operatorB, blockVectorY = Y,
+ operatorT = precond,
+ residualTolerance = 1e-4, maxIterations = 40,
+ largest = False, verbosityLevel = 1 )
+ print 'solution time:', time.clock() - tt
+ print eigs
+
+Usage notes:
+
+Notation: n - matrix size, m - number of required eigenvalues (smallest or
+largest)
+
+1) The LOBPCG code internally solves eigenproblems of the size 3m on every
+ iteration by calling the"standard" eigensolver, so if m is not small enough
+ compared to n, it does not make sense to call the LOBPCG code, but rather
+ one should use the "standard" eigensolver, e.g. symeig function in this
+ case. If one calls the LOBPCG algorithm for 5m>n, it will most likely break
+ internally, so the code tries to call symeig instead.
+
+ It is not that n should be large for the LOBPCG to work, but rather the
+ ratio n/m should be large. It you call the LOBPCG code with m=1 and n=10, it
+ should work, though n is small. The method is intended for extremely large
+ n/m, see e.g., reference [28] in http://arxiv.org/abs/0705.2626
+
+2) The convergence speed depends basically on two factors:
+
+ a) how well relatively separated the seeking eigenvalues are from the rest of
+ the eigenvalues. One can try to vary m to make this better.
+
+ b) how "well conditioned" the problem is. This can be changed by using proper
+ preconditioning. For example, a rod vibration test problem (under tests
+ directory) is ill-conditioned
+ for large n, so convergence will be slow, unless efficient preconditioning is
+ used. For this specific problem, a good simple preconditioner function would
+ be a linear solve for A, which is easy to code since A is tridiagonal.
+
"""
postpone_import = 1
Modified: trunk/Lib/sandbox/lobpcg/lobpcg.py
===================================================================
--- trunk/Lib/sandbox/lobpcg/lobpcg.py 2007-05-24 08:08:11 UTC (rev 3038)
+++ trunk/Lib/sandbox/lobpcg/lobpcg.py 2007-05-24 12:30:05 UTC (rev 3039)
@@ -5,7 +5,12 @@
License: BSD
+Depends upon symeig (http://mdp-toolkit.sourceforge.net/symeig.html) for the
+moment, as the symmetric eigenvalue solvers were not available in scipy.
+
(c) Robert Cimrman, Andrew Knyazev
+
+Examples in tests directory contributed by Nils Wagner.
"""
import numpy as nm
@@ -25,6 +30,10 @@
##
# 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!
@@ -35,11 +44,23 @@
##
# 05.04.2007, c
# 10.04.2007
+# 24.05.2007
def makeOperator( operatorInput, expectedShape ):
+ """
+ Internal. Takes a dense numpy array or a sparse matrix or a function and
+ makes an operator performing matrix * vector product.
+
+ :Example:
+
+ operatorA = makeOperator( arrayA, (n, n) )
+ vectorB = operatorA( vectorX )
+ """
class Operator( object ):
def __call__( self, vec ):
return self.call( vec )
-
+ def asMatrix( self ):
+ return self._asMatrix( self )
+
operator = Operator()
operator.obj = operatorInput
@@ -55,23 +76,30 @@
if sp.issparse( out ):
out = out.toarray()
return as2d( out )
+ def asMatrix( op ):
+ return op.obj.toarray()
else:
def call( vec ):
return as2d( nm.asarray( sc.dot( operator.obj, vec ) ) )
+ def asMatrix( op ):
+ return op.obj
operator.call = call
+ operator._asMatrix = asMatrix
+ operator.kind = 'matrix'
elif isinstance( operatorInput, types.FunctionType ) or \
isinstance( operatorInput, types.BuiltinFunctionType ):
operator.shape = expectedShape
operator.dtype = nm.float64
operator.call = operatorInput
+ operator.kind = 'function'
return operator
##
# 05.04.2007, c
def applyConstraints( blockVectorV, factYBY, blockVectorBY, blockVectorY ):
- """Changes blockVectorV in place."""
+ """Internal. Changes blockVectorV in place."""
gramYBV = sc.dot( blockVectorBY.T, blockVectorV )
tmp = la.cho_solve( factYBY, gramYBV )
blockVectorV -= sc.dot( blockVectorY, tmp )
@@ -80,7 +108,7 @@
# 05.04.2007, c
def b_orthonormalize( operatorB, blockVectorV,
blockVectorBV = None, retInvR = False ):
-
+ """Internal."""
if blockVectorBV is None:
if operatorB is not None:
blockVectorBV = operatorB( blockVectorV )
@@ -104,14 +132,68 @@
# 05.04.2007
# 06.04.2007
# 10.04.2007
+# 24.05.2007
def lobpcg( blockVectorX, operatorA,
operatorB = None, operatorT = None, blockVectorY = None,
residualTolerance = None, maxIterations = 20,
largest = True, verbosityLevel = 0,
retLambdaHistory = False, retResidualNormsHistory = False ):
+ """
+ LOBPCG solves symmetric partial eigenproblems using preconditioning.
- exitFlag = 0
+ Required input:
+ blockVectorX - initial approximation to eigenvectors, full or sparse matrix
+ n-by-blockSize
+
+ operatorA - the operator of the problem, can be given as a matrix or as an
+ M-file
+
+
+ Optional input:
+
+ operatorB - the second operator, if solving a generalized eigenproblem; by
+ default, or if empty, operatorB = I.
+
+ operatorT - preconditioner; by default, operatorT = I.
+
+
+ Optional constraints input:
+
+ blockVectorY - n-by-sizeY matrix of constraints, sizeY < n. The iterations
+ will be performed in the (operatorB-) orthogonal complement of the
+ column-space of blockVectorY. blockVectorY must be full rank.
+
+
+ Optional scalar input parameters:
+
+ residualTolerance - tolerance, by default, residualTolerance=n*sqrt(eps)
+
+ maxIterations - max number of iterations, by default, maxIterations =
+ min(n,20)
+
+ largest - when true, solve for the largest eigenvalues, otherwise for the
+ smallest
+
+ verbosityLevel - by default, verbosityLevel = 0.
+
+ retLambdaHistory - return eigenvalue history
+
+ retResidualNormsHistory - return history of residual norms
+
+ Output:
+
+ blockVectorX and lambda are computed blockSize eigenpairs, where
+ blockSize=size(blockVectorX,2) for the initial guess blockVectorX if it is
+ full rank.
+
+ If both retLambdaHistory and retResidualNormsHistory are True, the
+ return tuple has the flollowing order:
+
+ lambda, blockVectorX, lambda history, residual norms history
+ """
+ failureFlag = True
+
if blockVectorY is not None:
sizeY = blockVectorY.shape[1]
else:
@@ -133,6 +215,33 @@
if operatorB is not None:
operatorB = makeOperator( operatorB, (n, n) )
+ if (n - sizeY) < (5 * sizeX):
+ print 'The problem size is too small, compared to the block size, for LOBPCG to run.'
+ print 'Trying to use symeig instead, without preconditioning.'
+ if blockVectorY is not None:
+ print 'symeig does not support constraints'
+ raise ValueError
+
+ if largest:
+ lohi = (n - sizeX, n)
+ else:
+ lohi = (1, sizeX)
+
+ if operatorA.kind == 'function':
+ print 'symeig does not support matrix A given by function'
+
+ if operatorB is not None:
+ if operatorB.kind == 'function':
+ print 'symeig does not support matrix B given by function'
+
+ _lambda, eigBlockVector = symeig( operatorA.asMatrix(),
+ operatorB.asMatrix(),
+ range = lohi )
+ else:
+ _lambda, eigBlockVector = symeig( operatorA.asMatrix(),
+ range = lohi )
+ return _lambda, eigBlockVector
+
if operatorT is not None:
operatorT = makeOperator( operatorT, (n, n) )
## if n != operatorA.shape[0]:
@@ -254,9 +363,8 @@
previousBlockSize = currentBlockSize
ident = nm.eye( currentBlockSize, dtype = operatorA.dtype )
-
if currentBlockSize == 0:
- failureFlag = 0 # All eigenpairs converged.
+ failureFlag = False # All eigenpairs converged.
break
if verbosityLevel > 0:
@@ -377,6 +485,9 @@
_lambda = _lambda[ii].astype( nm.float64 )
eigBlockVector = nm.asarray( eigBlockVector[:,ii].astype( nm.float64 ) )
+
+ lambdaHistory.append( _lambda )
+
if verbosityLevel > 10:
print 'lambda:', _lambda
## # Normalize eigenvectors!
@@ -437,9 +548,17 @@
print 'final eigenvalue:', _lambda
print 'final residual norms:', residualNorms
+ if retLambdaHistory:
+ if retResidualNormsHistory:
+ return _lambda, eigBlockVectorX, lambdaHistory, residualNormsHistory
+ else:
+ return _lambda, eigBlockVectorX, lambdaHistory
+ else:
+ if retResidualNormsHistory:
+ return _lambda, eigBlockVectorX, residualNormsHistory
+ else:
+ return _lambda, eigBlockVectorX
- return _lambda, eigBlockVectorX
-
###########################################################################
if __name__ == '__main__':
from scipy.sparse import spdiags, speye
@@ -457,7 +576,7 @@
Y = nm.eye( n, 3 )
-## X = sc.rand( n, 3 )
+# X = sc.rand( n, 3 )
xfile = {100 : 'X.txt', 1000 : 'X2.txt', 10000 : 'X3.txt'}
X = nm.fromfile( xfile[n], dtype = nm.float64, sep = ' ' )
X.shape = (n, 3)
@@ -471,6 +590,8 @@
return as2d( y )
+# precond = spdiags( ivals, 0, n, n )
+
tt = time.clock()
eigs, vecs = lobpcg( X, operatorA, operatorB, blockVectorY = Y,
operatorT = precond,
Modified: trunk/Lib/sandbox/lobpcg/setup.py
===================================================================
--- trunk/Lib/sandbox/lobpcg/setup.py 2007-05-24 08:08:11 UTC (rev 3038)
+++ trunk/Lib/sandbox/lobpcg/setup.py 2007-05-24 12:30:05 UTC (rev 3039)
@@ -7,7 +7,7 @@
from numpy.distutils.system_info import get_info
config = Configuration('lobpcg',parent_package,top_path)
-# config.add_data_dir('tests')
+ config.add_data_dir('tests')
return config
Added: trunk/Lib/sandbox/lobpcg/tests/benchmark.py
===================================================================
--- trunk/Lib/sandbox/lobpcg/tests/benchmark.py 2007-05-24 08:08:11 UTC (rev 3038)
+++ trunk/Lib/sandbox/lobpcg/tests/benchmark.py 2007-05-24 12:30:05 UTC (rev 3039)
@@ -0,0 +1,64 @@
+from scipy import *
+from scipy.sandbox import lobpcg
+from symeig import symeig
+from pylab import plot, show, legend, xlabel, ylabel
+set_printoptions(precision=3,linewidth=90)
+import time
+
+def test(n):
+ x = arange(1,n+1)
+ B = diag(1./x)
+ y = arange(n-1,0,-1)
+ z = arange(2*n-1,0,-2)
+ A = diag(z)-diag(y,-1)-diag(y,1)
+ return A,B
+
+def as2d( ar ):
+ if ar.ndim == 2:
+ return ar
+ else: # Assume 1!
+ aux = nm.array( ar, copy = False )
+ aux.shape = (ar.shape[0], 1)
+ return aux
+
+def precond(x):
+ y= linalg.cho_solve((LorU, lower),x)
+ return as2d(y)
+
+m = 10 # Blocksize
+N = array(([128,256,512,1024,2048])) # Increasing matrix size
+
+data1=[]
+data2=[]
+
+for n in N:
+ print '******', n
+ A,B = test(n) # Mikota pair
+ X = rand(n,m)
+ X = linalg.orth(X)
+
+ tt = time.clock()
+ (LorU, lower) = linalg.cho_factor(A, lower=0, overwrite_a=0)
+ eigs,vecs = lobpcg.lobpcg(X,A,B,operatorT = precond,
+ residualTolerance = 1e-4, maxIterations = 40)
+ data1.append(time.clock()-tt)
+ eigs = sort(eigs)
+ print
+ print 'Results by LOBPCG'
+ print
+ print n,eigs
+
+ tt = time.clock()
+ w,v=symeig(A,B,range=(1,m))
+ data2.append(time.clock()-tt)
+ print
+ print 'Results by symeig'
+ print
+ print n, w
+
+xlabel(r'Size $n$')
+ylabel(r'Elapsed time $t$')
+plot(N,data1,label='LOBPCG')
+plot(N,data2,label='SYMEIG')
+legend()
+show()
Added: trunk/Lib/sandbox/lobpcg/tests/test_lobpcg.py
===================================================================
--- trunk/Lib/sandbox/lobpcg/tests/test_lobpcg.py 2007-05-24 08:08:11 UTC (rev 3038)
+++ trunk/Lib/sandbox/lobpcg/tests/test_lobpcg.py 2007-05-24 12:30:05 UTC (rev 3039)
@@ -0,0 +1,47 @@
+from scipy import *
+from scipy.sandbox import lobpcg
+from symeig import symeig
+from pylab import plot, show, legend, xlabel, ylabel
+set_printoptions(precision=3,linewidth=90)
+
+def test1(n):
+ L = 1.0
+ le=L/n
+ rho = 7.85e3
+ S = 1.e-4
+ E = 2.1e11
+ mass = rho*S*le/6.
+ k = E*S/le
+ A = k*(diag(r_[2.*ones(n-1),1])-diag(ones(n-1),1)-diag(ones(n-1),-1))
+ B = mass*(diag(r_[4.*ones(n-1),2])+diag(ones(n-1),1)+diag(ones(n-1),-1))
+ return A,B
+
+def test2(n):
+ x = arange(1,n+1)
+ B = diag(1./x)
+ y = arange(n-1,0,-1)
+ z = arange(2*n-1,0,-2)
+ A = diag(z)-diag(y,-1)-diag(y,1)
+ return A,B
+
+n = 100 # Dimension
+
+A,B = test1(n) # Fixed-free elastic rod
+A,B = test2(n) # Mikota pair acts as a nice test since the eigenvalues are the squares of the integers n, n=1,2,...
+
+m = 20
+V = rand(n,m)
+X = linalg.orth(V)
+
+eigs,vecs = lobpcg.lobpcg(X,A,B)
+eigs = sort(eigs)
+
+w,v=symeig(A,B)
+
+
+plot(arange(0,len(w[:m])),w[:m],'bx',label='Results by symeig')
+plot(arange(0,len(eigs)),eigs,'r+',label='Results by lobpcg')
+legend()
+xlabel(r'Eigenvalue $i$')
+ylabel(r'$\lambda_i$')
+show()
More information about the Scipy-svn
mailing list