[Scipy-svn] r3618 - in trunk/scipy/sandbox/multigrid: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Fri Dec 7 16:38:02 EST 2007
Author: wnbell
Date: 2007-12-07 15:37:58 -0600 (Fri, 07 Dec 2007)
New Revision: 3618
Modified:
trunk/scipy/sandbox/multigrid/tests/test_utils.py
trunk/scipy/sandbox/multigrid/utils.py
Log:
improved approximate_spectral_radius
added unittests
Modified: trunk/scipy/sandbox/multigrid/tests/test_utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_utils.py 2007-12-05 19:12:44 UTC (rev 3617)
+++ trunk/scipy/sandbox/multigrid/tests/test_utils.py 2007-12-07 21:37:58 UTC (rev 3618)
@@ -3,17 +3,41 @@
import numpy
import scipy
from numpy import matrix,array,diag,zeros,sqrt
+from scipy import rand
from scipy.sparse import csr_matrix
+from scipy.linalg import norm
-
set_package_path()
-from scipy.sandbox.multigrid.utils import infinity_norm, diag_sparse, \
+from scipy.sandbox.multigrid.utils import approximate_spectral_radius, \
+ infinity_norm, diag_sparse, \
symmetric_rescaling, \
expand_into_blocks
restore_path()
class TestUtils(NumpyTestCase):
+ def check_approximate_spectral_radius(self):
+ cases = []
+
+ cases.append( matrix([[-4]]) )
+ cases.append( array([[-4]]) )
+
+ cases.append( array([[2,0],[0,1]]) )
+ cases.append( array([[-2,0],[0,1]]) )
+
+ cases.append( array([[100,0,0],[0,101,0],[0,0,99]]) )
+
+ for i in range(1,5):
+ cases.append( rand(i,i) )
+
+ # method should be almost exact for small matrices
+ for A in cases:
+ Asp = csr_matrix(A)
+ assert_almost_equal( approximate_spectral_radius(A), norm(A,2) )
+ assert_almost_equal( approximate_spectral_radius(Asp), norm(A,2) )
+
+ #TODO test larger matrices
+
def check_infinity_norm(self):
A = matrix([[-4]])
assert_equal(infinity_norm(csr_matrix(A)),4)
Modified: trunk/scipy/sandbox/multigrid/utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/utils.py 2007-12-05 19:12:44 UTC (rev 3617)
+++ trunk/scipy/sandbox/multigrid/utils.py 2007-12-07 21:37:58 UTC (rev 3618)
@@ -5,48 +5,105 @@
import scipy
from scipy import ravel,arange,concatenate,tile,asarray,sqrt,diff, \
rand,zeros,empty,asmatrix,dot
-from scipy.linalg import norm,eig
+from scipy.linalg import norm,eigvals
from scipy.sparse import isspmatrix,isspmatrix_csr,isspmatrix_csc, \
csr_matrix,csc_matrix,extract_diagonal, \
coo_matrix
-def approximate_spectral_radius(A,tol=0.1,maxiter=8):
+def approximate_spectral_radius(A,tol=0.1,maxiter=10,symmetric=None):
+ """approximate the spectral radius of a matrix
+
+ *Parameters*:
+ A : dense or sparse matrix
+ E.g. csr_matrix, csc_matrix, ndarray, etc.
+
+ tol : {scalar}
+ Tolerance of approximation
+
+ maxiter : {integer}
+ Maximum number of iterations to perform
+
+ symmetric : {None,boolean}
+ True - if A is symmetric
+ Lanczos iteration is used (more efficient)
+ False - if A is non-symmetric
+ Arnoldi iteration is used (less efficient)
+ None - symmetry of A unknown
+ Method chosen automatically (default)
+ *Returns*:
+ An approximation to the spectral radius of A (scalar value)
+
"""
- Approximate the spectral radius of a matrix
- """
#from scipy.sandbox.arpack import eigen
#return norm(eigen(A, k=1, ncv=10, which='LM', maxiter=maxiter, tol=tol, return_eigenvectors=False))
+ if not isspmatrix(A):
+ A = asmatrix(A) #convert dense arrays to matrix type
+
+ if A.shape[0] != A.shape[1]:
+ raise ValueError,'expected square matrix'
+
+ #TODO make method adaptive
+
numpy.random.seed(0) #make results deterministic
- #TODO profile vs V -> V.T
- #TODO make algorithm adaptive
+ v0 = rand(A.shape[1],1)
+ v0 /= norm(v0)
- if not isspmatrix(A):
- #convert dense arrays to matrix type
- A = asmatrix(A)
+ H = zeros((maxiter+1,maxiter))
+ V = [v0]
- v0 = rand(A.shape[0])
- H = zeros((maxiter+1,maxiter))
- V = zeros((A.shape[0],maxiter+1))
-
- V[:,0]= v0/norm(v0)
+ #save past estimates
+ #estimates = []
+
for j in range(maxiter):
- w = A * V[:,j]
- for i in range(j+1):
- H[i,j] = dot(w,V[:,i])
- w -= H[i,j]*V[:,i]
- H[j+1,j] = norm(w)
- if (H[j+1,j] < 1e-12): break
- V[:,j+1] = (1.0/H[j+1,j]) * w
- # end
- m=j
- Heig,tmp = eig(H[0:m,0:m])
- return max( [norm(x) for x in Heig] )
+ w = A * V[-1]
+
+ if symmetric:
+ if j >= 1:
+ H[j-1,j] = beta
+ w -= beta * V[-2]
+ alpha = dot(ravel(w),ravel(V[-1]))
+ H[j,j] = alpha
+ w -= alpha * V[-1]
+
+ beta = norm(w)
+ if (H[j+1,j] < 1e-10): break
+
+ w /= beta
+ H[j+1,j] = beta
+ V.append(w)
+ V = V[-2:] #retain only last two vectors
+ else:
+ #orthogonalize against Vs
+ for i,v in enumerate(V):
+ H[i,j] = dot(ravel(w),ravel(v))
+ w -= H[i,j]*v
+ H[j+1,j] = norm(w)
+ if (H[j+1,j] < 1e-10): break
+
+ w /= H[j+1,j]
+ V.append(w)
+
+ # if upper 2x2 block of Hessenberg matrix H is almost symmetric,
+ # and the user has not explicitly specified symmetric=False,
+ # then switch to symmetric Lanczos algorithm
+ if symmetric is not False and j == 1:
+ if abs(H[1,0] - H[0,1]) < 1e-12:
+ symmetric = True
+ V = V[1:]
+ H[1,0] = H[0,1]
+ beta = H[2,1]
+
+ #estimates.append( max( [norm(x) for x in eigvals(H[:j+1,:j+1])] ) )
+ return norm(H[:j+1,:j+1],2)
+
+
+
def infinity_norm(A):
"""
Infinity norm of a sparse matrix (maximum absolute row sum). This serves
More information about the Scipy-svn
mailing list