Source code for pyamg.util.linalg

''' Linear Algebra Helper Routines '''

__docformat__ = "restructuredtext en"

import numpy
import scipy
import scipy.sparse as sparse
from scipy.linalg.lapack import get_lapack_funcs
from scipy.lib.lapack import calc_lwork

__all__ = ['approximate_spectral_radius', 'infinity_norm', 'norm',
           'residual_norm', 'condest', 'cond', 'ishermitian',
           'pinv_array']


[docs]def norm(x, pnorm='2'): """ 2-norm of a vector Parameters ---------- x : array_like Vector of complex or real values pnorm : string '2' calculates the 2-norm 'inf' calculates the infinity-norm Returns ------- n : float 2-norm of a vector Notes ----- - currently 1+ order of magnitude faster than scipy.linalg.norm(x), which calls sqrt(numpy.sum(real((conjugate(x)*x)),axis=0)) resulting in an extra copy - only handles the 2-norm and infinity-norm for vectors See Also -------- scipy.linalg.norm : scipy general matrix or vector norm """ # TODO check dimensions of x # TODO speedup complex case x = numpy.ravel(x) if pnorm == '2': return numpy.sqrt(numpy.inner(x.conj(), x).real) elif pnorm == 'inf': return numpy.max(numpy.abs(x)) else: raise ValueError('Only the 2-norm and infinity-norm are supported')
[docs]def infinity_norm(A): """ Infinity norm of a matrix (maximum absolute row sum). Parameters ---------- A : csr_matrix, csc_matrix, sparse, or numpy matrix Sparse or dense matrix Returns ------- n : float Infinity norm of the matrix Notes ----- - This serves as an upper bound on spectral radius. - csr and csc avoid a deep copy - dense calls scipy.linalg.norm See Also -------- scipy.linalg.norm : dense matrix norms Examples -------- >>> import numpy >>> from scipy.sparse import spdiags >>> from pyamg.util.linalg import infinity_norm >>> n=10 >>> e = numpy.ones((n,1)).ravel() >>> data = [ -1*e, 2*e, -1*e ] >>> A = spdiags(data,[-1,0,1],n,n) >>> print infinity_norm(A) 4.0 """ if sparse.isspmatrix_csr(A) or sparse.isspmatrix_csc(A): # avoid copying index and ptr arrays abs_A = A.__class__((numpy.abs(A.data), A.indices, A.indptr), shape=A.shape) return (abs_A * numpy.ones((A.shape[1]), dtype=A.dtype)).max() elif sparse.isspmatrix(A): return (abs(A) * numpy.ones((A.shape[1]), dtype=A.dtype)).max() else: return numpy.dot(numpy.abs(A), numpy.ones((A.shape[1],), dtype=A.dtype)).max()
[docs]def residual_norm(A, x, b): """Compute ||b - A*x||""" return norm(numpy.ravel(b) - A*numpy.ravel(x))
def axpy(x, y, a=1.0): """ Quick level-1 call to BLAS y = a*x+y Parameters ---------- x : array_like nx1 real or complex vector y : array_like nx1 real or complex vector a : float real or complex scalar Returns ------- y : array_like Input variable y is rewritten Notes ----- The call to get_blas_funcs automatically determines the prefix for the blas call. """ from scipy.linalg import get_blas_funcs fn = get_blas_funcs(['axpy'], [x, y])[0] fn(x, y, a) # def approximate_spectral_radius(A, tol=0.1, maxiter=10, symmetric=False): # """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 : {boolean} # True if A is symmetric, False otherwise (default) # # Returns # ------- # An approximation to the spectral radius of A # # """ # if symmetric: # method = eigen_symmetric # else: # method = eigen # # return norm( method(A, k=1, tol=0.1, which='LM', maxiter=maxiter, # return_eigenvectors=False) ) def _approximate_eigenvalues(A, tol, maxiter, symmetric=None, initial_guess=None): """Used by approximate_spectral_radius and condest Returns [W, E, H, V, breakdown_flag], where W and E are the eigenvectors and eigenvalues of the Hessenberg matrix H, respectively, and V is the Krylov space. breakdown_flag denotes whether Lanczos/Arnoldi suffered breakdown. E is therefore the approximate eigenvalues of A. To obtain approximate eigenvectors of A, compute V*W. """ from scipy.sparse.linalg import aslinearoperator A = aslinearoperator(A) # A could be dense or sparse, or something weird # Choose tolerance for deciding if break-down has occurred t = A.dtype.char eps = numpy.finfo(numpy.float).eps feps = numpy.finfo(numpy.single).eps geps = numpy.finfo(numpy.longfloat).eps _array_precision = {'f': 0, 'd': 1, 'g': 2, 'F': 0, 'D': 1, 'G': 2} breakdown = {0: feps*1e3, 1: eps*1e6, 2: geps*1e6}[_array_precision[t]] breakdown_flag = False if A.shape[0] != A.shape[1]: raise ValueError('expected square matrix') maxiter = min(A.shape[0], maxiter) if initial_guess is None: v0 = scipy.rand(A.shape[1], 1) if A.dtype == complex: v0 = v0 + 1.0j * scipy.rand(A.shape[1], 1) else: v0 = initial_guess v0 /= norm(v0) # Important to type H based on v0, so that a real nonsymmetric matrix, can # have an imaginary initial guess for its Arnoldi Krylov space H = numpy.zeros((maxiter+1, maxiter), dtype=numpy.find_common_type([v0.dtype, A.dtype], [])) V = [v0] beta = 0.0 for j in range(maxiter): w = A * V[-1] if symmetric: if j >= 1: H[j-1, j] = beta w -= beta * V[-2] alpha = numpy.dot(numpy.conjugate(w.ravel()), V[-1].ravel()) H[j, j] = alpha w -= alpha * V[-1] # axpy(V[-1],w,-alpha) beta = norm(w) H[j+1, j] = beta if (H[j+1, j] < breakdown): breakdown_flag = True break w /= 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] = numpy.dot(numpy.conjugate(v.ravel()), w.ravel()) w = w - H[i, j]*v H[j+1, j] = norm(w) if (H[j+1, j] < breakdown): breakdown_flag = True if H[j+1, j] != 0.0: w = w/H[j+1, j] V.append(w) break w = 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: # #print "using symmetric mode" # symmetric = True # V = V[1:] # H[1,0] = H[0,1] # beta = H[2,1] # print "Approximated spectral radius in %d iterations" % (j + 1) from scipy.linalg import eig Eigs, Vects = eig(H[:j+1, :j+1], left=False, right=True) return (Vects, Eigs, H, V, breakdown_flag)
[docs]def approximate_spectral_radius(A, tol=0.01, maxiter=15, restart=5, symmetric=None, initial_guess=None, return_vector=False): """ Approximate the spectral radius of a matrix Parameters ---------- A : {dense or sparse matrix} E.g. csr_matrix, csc_matrix, ndarray, etc. tol : {scalar} Relative tolerance of approximation, i.e., the error divided by the approximate spectral radius is compared to tol. maxiter : {integer} Maximum number of iterations to perform restart : {integer} Number of restarted Arnoldi processes. For example, a value of 0 will run Arnoldi once, for maxiter iterations, and a value of 1 will restart Arnoldi once, using the maximal eigenvector from the first Arnoldi process as the initial guess. symmetric : {boolean} True - if A is symmetric Lanczos iteration is used (more efficient) False - if A is non-symmetric (default Arnoldi iteration is used (less efficient) initial_guess : {array|None} If n x 1 array, then use as initial guess for Arnoldi/Lanczos. If None, then use a random initial guess. return_vector : {boolean} True - return an approximate dominant eigenvector, in addition to the spectral radius. False - Do not return the approximate dominant eigenvector Returns ------- An approximation to the spectral radius of A, and if return_vector=True, then also return the approximate dominant eigenvector Notes ----- The spectral radius is approximated by looking at the Ritz eigenvalues. Arnoldi iteration (or Lanczos) is used to project the matrix A onto a Krylov subspace: H = Q* A Q. The eigenvalues of H (i.e. the Ritz eigenvalues) should represent the eigenvalues of A in the sense that the minimum and maximum values are usually well matched (for the symmetric case it is true since the eigenvalues are real). References ---------- .. [1] Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. van der Vorst, editors. "Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide", SIAM, Philadelphia, 2000. Examples -------- >>> from pyamg.util.linalg import approximate_spectral_radius >>> import numpy >>> from scipy.linalg import eigvals, norm >>> A = numpy.array([[1.,0.],[0.,1.]]) >>> print approximate_spectral_radius(A,maxiter=3) 1.0 >>> print max([norm(x) for x in eigvals(A)]) 1.0 """ if not hasattr(A, 'rho') or return_vector: # somehow more restart causes a nonsymmetric case to fail...look at # this what about A.dtype=int? convert somehow? # The use of the restart vector v0 requires that the full Krylov # subspace V be stored. So, set symmetric to False. symmetric = False if maxiter < 1: raise ValueError('expected maxiter > 0') if restart < 0: raise ValueError('expected restart >= 0') if A.dtype == int: raise ValueError('expected A to be float (complex or real)') if A.shape[0] != A.shape[1]: raise ValueError('expected square A') if initial_guess is None: v0 = scipy.rand(A.shape[1], 1) if A.dtype == complex: v0 = v0 + 1.0j * scipy.rand(A.shape[1], 1) else: if initial_guess.shape[0] != A.shape[0]: raise ValueError('initial_guess and A must have same shape') if (len(initial_guess.shape) > 1) and (initial_guess.shape[1] > 1): raise ValueError('initial_guess must be an (n,1) or\ (n,) vector') v0 = initial_guess.reshape(-1, 1) v0 = numpy.array(v0, dtype=A.dtype) for j in range(restart+1): [evect, ev, H, V, breakdown_flag] =\ _approximate_eigenvalues(A, tol, maxiter, symmetric, initial_guess=v0) # Calculate error in dominant eigenvector nvecs = ev.shape[0] max_index = numpy.abs(ev).argmax() error = H[nvecs, nvecs-1]*evect[-1, max_index] # error is a fast way of calculating the following line # error2 = ( A - ev[max_index]*scipy.mat( # scipy.eye(A.shape[0],A.shape[1])) )*\ # ( scipy.mat(scipy.hstack(V[:-1]))*\ # evect[:,max_index].reshape(-1,1) ) # print str(error) + " " + str(scipy.linalg.norm(e2)) if (numpy.abs(error)/numpy.abs(ev[max_index]) < tol) or\ breakdown_flag: # halt if below relative tolerance v0 = numpy.dot(numpy.hstack(V[:-1]), evect[:, max_index].reshape(-1, 1)) break else: v0 = numpy.dot(numpy.hstack(V[:-1]), evect[:, max_index].reshape(-1, 1)) # end j-loop rho = numpy.abs(ev[max_index]) if sparse.isspmatrix(A): A.rho = rho if return_vector: return (rho, v0) else: return rho else: return A.rho
[docs]def condest(A, tol=0.1, maxiter=25, symmetric=False): """Estimates the condition number of A Parameters ---------- A : {dense or sparse matrix} e.g. array, matrix, csr_matrix, ... tol : {float} Approximation tolerance, currently not used maxiter: {int} Max number of Arnoldi/Lanczos iterations symmetric : {bool} If symmetric use the far more efficient Lanczos algorithm, Else use Arnoldi Returns ------- Estimate of cond(A) with \|lambda_max\| / \|lambda_min\| through the use of Arnoldi or Lanczos iterations, depending on the symmetric flag Notes ----- The condition number measures how large of a change in the the problems solution is caused by a change in problem's input. Large condition numbers indicate that small perturbations and numerical errors are magnified greatly when solving the system. Examples -------- >>> import numpy >>> from pyamg.util.linalg import condest >>> c = condest(numpy.array([[1.,0.],[0.,2.]])) >>> print c 2.0 """ [evect, ev, H, V, breakdown_flag] =\ _approximate_eigenvalues(A, tol, maxiter, symmetric) return numpy.max([norm(x) for x in ev])/min([norm(x) for x in ev])
[docs]def cond(A): """Returns condition number of A Parameters ---------- A : {dense or sparse matrix} e.g. array, matrix, csr_matrix, ... Returns ------- 2-norm condition number through use of the SVD Use for small to moderate sized dense matrices. For large sparse matrices, use condest. Notes ----- The condition number measures how large of a change in the problems solution is caused by a change in problem's input. Large condition numbers indicate that small perturbations and numerical errors are magnified greatly when solving the system. Examples -------- >>> import numpy >>> from pyamg.util.linalg import condest >>> c = condest(numpy.array([[1.0,0.],[0.,2.0]])) >>> print c 2.0 """ if A.shape[0] != A.shape[1]: raise ValueError('expected square matrix') if sparse.isspmatrix(A): A = A.todense() # 2-Norm Condition Number from scipy.linalg import svd U, Sigma, Vh = svd(A) return numpy.max(Sigma)/min(Sigma)
[docs]def ishermitian(A, fast_check=True, tol=1e-6, verbose=False): """Returns True if A is Hermitian to within tol Parameters ---------- A : {dense or sparse matrix} e.g. array, matrix, csr_matrix, ... fast_check : {bool} If True, use the heuristic < Ax, y> = < x, Ay> for random vectors x and y to check for conjugate symmetry. If False, compute A - A.H. tol : {float} Symmetry tolerance verbose: {bool} prints max( \|A - A.H\| ) if nonhermitian and fast_check=False abs( <Ax, y> - <x, Ay> ) if nonhermitian and fast_check=False Returns ------- True if hermitian False if nonhermitian Notes ----- This function applies a simple test of conjugate symmetry Examples -------- >>> import numpy >>> from pyamg.util.linalg import ishermitian >>> ishermitian(numpy.array([[1,2],[1,1]])) False >>> from pyamg.gallery import poisson >>> ishermitian(poisson((10,10))) True """ # convert to matrix type if not sparse.isspmatrix(A): A = numpy.asmatrix(A) if fast_check: x = scipy.rand(A.shape[0], 1) y = scipy.rand(A.shape[0], 1) if A.dtype == complex: x = x + 1.0j*scipy.rand(A.shape[0], 1) y = y + 1.0j*scipy.rand(A.shape[0], 1) xAy = numpy.dot((A*x).conjugate().T, y) xAty = numpy.dot(x.conjugate().T, A*y) diff = float(numpy.abs(xAy - xAty) / numpy.sqrt(numpy.abs(xAy*xAty))) else: # compute the difference, A - A.H if sparse.isspmatrix(A): diff = numpy.ravel((A - A.H).data) else: diff = numpy.ravel(A - A.H) if numpy.max(diff.shape) == 0: diff = 0 else: diff = numpy.max(numpy.abs(diff)) if diff < tol: diff = 0 return True else: if verbose: print diff return False return diff
[docs]def pinv_array(a, cond=None): """Calculate the Moore-Penrose pseudo inverse of each block of the three dimensional array a. Parameters ---------- a : {dense array} Is of size (n, m, m) cond : {float} Used by *gelss to filter numerically zeros singular values. If None, a suitable value is chosen for you. Returns ------- Nothing, a is modified in place so that a[k] holds the pseudoinverse of that block. Notes ----- By using lapack wrappers, this can be much faster for large n, than directly calling pinv2 Examples -------- >>> import numpy >>> from pyamg.util.linalg import pinv_array >>> a = numpy.array([[[1.,2.],[1.,1.]], [[1.,1.],[3.,3.]]]) >>> ac = a.copy() >>> # each block of a is inverted in-place >>> pinv_array(a) """ n = a.shape[0] m = a.shape[1] if m == 1: # Pseudo-inverse of 1 x 1 matrices is trivial zero_entries = (a == 0.0).nonzero()[0] a[zero_entries] = 1.0 a[:] = 1.0/a a[zero_entries] = 0.0 del zero_entries else: # The block size is greater than 1 # Create necessary arrays and function pointers for calculating pinv gelss, = get_lapack_funcs(('gelss',), (numpy.ones((1,), dtype=a.dtype))) RHS = numpy.eye(m, dtype=a.dtype) lwork = calc_lwork.gelss(gelss.prefix, m, m, m)[1] # Choose tolerance for which singular values are zero in *gelss below if cond is None: t = a.dtype.char eps = numpy.finfo(numpy.float).eps feps = numpy.finfo(numpy.single).eps geps = numpy.finfo(numpy.longfloat).eps _array_precision = {'f': 0, 'd': 1, 'g': 2, 'F': 0, 'D': 1, 'G': 2} cond = {0: feps*1e3, 1: eps*1e6, 2: geps*1e6}[_array_precision[t]] # Invert each block of a for kk in xrange(n): gelssoutput = gelss(a[kk], RHS, cond=cond, lwork=lwork, overwrite_a=True, overwrite_b=False) a[kk] = gelssoutput[1]