From fa0ecc74d78375b92c2779091bb0de0e9e9b6d9c Mon Sep 17 00:00:00 2001 From: Venkataramana Ganesh Date: Wed, 16 Sep 2020 00:55:08 -0700 Subject: [PATCH] refactoring to pep8 format #1 --- python/cuml/sparse/linalg/_cp_linalg_utils.py | 11 +- python/cuml/sparse/linalg/blopex_lobpcg.py | 97 ++++++------ python/cuml/sparse/linalg/lobpcg.py | 50 ++++-- python/cuml/sparse/linalg/robust_lobpcg.py | 80 ++++++---- python/cuml/test/test_lobpcg.py | 148 ++++++++++-------- 5 files changed, 218 insertions(+), 168 deletions(-) diff --git a/python/cuml/sparse/linalg/_cp_linalg_utils.py b/python/cuml/sparse/linalg/_cp_linalg_utils.py index 8a9285db6a..0c12b4c03f 100644 --- a/python/cuml/sparse/linalg/_cp_linalg_utils.py +++ b/python/cuml/sparse/linalg/_cp_linalg_utils.py @@ -19,7 +19,8 @@ import cupyx.scipy.sparse """ -A Set of utility methods that is called by the `robust_lobpcg.py` implementation. +A Set of utility methods that is called by the +`robust_lobpcg.py` implementation. """ @@ -51,10 +52,12 @@ def matmul(A, B): return A.multiply(B) return cp.matmul(A, B) + def transpose(A): ndim = len(A.shape) return A.transpose(ndim - 1, ndim - 2) + def qform(A, S): return matmul(transpose(S), matmul(A, S)) @@ -65,15 +68,15 @@ def basis(A): def symeig(A, largest=False): - stride=1 + stride = 1 if largest is None: largest = False assert len(A.shape) == 2, "only accepts 2D matrix!" E, Z = cp.linalg.eigh(A) idx = cp.argsort(E) if largest: - stride=-1 + stride = -1 idx = idx[::stride] E = E[idx] Z = Z[:, idx] - return E, Z \ No newline at end of file + return E, Z diff --git a/python/cuml/sparse/linalg/blopex_lobpcg.py b/python/cuml/sparse/linalg/blopex_lobpcg.py index 697a1eb9e3..5670ee7060 100644 --- a/python/cuml/sparse/linalg/blopex_lobpcg.py +++ b/python/cuml/sparse/linalg/blopex_lobpcg.py @@ -31,9 +31,9 @@ """ import cupy as cp -from cupy.linalg import (solve,cholesky, inv, eigh) +from cupy.linalg import (solve, cholesky, inv, eigh) from cupyx.scipy.sparse import (spmatrix, issparse) -import warnings +# import warnings from warnings import warn __all__ = ['lobpcg'] @@ -68,6 +68,7 @@ def _as2d(ar): aux.shape = (ar.shape[0], 1) return aux + def _applyConstraints(blockVectorV, YBY, blockVectorBY, blockVectorY): """Changes blockVectorV in place.""" YBV = cp.dot(blockVectorBY.T.conj(), blockVectorV) @@ -82,7 +83,7 @@ def _b_orthonormalize(B, blockVectorV, blockVectorBV=None, retInvR=False): blockVectorV = blockVectorV / normalization if blockVectorBV is None: if B is not None: - blockVectorBV = cp.matmul(B,blockVectorV) + blockVectorBV = cp.matmul(B, blockVectorV) else: blockVectorBV = blockVectorV # Shared data!!! else: @@ -90,15 +91,16 @@ def _b_orthonormalize(B, blockVectorV, blockVectorBV=None, retInvR=False): VBV = cp.matmul(blockVectorV.T.conj(), blockVectorBV) try: # VBV is a Cholesky factor from now on... - VBV = cp.transpose(cholesky(VBV)) #cupy cholesky module returns lower triangular matrix!!! + # cupy cholesky module returns lower triangular matrix!!! + VBV = cp.transpose(cholesky(VBV)) VBV = inv(VBV) blockVectorV = cp.matmul(blockVectorV, VBV) if B is not None: blockVectorBV = cp.matmul(blockVectorBV, VBV) else: blockVectorBV = None - except: - #raise ValueError('Cholesky has failed') + except Exception as e: + # raise ValueError('Cholesky has failed') blockVectorV = None blockVectorBV = None VBV = None @@ -121,30 +123,35 @@ def _get_indx(_lambda, num, largest): def _genEigh(A, B): """ - Helper function for converting a generalized eigenvalue problem AX = lambdaBX to standard using cholesky. - This is because cupy does not have a functional api to solve generalized eigenvalue problem. + Helper function for converting a generalized eigenvalue problem + AX = lambdaBX to standard using cholesky. + This is because cupy does not have a functional api to solve + generalized eigenvalue problem. Factorizing B = R^TR. Let F = (R^T)^-1 A R^-1 - Equivalent Standard form: Fy = lambda(y), where our required eigvec x = R^-1 y + Equivalent Standard form: Fy = lambda(y), where our required + eigvec x = R^-1 y """ - R = cp.transpose(cholesky(B)) #we transpose lower triangular matrix to get upper triangular matrix + # we transpose lower triangular matrix to get upper triangular matrix + R = cp.transpose(cholesky(B)) RTi = inv(cp.transpose(R)) Ri = inv(R) - F = cp.matmul(RTi, cp.matmul(A,Ri)) + F = cp.matmul(RTi, cp.matmul(A, Ri)) vals, vecs = eigh(F) - eigVec = cp.matmul(Ri,vecs) + eigVec = cp.matmul(Ri, vecs) return vals, eigVec + def _bmat(list_obj): """ - Helper function to create a block matrix in cupy from a list of smaller 2D matrices - using iterated vertical and horizontal stacking + Helper function to create a block matrix in cupy from a list + of smaller 2D matrices using iterated vertical and horizontal stacking """ arr_rows = cp.array([]) - for row in list_obj: #row & list_obj are a list of cupy arrays + for row in list_obj: # row & list_obj are a list of cupy arrays arr_cols = cp.array([]) - for col in row: #col is a cupy ndarray + for col in row: # col is a cupy ndarray if col.ndim == 0: col = cp.array([col]) if(arr_cols.size == 0): @@ -183,7 +190,7 @@ def lobpcg(A, Parameters ---------- A : {cupy.ndarray} - The symmetric linear operator of the problem, Is treated as a 2D matrix. + The symmetric linear operator of the problem, Is treated as a 2D matrix Often called the "stiffness matrix". X : cupy.ndarray, float32 or float64 Initial approximation to the ``k`` eigenvectors (non-sparse). If `A` @@ -197,7 +204,7 @@ def lobpcg(A, Y : ndarray, float32 or float64, optional n-by-sizeY matrix of constraints (non-sparse), sizeY < n The iterations will be performed in the B-orthogonal complement - of the column-space of Y. Y must be full rank. Currently exposed to the user. + of the column-space of Y. Y must be full rank. tol : scalar, optional Solver tolerance (stopping criterion). The default is ``tol=n*sqrt(eps)``. @@ -223,7 +230,7 @@ def lobpcg(A, The history of residual norms, if `retResidualNormsHistory` is True. """ - #print('A: {}\nX:{}\nmaxiter:{}\n'.format(A,X,maxiter)) + blockVectorX = X blockVectorY = Y residualTolerance = tol @@ -262,10 +269,9 @@ def lobpcg(A, aux += "%d constraint\n\n" % sizeY print(aux) - if (n - sizeY) < (5 * sizeX): - warn('The problem size is small compared to the block size!' \ - ' Using the standard eigen solver instead of LOBPCG.') + warn('The problem size is small compared to the block size!' + ' Using the standard eigen solver instead of LOBPCG.') sizeX = min(sizeX, n) @@ -281,7 +287,7 @@ def lobpcg(A, vals = vals[::-1] vecs = vecs[:, ::-1] - return vals[:sizeX], vecs[:,:sizeX] + return vals[:sizeX], vecs[:, :sizeX] if (residualTolerance is None) or (residualTolerance <= 0.0): residualTolerance = cp.sqrt(1e-15) * n @@ -290,7 +296,7 @@ def lobpcg(A, if blockVectorY is not None: if B is not None: - blockVectorBY = cp.matmul(B,blockVectorY) + blockVectorBY = cp.matmul(B, blockVectorY) else: blockVectorBY = blockVectorY @@ -305,7 +311,7 @@ def lobpcg(A, ## # Compute the initial Ritz vectors: solve the eigenproblem. - blockVectorAX = cp.matmul(A,blockVectorX) + blockVectorAX = cp.matmul(A, blockVectorX) gramXAX = cp.dot(blockVectorX.T.conj(), blockVectorAX) _lambda, eigBlockVector = eigh(gramXAX) @@ -359,7 +365,7 @@ def lobpcg(A, ii = cp.where(residualNorms > residualTolerance, True, False) activeMask = activeMask & ii if verbosityLevel > 2: - print('the mask depending on var largest:',activeMask) + print('the mask depending on var largest:', activeMask) currentBlockSize = activeMask.sum() currentBlockSize = int(currentBlockSize) @@ -376,7 +382,7 @@ def lobpcg(A, print('eigenvalue:', _lambda) print('residual norms:', residualNorms) if verbosityLevel > 10: - print('Eigen Block Vector:',eigBlockVector) + print('Eigen Block Vector:', eigBlockVector) activeBlockVectorR = _as2d(blockVectorR[:, activeMask]) @@ -388,7 +394,7 @@ def lobpcg(A, if M is not None: # Apply preconditioner T to the active residuals. - activeBlockVectorR = cp.matmul(M,activeBlockVectorR) + activeBlockVectorR = cp.matmul(M, activeBlockVectorR) ## # Apply constraints to the preconditioned residuals. @@ -412,7 +418,7 @@ def lobpcg(A, aux = _b_orthonormalize(B, activeBlockVectorR) activeBlockVectorR, activeBlockVectorBR = aux - activeBlockVectorAR = cp.matmul(A,activeBlockVectorR) + activeBlockVectorAR = cp.matmul(A, activeBlockVectorR) if iterationNumber > 0: if B is not None: @@ -474,7 +480,6 @@ def lobpcg(A, assert isinstance(currentBlockSize, (int)) gramXBR = cp.zeros((sizeX, currentBlockSize), dtype=A.dtype) - if not restart: gramXAP = cp.dot(blockVectorX.T.conj(), activeBlockVectorAP) gramRAP = cp.dot(activeBlockVectorR.T.conj(), activeBlockVectorAP) @@ -489,21 +494,19 @@ def lobpcg(A, gramPBP = ident gramA = _bmat([[gramXAX, gramXAR, gramXAP], - [gramXAR.T.conj(), gramRAR, gramRAP], - [gramXAP.T.conj(), - gramRAP.T.conj(), gramPAP]]) + [gramXAR.T.conj(), gramRAR, gramRAP], + [gramXAP.T.conj(), + gramRAP.T.conj(), gramPAP]]) gramB = _bmat([[gramXBX, gramXBR, gramXBP], - [gramXBR.T.conj(), gramRBR, gramRBP], - [gramXBP.T.conj(), - gramRBP.T.conj(), gramPBP]]) - - + [gramXBR.T.conj(), gramRBR, gramRBP], + [gramXBP.T.conj(), + gramRBP.T.conj(), gramPBP]]) _handle_gramA_gramB_verbosity(gramA, gramB, verbosityLevel) try: _lambda, eigBlockVector = _genEigh(gramA, - gramB) + gramB) except: # try again after dropping the direction vectors P from RR restart = True @@ -516,15 +519,16 @@ def lobpcg(A, try: _lambda, eigBlockVector = _genEigh(gramA, - gramB) + gramB) - except: - raise ValueError('eigh has failed in lobpcg iterations') + except Exception as e: + raise ValueError('''eigh has failed in lobpcg iterations with + with exception {}\n'''.format(e)) ii = _get_indx(_lambda, sizeX, largest) if verbosityLevel > 10: - print('indices',ii) - print('lambda',_lambda) + print('indices', ii) + print('lambda', _lambda) _lambda = _lambda[ii] eigBlockVector = eigBlockVector[:, ii] @@ -534,7 +538,6 @@ def lobpcg(A, if verbosityLevel > 10: print('order-updated lambda:', _lambda) - # # Normalize eigenvectors! # aux = cp.sum( eigBlockVector.conj() * eigBlockVector, 0 ) # eigVecNorms = cp.sqrt( aux ) @@ -542,7 +545,7 @@ def lobpcg(A, # eigBlockVector, aux = _b_orthonormalize( B, eigBlockVector ) if verbosityLevel > 10: - print('eigenblockvector so far:',eigBlockVector) + print('eigenblockvector so far:', eigBlockVector) # Compute Ritz vectors. if B is not None: @@ -639,4 +642,4 @@ def lobpcg(A, return _lambda, blockVectorX, residualNormsHistory else: print("val: {}\n vec:{}\n".format(_lambda, blockVectorX)) - return _lambda, blockVectorX \ No newline at end of file + return _lambda, blockVectorX diff --git a/python/cuml/sparse/linalg/lobpcg.py b/python/cuml/sparse/linalg/lobpcg.py index 5124925b35..1f9fcffa0d 100755 --- a/python/cuml/sparse/linalg/lobpcg.py +++ b/python/cuml/sparse/linalg/lobpcg.py @@ -18,13 +18,14 @@ from cuml.sparse.linalg.robust_lobpcg import lobpcg as robust_lobpcg from cuml.sparse.linalg.blopex_lobpcg import lobpcg as blopex_lobpcg + def lobpcg(A, X=None, B=None, M=None, k=None, maxiter=None, - method= 'blopex', + method='blopex', tol=None, largest=True, verbosityLevel=0, @@ -34,7 +35,6 @@ def lobpcg(A, retLambdaHistory=False, retResidualNormsHistory=False ): - """ Find the k largest (or smallest) eigenvalues and the corresponding eigenvectors of a symmetric positive defined generalized @@ -171,12 +171,21 @@ def lobpcg(A, ``(lambda, V, lambda history, residual norms history)``. In the following ``n`` denotes the matrix size and ``k`` the number of required eigenvalues (smallest or largest). - The LOBPCG code internally solves eigenproblems of the size ``3k`` on every iteration by calling the "standard" dense eigensolver, so if ``k`` 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. + The LOBPCG code internally solves eigenproblems of the size ``3k`` on every + iteration by calling the "standard" dense eigensolver, so if ``k`` 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. In the BLOPEX implementation, (i.e, method="blopex"), - if one calls the LOBPCG algorithm for ``5k > n``, it will most likely break internally, so the code tries to call the standard function instead. - It is not that ``n`` should be large for the LOBPCG to work, but rather the ratio ``n / k`` should be large. It you call LOBPCG with ``k=1`` and ``n=10``, it works though ``n`` is small. The method is intended for extremely large ``n / k`` + if one calls the LOBPCG algorithm for ``5k > n``, + it will most likely break internally, so the code tries to call the + standard function instead. + It is not that ``n`` should be large for the LOBPCG to work, but rather the + ratio ``n / k`` should be large. If LOBPCG with ``k=1`` and ``n=10``, + it works though ``n`` is small. + The method is intended for extremely large ``n / k`` - In the other 2 implementations, if ``3k > n``, it Raises an Error for the same reason as stated above. + In the other 2 implementations, if ``3k > n``, it Raises an Error for + the same reason as stated above. The convergence speed depends basically on two factors: 1. How well relatively separated the seeking eigenvalues @@ -188,15 +197,26 @@ def lobpcg(A, """ if method == 'ortho': - return robust_lobpcg(A, k= k, B=B, X=X, niter=maxiter, iK=M, tol=tol, largest=largest, method='ortho',\ - verbosityLevel = verbosityLevel, ortho_iparams = ortho_iparams, ortho_fparams = ortho_fparams,\ - ortho_bparams= ortho_bparams, retLambdaHistory = retLambdaHistory, retResidualNormsHistory= retResidualNormsHistory) + return robust_lobpcg(A, k=k, B=B, X=X, niter=maxiter, + iK=M, tol=tol, largest=largest, method='ortho', + verbosityLevel=verbosityLevel, + ortho_iparams=ortho_iparams, + ortho_fparams=ortho_fparams, + ortho_bparams=ortho_bparams, + retLambdaHistory=retLambdaHistory, + retResidualNormsHistory=retResidualNormsHistory) elif method == 'basic': - return robust_lobpcg(A, k= k, B=B, X=X, niter=maxiter, iK=M, tol=tol, largest=largest, method='basic',\ - verbosityLevel = verbosityLevel, ortho_iparams = ortho_iparams, ortho_fparams = ortho_fparams,\ - ortho_bparams= ortho_bparams, retLambdaHistory = retLambdaHistory, retResidualNormsHistory= retResidualNormsHistory) + return robust_lobpcg(A, k=k, B=B, X=X, niter=maxiter, + iK=M, tol=tol, largest=largest, method='basic', + verbosityLevel=verbosityLevel, + ortho_iparams=ortho_iparams, + ortho_fparams=ortho_fparams, + ortho_bparams=ortho_bparams, + retLambdaHistory=retLambdaHistory, + retResidualNormsHistory=retResidualNormsHistory) elif method == 'blopex': X = cp.random.randn(A, k, dtype=A.dtype) if X is None else X - return blopex_lobpcg(A, X, B=B, M=M, tol=tol, maxiter=maxiter, largest=largest, verbosityLevel=verbosityLevel,\ - retLambdaHistory=retLambdaHistory, retResidualNormsHistory=retResidualNormsHistory) - + return blopex_lobpcg(A, X, B=B, M=M, tol=tol, maxiter=maxiter, + largest=largest, verbosityLevel=verbosityLevel, + retLambdaHistory=retLambdaHistory, + retResidualNormsHistory=retResidualNormsHistory) diff --git a/python/cuml/sparse/linalg/robust_lobpcg.py b/python/cuml/sparse/linalg/robust_lobpcg.py index 713e232c42..e62a061618 100644 --- a/python/cuml/sparse/linalg/robust_lobpcg.py +++ b/python/cuml/sparse/linalg/robust_lobpcg.py @@ -21,8 +21,9 @@ __all__ = ['lobpcg'] + def isNaN(num): - return num!=num + return num != num def lobpcg(A, @@ -38,7 +39,7 @@ def lobpcg(A, ortho_iparams=None, ortho_fparams=None, ortho_bparams=None, - verbosityLevel = 0, + verbosityLevel=0, retLambdaHistory=False, retResidualNormsHistory=False ): @@ -134,7 +135,7 @@ def lobpcg(A, # A and B must have the same shapes: assert A.shape == B.shape, (A.shape, B.shape) - #dtype = _utils.get_floating_dtype(A) + # dtype = _utils.get_floating_dtype(A) dtype = cp.float32 if A.dtype not in (cp.float32, cp.float64) else A.dtype if tol is None: feps = 2.23e-16 if dtype == cp.float64 else 1.2e-07 @@ -146,7 +147,7 @@ def lobpcg(A, if (m < 3 * n): raise ValueError( - 'LPBPCG algorithm is not applicable when the number of A rows (={})' + 'LPBPCG is not applicable when the number of A rows (={})' ' is smaller than 3 x the number of requested eigenpairs (={})' .format(m, n)) @@ -195,9 +196,11 @@ def lobpcg(A, A_ = bA[i] B_ = bB[i] if bB is not None else None X_ = cp.random.randn(m, n, dtype=dtype) if bX is None else bX[i] - assert len(X_.shape) == 2 and X_.shape == (m, n), (X_.shape, (m, n)) + assert len(X_.shape) == 2 and X_.shape == ( + m, n), (X_.shape, (m, n)) iparams['batch_index'] = i - worker = LOBPCG(A_, B_, X_, iK, iparams, fparams, bparams, method, verbosityLevel) + worker = LOBPCG(A_, B_, X_, iK, iparams, fparams, + bparams, method, verbosityLevel) worker.run() bE[i] = worker.E[:k] bXret[i] = worker.X[:, :k] @@ -206,34 +209,39 @@ def lobpcg(A, if retLambdaHistory: if retResidualNormsHistory: - return bE.reshape(A.shape[:-2] + (k,)), bXret.reshape(A.shape[:-2] + (m, k)), bLambdaHistory, bResidualNormsHistory + return bE.reshape(A.shape[:-2] + (k,)), + bXret.reshape(A.shape[:-2] + (m, k)), + bLambdaHistory, bResidualNormsHistory else: - return bE.reshape(A.shape[:-2] + (k,)), bXret.reshape(A.shape[:-2] + (m, k)), bLambdaHistory + return bE.reshape(A.shape[:-2] + (k,)), + bXret.reshape(A.shape[:-2] + (m, k)), bLambdaHistory else: if retResidualNormsHistory: - return bE.reshape(A.shape[:-2] + (k,)), bXret.reshape(A.shape[:-2] + (m, k)), bResidualNormsHistory + return bE.reshape(A.shape[:-2] + (k,)), + bXret.reshape(A.shape[:-2] + (m, k)), bResidualNormsHistory else: - return bE.reshape(A.shape[:-2] + (k,)), bXret.reshape(A.shape[:-2] + (m, k)) + return bE.reshape(A.shape[:-2] + (k,)), + bXret.reshape(A.shape[:-2] + (m, k)) X = cp.random.randn(m, n, dtype=dtype) if X is None else X assert len(X.shape) == 2 and X.shape == (m, n), (X.shape, (m, n)) - worker = LOBPCG(A, B, X, iK, iparams, fparams, bparams, method,verbosityLevel) + worker = LOBPCG(A, B, X, iK, iparams, fparams, + bparams, method, verbosityLevel) worker.run() if retLambdaHistory: if retResidualNormsHistory: - return worker.E[:k], worker.X[:,:k], worker.lambdaHistory, worker.residualNormsHistory + return worker.E[:k], worker.X[:, :k], + worker.lambdaHistory, worker.residualNormsHistory else: - return worker.E[:k], worker.X[:,:k], worker.lambdaHistory + return worker.E[:k], worker.X[:, :k], worker.lambdaHistory else: if retResidualNormsHistory: - return worker.E[:k], worker.X[:,:k], worker.residualNormsHistory + return worker.E[:k], worker.X[:, :k], worker.residualNormsHistory else: - return worker.E[:k], worker.X[:,:k] - - + return worker.E[:k], worker.X[:, :k] class LOBPCG(object): @@ -269,10 +277,10 @@ def __init__(self, self.E = cp.zeros((n, ), dtype=X.dtype) self.R = cp.zeros((m, n), dtype=X.dtype) self.S = cp.zeros((m, 3 * n), dtype=X.dtype) - self.tvars = {} # type: Dict[str, ndarray] - self.ivars = {'istep': 0} # type: Dict[str, int] - self.fvars = {'_': 0.0} # type: Dict[str, float] - self.bvars = {'_': False} # type: Dict[str, bool] + self.tvars = {} + self.ivars = {'istep': 0} + self.fvars = {'_': 0.0} + self.bvars = {'_': False} self.lambdaHistory = [] self.residualNormsHistory = [] @@ -301,8 +309,10 @@ def update(self): if self.ivars['istep'] == 0: X_norm = float(cp.linalg.norm(self.X)) iX_norm = X_norm ** -1 - A_norm = float(cp.linalg.norm(_utils.matmul(self.A, self.X))) * iX_norm - B_norm = float(cp.linalg.norm(_utils.matmul(self.B, self.X))) * iX_norm + A_norm = float(cp.linalg.norm( + _utils.matmul(self.A, self.X))) * iX_norm + B_norm = float(cp.linalg.norm( + _utils.matmul(self.B, self.X))) * iX_norm self.fvars['X_norm'] = X_norm self.fvars['A_norm'] = A_norm self.fvars['B_norm'] = B_norm @@ -334,13 +344,15 @@ def update_converged_count(self): convergence criterion, see discussion in Sec 4.3 of [DuerschEtal2018]. Users may redefine this method for custom convergence criteria. """ - # (...) -> int + prev_count = self.ivars['converged_count'] tol = self.fparams['tol'] A_norm = self.fvars['A_norm'] B_norm = self.fvars['B_norm'] E, X, R = self.E, self.X, self.R - rerr = cp.linalg.norm(R, 2, (0, )) * (cp.linalg.norm(X, 2, (0, )) * (A_norm + E[:X.shape[-1]] * B_norm)) ** -1 + rerr = cp.linalg.norm(R, 2, (0, ))\ + * (cp.linalg.norm(X, 2, (0, )) + * (A_norm + E[:X.shape[-1]] * B_norm)) ** -1 converged = rerr < tol count = 0 for b in converged: @@ -377,8 +389,6 @@ def run(self): if(self.verbosityLevel > 0): print(self.__str__()) - - def _update_basic(self): """ Update or initialize iteration variables when `method == "basic"`. @@ -453,7 +463,10 @@ def _update_ortho(self): # Update E, X, P self.X[:, nc:] = mm(S_, Z[:, :n - nc]) self.E[nc:] = E_[:n - nc] - P = mm(S_, mm(Z[:, n - nc:], _utils.basis(_utils.transpose(Z[:n - nc, n - nc:])))) + P = mm( + S_, + mm(Z[:, n - nc:], + _utils.basis(_utils.transpose(Z[:n - nc, n - nc:])))) np = P.shape[-1] # check convergence @@ -507,12 +520,12 @@ def _get_rayleigh_ritz_transform(self, S): :math:`(n, n)`. """ B = self.B - mm = cp.matmul SBS = _utils.qform(B, S) d_row = SBS.diagonal(0, -2, -1) ** -0.5 d_col = d_row.reshape(d_row.shape[0], 1) - ch_inp= (SBS * d_row) * d_col - R = cp.linalg.cholesky((SBS * d_row) * d_col).transpose()#cupy linalg cholesky returns the lower triangular matrix. we transpose to get the upper one. + # cupy linalg cholesky returns the lower triangular matrix. + # we transpose to get the upper one. + R = cp.linalg.cholesky((SBS * d_row) * d_col).transpose() Rinv = cp.linalg.inv(R) return Rinv * d_col @@ -521,7 +534,7 @@ def _get_svqb(self, drop, # bool tau # float ): - # type: (ndarray, bool, float) -> ndarray + """Return B-orthonormal U. .. note:: When `drop` is `False` then `svqb` is based on the Algorithm 4 from [DuerschPhD2015] that is a slight @@ -624,7 +637,6 @@ def _get_ortho(self, U, V): BU = mm_B(self.B, U) VBU = mm(_utils.transpose(V), BU) i = j = 0 - stats = '' for i in range(i_max): U = U - mm(V, VBU) drop = False @@ -646,7 +658,7 @@ def _get_ortho(self, U, V): U_norm = cp.linalg.norm(U) BU_norm = cp.linalg.norm(BU) R = UBU - cp.eye(UBU.shape[-1], - dtype=UBU.dtype) + dtype=UBU.dtype) R_norm = cp.linalg.norm(R) rerr = float(R_norm) * float(BU_norm * U_norm) ** -1 vkey = 'ortho_UBUmI_rerr[{}, {}]'.format(i, j) diff --git a/python/cuml/test/test_lobpcg.py b/python/cuml/test/test_lobpcg.py index 127a8b9ec5..5a37c7cc5e 100644 --- a/python/cuml/test/test_lobpcg.py +++ b/python/cuml/test/test_lobpcg.py @@ -15,19 +15,16 @@ # import cupy as cp -from cupy import prof import pytest -import cuml.sparse.linalg.lobpcg from cuml.sparse.linalg.lobpcg import lobpcg as cp_lobpcg -#from cuml.common.input_utils import sparse_scipy_to_cp -from numpy.testing import assert_allclose +# from numpy.testing import assert_allclose -import sys +# import sys -import time as t +# import time as t import math @@ -35,52 +32,51 @@ import warnings -import scipy +# import scipy from scipy.sparse.linalg import lobpcg as scipy_lobpcg -@pytest.mark.parametrize("dtype", [\ - cp.float32,\ - cp.float64\ - ]) -@pytest.mark.parametrize("k,n", [\ - (1,10), (2,10), \ - (1, 100), (2, 100), (5, 100), \ - (2,1000), (10,1000), \ - #(2,10000)\ - ]) - -#TODO: debug and run tests for method ortho and basic, currently they aren't that accurate -@pytest.mark.parametrize("method", [\ - #"ortho",\ - #"basic",\ - 'blopex'\ - ]) -@pytest.mark.parametrize("largest", [True, \ - False\ - ]) -@pytest.mark.parametrize("maxiter", [\ - 20,\ - #200,\ #Not practically required - #2000\ - ]) -@pytest.mark.parametrize("isB", [True, False]) - - +@pytest.mark.parametrize("dtype", [ + cp.float32, + cp.float64 +]) +@pytest.mark.parametrize("k,n", [ + (1, 10), (2, 10), + (1, 100), (2, 100), (5, 100), + (2, 1000), (10, 1000), + # (2,10000) +]) +# TODO: debug and run tests for method ortho and basic +# currently they aren't that accurate +@pytest.mark.parametrize("method", [ + # "ortho", + # "basic", + 'blopex' +]) +@pytest.mark.parametrize("largest", [True, + False + ]) +@pytest.mark.parametrize("maxiter", [ + 20, + # 200, #Not practically required + # 2000 +]) +@pytest.mark.parametrize("isB", [True, False]) def test_lobpcg(dtype, k, n, method, largest, maxiter, isB): - #settings default value of Error flags + # settings default value of Error flags cp_error_flag = False scipy_error_flag = False - A = np.random.rand(n,n) - A = np.matmul(A, A.transpose()) #making A symmetric positive definite + A = np.random.rand(n, n) + A = np.matmul(A, A.transpose()) # making A symmetric positive definite B = np.eye(n) - X = np.random.randn(n,k) + X = np.random.randn(n, k) while(isB): - B = np.random.randn(n,n) - B = np.matmul(B,B.transpose()) #making B Symmetric positive definite - if(not math.isclose(np.linalg.det(B),0, abs_tol=1e-04)): #making sure B is not singular + B = np.random.randn(n, n) + B = np.matmul(B, B.transpose()) # making B Symmetric positive definite + # making sure B is not singular + if(not math.isclose(np.linalg.det(B), 0, abs_tol=1e-04)): break A_gpu_array = cp.asarray(A) B_gpu_array = cp.asarray(B) @@ -88,45 +84,57 @@ def test_lobpcg(dtype, k, n, method, largest, maxiter, isB): """running cupy implementation: """ - try:#(please work :3) + try: # (please work :3) with cp.prof.time_range(message='start', color_id=10): - cp_val, cp_vec = cp_lobpcg(A_gpu_array,X=X_gpu_array, B=B_gpu_array, \ - maxiter=maxiter, largest=largest, method=method) + cp_val, cp_vec = cp_lobpcg(A_gpu_array, X=X_gpu_array, + B=B_gpu_array, maxiter=maxiter, + largest=largest, method=method) cp.cuda.Stream.null.synchronize() - - isnan_and_reduc_kernel = cp.ReductionKernel('T x', 'T y', 'x == x', 'a and b', 'a', '1', 'isnan_and_reduc') + isnan_and_reduc_kernel = cp.ReductionKernel( + 'T x', 'T y', 'x == x', 'a and b', 'a', '1', 'isnan_and_reduc') is_vec_ok = isnan_and_reduc_kernel(cp_vec) if(not is_vec_ok): cp_error_flag = True print("cupy: Nan values encountered in output vec method: {}!!") - except Exception as e: - print("{} occured in running the cupy lobpcg method for method: {} k:{} and n:{} with maxiter:{}!\n".format(e,method,k,n,maxiter)) + print('''{} occured in running the cupy lobpcg method for method: {} k: {} + and n: {} with maxiter: {}!'''.format( + e, + method, + k, + n, + maxiter)) cp_error_flag = True - """running scipy implementation: """ try: - scipy_val, scipy_vec = scipy_lobpcg(A, X, B = B, maxiter=maxiter, largest=largest) - except: - print("error occured in running the scipy lobpcg method for k:{} and n:{} with maxiter:{}!\n".format(k,n,maxiter)) + scipy_val, scipy_vec = scipy_lobpcg( + A, X, B=B, maxiter=maxiter, largest=largest) + except BaseException: + print( + '''error occured in running the scipy lobpcg method for + k:{} and n:{} with maxiter:{}!\n'''.format( + k, + n, + maxiter)) scipy_error_flag = True - if(scipy_error_flag or cp_error_flag): print("cholesky failed") - warnings.warn("cholesky failed for k:{}, n:{}, maxiter:{}".format(k,n,maxiter)) + warnings.warn( + "cholesky failed for k:{}, n:{}, maxiter:{}".format(k, n, maxiter)) return - #taking absolutes of eigen vectors as they can toggle between positive and negative - abs_cp_vec = cp.absolute(cp_vec) - abs_scipy_vec = np.absolute(scipy_vec) + # taking absolutes of eigen vectors as they can toggle between positive + # and negative + # abs_cp_vec = cp.absolute(cp_vec) + # abs_scipy_vec = np.absolute(scipy_vec) """ generalized eigen value problem: AX = LBX where, @@ -138,28 +146,32 @@ def test_lobpcg(dtype, k, n, method, largest, maxiter, isB): where, LHS = AX and RHS = LBX further, we take the absolutes of these residuals """ - cp_lhs = cp.asnumpy(cp.matmul(A_gpu_array,cp_vec)) - cp_rhs = cp.asnumpy(cp.multiply(cp_val, cp.matmul(B_gpu_array,cp_vec))) + cp_lhs = cp.asnumpy(cp.matmul(A_gpu_array, cp_vec)) + cp_rhs = cp.asnumpy(cp.multiply(cp_val, cp.matmul(B_gpu_array, cp_vec))) - scipy_lhs = np.matmul(A,scipy_vec) + scipy_lhs = np.matmul(A, scipy_vec) scipy_rhs = np.multiply(scipy_val, np.matmul(B, scipy_vec)) - #absolute of the residuals + # absolute of the residuals abs_cp_res = abs(cp_lhs - cp_rhs) abs_scipy_res = abs(scipy_lhs - scipy_rhs) - #taking a eigenvalue-weighted average residues across axis=1 - cp_weighted_res = np.average(abs_cp_res, axis=1, weights=cp.asnumpy(cp_val)) + # taking a eigenvalue-weighted average residues across axis=1 + cp_weighted_res = np.average( + abs_cp_res, axis=1, weights=cp.asnumpy(cp_val)) scipy_weighted_res = np.average(abs_scipy_res, axis=1, weights=scipy_val) - #taking norm of vector + # taking norm of vector cp_norm = np.linalg.norm(cp_weighted_res) scipy_norm = np.linalg.norm(scipy_weighted_res) """ TEST-STRATEGY-------------------- - The following is a comparison after norm-reduction (ie, (n,1) -> 1) of eigen-value-weighted-reduction (ie, (n,k) -> (n,1)) + The following is a comparison after norm-reduction + (ie, (n,1) -> 1) of eigen-value-weighted-reduction (ie, (n,k) -> (n,1)) of absolutes of resulting residuals """ - if cp_norm > scipy_norm: #most of the times, cupy cuml.lobpcg is more accurate than scipy - assert math.isclose(cp_norm, scipy_norm, abs_tol=1e-2, rel_tol=1e-2), "cp_norm:{} scipy_norm:{}".format(cp_norm, scipy_norm) \ No newline at end of file + if cp_norm > scipy_norm: + assert math.isclose(cp_norm, scipy_norm, abs_tol=1e-2, + rel_tol=1e-2), "cp_norm:{} scipy_norm:{}".format( + cp_norm, scipy_norm)