scikit-learn · glemaitre · Sep 27, 2021 · Sep 27, 2021
diff --git a/sklearn/externals/_lobpcg.py b/sklearn/externals/_lobpcg.py
@@ -1,5 +1,5 @@
 """
-scikit-learn copy of scipy/sparse/linalg/eigen/lobpcg/lobpcg.py v1.3.0
+scikit-learn copy of scipy/sparse/linalg/eigen/lobpcg/lobpcg.py v1.7.1
 to be deleted after scipy 1.3.0 becomes a dependency in scikit-lean
 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 Locally Optimal Block Preconditioned Conjugate Gradient Method (LOBPCG).
@@ -10,41 +10,28 @@
        Toward the Optimal Preconditioned Eigensolver: Locally Optimal
        Block Preconditioned Conjugate Gradient Method.
        SIAM Journal on Scientific Computing 23, no. 2,
-       pp. 517-541. http://dx.doi.org/10.1137/S1064827500366124
+       pp. 517-541. :doi:`10.1137/S1064827500366124`
 
 .. [2] A. V. Knyazev, I. Lashuk, M. E. Argentati, and E. Ovchinnikov (2007),
        Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX)
-       in hypre and PETSc.  https://arxiv.org/abs/0705.2626
+       in hypre and PETSc.  :arxiv:`0705.2626`
 
 .. [3] A. V. Knyazev's C and MATLAB implementations:
-       https://bitbucket.org/joseroman/blopex
+       https://github.com/lobpcg/blopex
 """
 
-from __future__ import division, print_function, absolute_import
 import numpy as np
-from scipy.linalg import (inv, eigh, cho_factor, cho_solve, cholesky, orth,
+from scipy.linalg import (inv, eigh, cho_factor, cho_solve, cholesky,
                           LinAlgError)
 from scipy.sparse.linalg import aslinearoperator
+from numpy import block as bmat
 
 __all__ = ['lobpcg']
 
 
-def bmat(*args, **kwargs):
-    import warnings
-    with warnings.catch_warnings(record=True):
-        warnings.filterwarnings(
-            'ignore', '.*the matrix subclass is not the recommended way.*')
-        return np.bmat(*args, **kwargs)
-
-
-def _save(ar, fileName):
-    # Used only when verbosity level > 10.
-    np.savetxt(fileName, ar)
-
-
 def _report_nonhermitian(M, name):
     """
-    Report if `M` is not a hermitian matrix given its type.
+    Report if `M` is not a Hermitian matrix given its type.
     """
     from scipy.linalg import norm
 
@@ -118,7 +105,7 @@ def _b_orthonormalize(B, blockVectorV, blockVectorBV=None, retInvR=False):
         else:
             blockVectorBV = None
     except LinAlgError:
-        # raise ValueError('Cholesky has failed')
+        #raise ValueError('Cholesky has failed')
         blockVectorV = None
         blockVectorBV = None
         VBV = None
@@ -142,7 +129,7 @@ def _get_indx(_lambda, num, largest):
 
 def lobpcg(A, X,
            B=None, M=None, Y=None,
-           tol=None, maxiter=20,
+           tol=None, maxiter=None,
            largest=True, verbosityLevel=0,
            retLambdaHistory=False, retResidualNormsHistory=False):
     """Locally Optimal Block Preconditioned Conjugate Gradient Method (LOBPCG)
@@ -172,7 +159,7 @@ def lobpcg(A, X,
         Solver tolerance (stopping criterion).
         The default is ``tol=n*sqrt(eps)``.
     maxiter : int, optional
-        Maximum number of iterations.  The default is ``maxiter=min(n, 20)``.
+        Maximum number of iterations.  The default is ``maxiter = 20``.
     largest : bool, optional
         When True, solve for the largest eigenvalues, otherwise the smallest.
     verbosityLevel : int, optional
@@ -213,8 +200,7 @@ def lobpcg(A, X,
     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 LOBPCG with ``m=1``
     and ``n=10``, it works though ``n`` is small. The method is intended
-    for extremely large ``n / m``, see e.g., reference [28] in
-    https://arxiv.org/abs/0705.2626
+    for extremely large ``n / m`` [4]_.
 
     The convergence speed depends basically on two factors:
 
@@ -234,15 +220,21 @@ def lobpcg(A, X,
            Toward the Optimal Preconditioned Eigensolver: Locally Optimal
            Block Preconditioned Conjugate Gradient Method.
            SIAM Journal on Scientific Computing 23, no. 2,
-           pp. 517-541. http://dx.doi.org/10.1137/S1064827500366124
+           pp. 517-541. :doi:`10.1137/S1064827500366124`
 
     .. [2] A. V. Knyazev, I. Lashuk, M. E. Argentati, and E. Ovchinnikov
            (2007), Block Locally Optimal Preconditioned Eigenvalue Xolvers
-           (BLOPEX) in hypre and PETSc. https://arxiv.org/abs/0705.2626
+           (BLOPEX) in hypre and PETSc. :arxiv:`0705.2626`
 
     .. [3] A. V. Knyazev's C and MATLAB implementations:
            https://bitbucket.org/joseroman/blopex
 
+    .. [4] S. Yamada, T. Imamura, T. Kano, and M. Machida (2006),
+           High-performance computing for exact numerical approaches to
+           quantum many-body problems on the earth simulator. In Proceedings
+           of the 2006 ACM/IEEE Conference on Supercomputing.
+           :doi:`10.1145/1188455.1188504`
+
     Examples
     --------
 
@@ -270,7 +262,8 @@ def lobpcg(A, X,
     Initial guess for eigenvectors, should have linearly independent
     columns. Column dimension = number of requested eigenvalues.
 
-    >>> X = np.random.rand(n, 3)
+    >>> rng = np.random.default_rng()
+    >>> X = rng.random((n, 3))
 
     Preconditioner in the inverse of A in this example:
 
@@ -302,7 +295,8 @@ def lobpcg(A, X,
     blockVectorX = X
     blockVectorY = Y
     residualTolerance = tol
-    maxIterations = maxiter
+    if maxiter is None:
+        maxiter = 20
 
     if blockVectorY is not None:
         sizeY = blockVectorY.shape[1]
@@ -429,7 +423,7 @@ def lobpcg(A, X,
     iterationNumber = -1
     restart = True
     explicitGramFlag = False
-    while iterationNumber < maxIterations:
+    while iterationNumber < maxiter:
         iterationNumber += 1
         if verbosityLevel > 0:
             print('iteration %d' % iterationNumber)
@@ -487,15 +481,13 @@ def lobpcg(A, X,
         ##
         # B-orthogonalize the preconditioned residuals to X.
         if B is not None:
-            activeBlockVectorR = activeBlockVectorR - \
-                np.matmul(blockVectorX,
-                          np.matmul(blockVectorBX.T.conj(),
-                                    activeBlockVectorR))
+            activeBlockVectorR = activeBlockVectorR - np.matmul(blockVectorX,
+                                 np.matmul(blockVectorBX.T.conj(),
+                                 activeBlockVectorR))
         else:
-            activeBlockVectorR = activeBlockVectorR - \
-                np.matmul(blockVectorX,
-                          np.matmul(blockVectorX.T.conj(),
-                                    activeBlockVectorR))
+            activeBlockVectorR = activeBlockVectorR - np.matmul(blockVectorX,
+                                 np.matmul(blockVectorX.T.conj(),
+                                 activeBlockVectorR))
 
         ##
         # B-orthonormalize the preconditioned residuals.