def _projection_matrix(self, r, sigma, b, c, projection):
        fom = self.fom
        if self.version == 'V':
            V = fom.A.source.empty(reserve=r)
        else:
            W = fom.A.source.empty(reserve=r)
        for i in range(r):
            if sigma[i].imag == 0:
                sEmA = sigma[i].real * self.fom.E - self.fom.A
                if self.version == 'V':
                    Bb = fom.B.apply(b.real[i])
                    V.append(sEmA.apply_inverse(Bb))
                else:
                    CTc = fom.C.apply_adjoint(c.real[i])
                    W.append(sEmA.apply_inverse_adjoint(CTc))
            elif sigma[i].imag > 0:
                sEmA = sigma[i] * self.fom.E - self.fom.A
                if self.version == 'V':
                    Bb = fom.B.apply(b[i])
                    v = sEmA.apply_inverse(Bb)
                    V.append(v.real)
                    V.append(v.imag)
                else:
                    CTc = fom.C.apply_adjoint(c[i].conj())
                    w = sEmA.apply_inverse_adjoint(CTc)
                    W.append(w.real)
                    W.append(w.imag)

        if self.version == 'V':
            self.V = gram_schmidt(V, atol=0, rtol=0, product=None if projection == 'orth' else fom.E)
        else:
            self.V = gram_schmidt(W, atol=0, rtol=0, product=None if projection == 'orth' else fom.E)

        self._pg_reductor = LTIPGReductor(fom, self.V, self.V, projection == 'Eorth')

def extend_basis(U, basis, product=None, method='gram_schmidt', pod_modes=1, pod_orthonormalize=True, copy_U=True):
    assert method in ('trivial', 'gram_schmidt', 'pod')

    basis_length = len(basis)

    if method == 'trivial':
        remove = set()
        for i in range(len(U)):
            if np.any(almost_equal(U[i], basis)):
                remove.add(i)
        basis.append(U[[i for i in range(len(U)) if i not in remove]],
                     remove_from_other=(not copy_U))
    elif method == 'gram_schmidt':
        basis.append(U, remove_from_other=(not copy_U))
        gram_schmidt(basis, offset=basis_length, product=product, copy=False, check=False)
    elif method == 'pod':
        U_proj_err = U - basis.lincomb(U.inner(basis, product))

        basis.append(pod(U_proj_err, modes=pod_modes, product=product, orthonormalize=False)[0])

        if pod_orthonormalize:
            gram_schmidt(basis, offset=basis_length, product=product, copy=False, check=False)

    if len(basis) <= basis_length:
        raise ExtensionError

def test_gram_schmidt():
    for i in (1, 32):
        b = NumpyVectorArray(np.identity(i, dtype=np.float))
        a = gram_schmidt(b)
        assert np.all(almost_equal(b, a))
    c = NumpyVectorArray([[1.0, 0], [0., 0]])
    a = gram_schmidt(c)
    assert (a.data == np.array([[1.0, 0]])).all()

    def reduce(self, r, projection='bfsr'):
        """Reduce using SOBTfv.

        Parameters
        ----------
        r
            Order of the reduced model.
        projection
            Projection method used:

                - `'sr'`: square root method
                - `'bfsr'`: balancing-free square root method (default,
                    since it avoids scaling by singular values and
                    orthogonalizes the projection matrices, which might
                    make it more accurate than the square root method)
                - `'biorth'`: like the balancing-free square root
                    method, except it biorthogonalizes the projection
                    matrices

        Returns
        -------
        rd
            Reduced system.
        """
        assert 0 < r < self.d.n
        assert projection in ('sr', 'bfsr', 'biorth')

        # compute all necessary Gramian factors
        pcf = self.d.gramian('pcf')
        pof = self.d.gramian('pof')

        if r > min(len(pcf), len(pof)):
            raise ValueError('r needs to be smaller than the sizes of Gramian factors.')

        # find necessary SVDs
        _, sp, Vp = spla.svd(pof.inner(pcf))

        # compute projection matrices and find the reduced model
        self.V = pcf.lincomb(Vp[:r])
        if projection == 'sr':
            alpha = 1 / np.sqrt(sp[:r])
            self.V.scal(alpha)
            self.bases_are_biorthonormal = False
        elif projection == 'bfsr':
            self.V = gram_schmidt(self.V, atol=0, rtol=0)
            self.bases_are_biorthonormal = False
        elif projection == 'biorth':
            self.V = gram_schmidt(self.V, product=self.d.M, atol=0, rtol=0)
            self.bases_are_biorthonormal = True

        self.W = self.V

        self.pg_reductor = GenericPGReductor(self.d, self.W, self.V, projection == 'biorth', product=self.d.M)

        rd = self.pg_reductor.reduce()

        return rd

    def reduce(self, r=None, tol=None, projection='bfsr'):
        """Generic Balanced Truncation.

        Parameters
        ----------
        r
            Order of the reduced model if `tol` is `None`, maximum order if `tol` is specified.
        tol
            Tolerance for the error bound if `r` is `None`.
        projection
            Projection method used:

            - `'sr'`: square root method
            - `'bfsr'`: balancing-free square root method (default, since it avoids scaling by
              singular values and orthogonalizes the projection matrices, which might make it more
              accurate than the square root method)
            - `'biorth'`: like the balancing-free square root method, except it biorthogonalizes the
              projection matrices (using :func:`~pymor.algorithms.gram_schmidt.gram_schmidt_biorth`)

        Returns
        -------
        rom
            Reduced-order model.
        """
        assert r is not None or tol is not None
        assert r is None or 0 < r < self.fom.order
        assert projection in ('sr', 'bfsr', 'biorth')

        cf, of = self._gramians()
        sv, sU, sV = self._sv_U_V()

        # find reduced order if tol is specified
        if tol is not None:
            error_bounds = self.error_bounds()
            r_tol = np.argmax(error_bounds <= tol) + 1
            r = r_tol if r is None else min(r, r_tol)
        if r > min(len(cf), len(of)):
            raise ValueError('r needs to be smaller than the sizes of Gramian factors.')

        # compute projection matrices
        self.V = cf.lincomb(sV[:r])
        self.W = of.lincomb(sU[:r])
        if projection == 'sr':
            alpha = 1 / np.sqrt(sv[:r])
            self.V.scal(alpha)
            self.W.scal(alpha)
        elif projection == 'bfsr':
            self.V = gram_schmidt(self.V, atol=0, rtol=0)
            self.W = gram_schmidt(self.W, atol=0, rtol=0)
        elif projection == 'biorth':
            self.V, self.W = gram_schmidt_biorth(self.V, self.W, product=self.fom.E)

        # find reduced-order model
        self._pg_reductor = LTIPGReductor(self.fom, self.W, self.V, projection in ('sr', 'biorth'))
        rom = self._pg_reductor.reduce()
        return rom

def test_gram_schmidt_with_product(operator_with_arrays_and_products):
    _, _, U, _, p, _ = operator_with_arrays_and_products

    V = U.copy()
    onb = gram_schmidt(U, product=p, copy=True)
    assert np.all(almost_equal(U, V))
    assert np.allclose(p.apply2(onb, onb), np.eye(len(onb)))
    assert np.all(almost_equal(U, onb.lincomb(p.apply2(U, onb)), rtol=1e-13))

    onb2 = gram_schmidt(U, product=p, copy=False)
    assert np.all(almost_equal(onb, onb2))
    assert np.all(almost_equal(onb, U))

def test_gram_schmidt(vector_array):
    U = vector_array

    V = U.copy()
    onb = gram_schmidt(U, copy=True)
    assert np.all(almost_equal(U, V))
    assert np.allclose(onb.dot(onb), np.eye(len(onb)))
    assert np.all(almost_equal(U, onb.lincomb(U.dot(onb)), rtol=1e-13))

    onb2 = gram_schmidt(U, copy=False)
    assert np.all(almost_equal(onb, onb2))
    assert np.all(almost_equal(onb, U))

    def _projection_matrices(self, rom, projection):
        fom = self.fom
        self.V, self.W = solve_sylv_schur(fom.A, rom.A,
                                          E=fom.E, Er=rom.E,
                                          B=fom.B, Br=rom.B,
                                          C=fom.C, Cr=rom.C)
        if projection == 'orth':
            self.V = gram_schmidt(self.V, atol=0, rtol=0)
            self.W = gram_schmidt(self.W, atol=0, rtol=0)
        elif projection == 'biorth':
            self.V, self.W = gram_schmidt_biorth(self.V, self.W, product=fom.E)

        self._pg_reductor = LTIPGReductor(fom, self.W, self.V, projection == 'biorth')

def test_gram_schmidt_with_R(vector_array):
    U = vector_array

    V = U.copy()
    onb, R = gram_schmidt(U, return_R=True, copy=True)
    assert np.all(almost_equal(U, V))
    assert np.allclose(onb.dot(onb), np.eye(len(onb)))
    assert np.all(almost_equal(U, onb.lincomb(U.dot(onb)), rtol=1e-13))
    assert np.all(almost_equal(V, onb.lincomb(R.T)))

    onb2, R2 = gram_schmidt(U, return_R=True, copy=False)
    assert np.all(almost_equal(onb, onb2))
    assert np.all(R == R2)
    assert np.all(almost_equal(onb, U))

def projection_shifts_init(A, E, B, shift_options):
    """Find starting shift parameters for low-rank ADI iteration using
    Galerkin projection on spaces spanned by LR-ADI iterates.

    See [PK16]_, pp. 92-95.

    Parameters
    ----------
    A
        The |Operator| A from the corresponding Lyapunov equation.
    E
        The |Operator| E from the corresponding Lyapunov equation.
    B
        The |VectorArray| B from the corresponding Lyapunov equation.
    shift_options
        The shift options to use (see :func:`lyap_solver_options`).

    Returns
    -------
    shifts
        A |NumPy array| containing a set of stable shift parameters.
    """
    for i in range(shift_options['init_maxiter']):
        Q = gram_schmidt(B, atol=0, rtol=0)
        shifts = spla.eigvals(A.apply2(Q, Q), E.apply2(Q, Q))
        shifts = shifts[np.real(shifts) < 0]
        if shifts.size == 0:
            # use random subspace instead of span{B} (with same dimensions)
            if shift_options['init_seed'] is not None:
                np.random.seed(shift_options['init_seed'])
                np.random.seed(np.random.random() + i)
            B = B.space.make_array(np.random.rand(len(B), B.space.dim))
        else:
            return shifts
    raise RuntimeError('Could not generate initial shifts for low-rank ADI iteration.')

def test_project_array():
    np.random.seed(123)
    U = NumpyVectorSpace.from_numpy(np.random.random((2, 10)))
    basis = NumpyVectorSpace.from_numpy(np.random.random((3, 10)))
    U_p = project_array(U, basis, orthonormal=False)
    onb = gram_schmidt(basis)
    U_p2 = project_array(U, onb, orthonormal=True)
    assert np.all(relative_error(U_p, U_p2) < 1e-10)

def gram_schmidt_basis_extension(basis, U, product=None, copy_basis=True, copy_U=True):
    """Extend basis using Gram-Schmidt orthonormalization.

    Parameters
    ----------
    basis
        |VectorArray| containing the basis to extend.
    U
        |VectorArray| containing the new basis vectors.
    product
        The scalar product w.r.t. which to orthonormalize; if `None`, the Euclidean
        product is used.
    copy_basis
        If `copy_basis` is `False`, the old basis is extended in-place.
    copy_U
        If `copy_U` is `False`, the new basis vectors are removed from `U`.

    Returns
    -------
    new_basis
        The extended basis.
    extension_data
        Dict containing the following fields:

            :hierarchic: `True` if `new_basis` contains `basis` as its first vectors.

    Raises
    ------
    ExtensionError
        Gram-Schmidt orthonormalization fails. This is the case when no
        vector in `U` is linearly independent from the basis.
    """
    if basis is None:
        basis = U.empty(reserve=len(U))

    basis_length = len(basis)

    new_basis = basis.copy() if copy_basis else basis
    new_basis.append(U, remove_from_other=(not copy_U))
    gram_schmidt(new_basis, offset=basis_length, product=product, copy=False)

    if len(new_basis) <= basis_length:
        raise ExtensionError

    return new_basis, {'hierarchic': True}

def test_project_array_with_product():
    np.random.seed(123)
    U = NumpyVectorSpace.from_numpy(np.random.random((1, 10)))
    basis = NumpyVectorSpace.from_numpy(np.random.random((3, 10)))
    product = np.random.random((10, 10))
    product = NumpyMatrixOperator(product.T.dot(product))
    U_p = project_array(U, basis, product=product, orthonormal=False)
    onb = gram_schmidt(basis, product=product)
    U_p2 = project_array(U, onb, product=product, orthonormal=True)
    assert np.all(relative_error(U_p, U_p2, product) < 1e-10)

def projection_shifts(A, E, V, prev_shifts):
    """Find further shift parameters for low-rank ADI iteration using
    Galerkin projection on spaces spanned by LR-ADI iterates.

    See [PK16]_, pp. 92-95.

    Parameters
    ----------
    A
        The |Operator| A from the corresponding Lyapunov equation.
    E
        The |Operator| E from the corresponding Lyapunov equation.
    V
        A |VectorArray| representing the currently computed iterate.
    prev_shifts
        A |NumPy array| containing the set of all previously used shift
        parameters.

    Returns
    -------
    shifts
        A |NumPy array| containing a set of stable shift parameters.
    """
    if prev_shifts[-1].imag != 0:
        Q = gram_schmidt(cat_arrays([V.real, V.imag]), atol=0, rtol=0)
    else:
        Q = gram_schmidt(V, atol=0, rtol=0)

    Ap = A.apply2(Q, Q)
    Ep = E.apply2(Q, Q)

    shifts = spla.eigvals(Ap, Ep)
    shifts.imag[abs(shifts.imag) < np.finfo(float).eps] = 0
    shifts = shifts[np.real(shifts) < 0]
    if shifts.size == 0:
        return prev_shifts
    else:
        if np.any(shifts.imag != 0):
            shifts = shifts[np.abs(shifts).argsort()]
        else:
            shifts.sort()
        return shifts

def rrf(A, source_product=None, range_product=None, q=2, l=8, iscomplex=False):
    """Randomized range approximation of `A`.

    This is an implementation of Algorithm 4.4 in [HMT11]_.

    Given the |Operator| `A`, the return value of this method is the |VectorArray|
    `Q` whose vectors form an approximate orthonomal basis for the range of `A`.

    Parameters
    ----------
    A
        The |Operator| A.
    source_product
        Inner product |Operator| of the source of A.
    range_product
        Inner product |Operator| of the range of A.
    q
        The number of power iterations.
    l
        The block size of the normalized power iterations.
    iscomplex
        If `True`, the random vectors are chosen complex.

    Returns
    -------
    Q
        |VectorArray| which contains the basis, whose span approximates the range of A.
    """

    assert source_product is None or isinstance(source_product, OperatorInterface)
    assert range_product is None or isinstance(range_product, OperatorInterface)
    assert isinstance(A, OperatorInterface)

    R = A.source.random(l, distribution='normal')
    if iscomplex:
        R += 1j*A.source.random(l, distribution='normal')
    Q = A.apply(R)
    gram_schmidt(Q, range_product, atol=0, rtol=0, copy=False)

    for i in range(q):
        Q = A.apply_adjoint(Q)
        gram_schmidt(Q, source_product, atol=0, rtol=0, copy=False)
        Q = A.apply(Q)
        gram_schmidt(Q, range_product, atol=0, rtol=0, copy=False)

    return Q

    def reduce(self, sigma, b, c, projection='orth'):
        """Bitangential Hermite interpolation.

        Parameters
        ----------
        sigma
            Interpolation points (closed under conjugation), list of
            length `r`.
        b
            Right tangential directions, |VectorArray| of length `r`
            from `self._B_source`.
        c
            Left tangential directions, |VectorArray| of length `r` from
            `self._C_range`.
        projection
            Projection method:

                - `'orth'`: projection matrices are orthogonalized with
                    respect to the Euclidean inner product
                - `'biorth'`: projection matrices are biorthogolized
                    with respect to the E product

        Returns
        -------
        rd
            Reduced discretization.
        """
        r = len(sigma)
        assert b in self._B_source and len(b) == r
        assert c in self._C_range and len(c) == r
        assert projection in ('orth', 'biorth')

        # rescale tangential directions (to avoid overflow or underflow)
        if b.dim > 1:
            b.scal(1 / b.l2_norm())
        else:
            b = self._B_source.from_numpy(np.ones((r, 1)))
        if c.dim > 1:
            c.scal(1 / c.l2_norm())
        else:
            c = self._C_range.from_numpy(np.ones((r, 1)))

        # compute projection matrices
        self.V = self._K_source.empty(reserve=r)
        self.W = self._K_source.empty(reserve=r)
        for i in range(r):
            if sigma[i].imag == 0:
                Bb = self._B_apply(sigma[i].real, b.real[i])
                self.V.append(self._K_apply_inverse(sigma[i].real, Bb))

                CTc = self._C_apply_adjoint(sigma[i].real, c.real[i])
                self.W.append(self._K_apply_inverse_adjoint(sigma[i].real, CTc))
            elif sigma[i].imag > 0:
                Bb = self._B_apply(sigma[i], b[i])
                v = self._K_apply_inverse(sigma[i], Bb)
                self.V.append(v.real)
                self.V.append(v.imag)

                CTc = self._C_apply_adjoint(sigma[i], c[i].conj())
                w = self._K_apply_inverse_adjoint(sigma[i], CTc)
                self.W.append(w.real)
                self.W.append(w.imag)

        if projection == 'orth':
            self.V = gram_schmidt(self.V, atol=0, rtol=0)
            self.W = gram_schmidt(self.W, atol=0, rtol=0)
        elif projection == 'biorth':
            self.V, self.W = gram_schmidt_biorth(self.V, self.W, product=self._product)

        self.pg_reductor = GenericPGReductor(self.d, self.W, self.V, projection == 'biorth', product=self._product)

        rd = self.pg_reductor.reduce()
        return rd

def estimate_image(operators=(), vectors=(),
                   domain=None, extends=False, orthonormalize=True, product=None,
                   riesz_representatives=False):
    """Estimate the image of given operators for all mu.

    Let `operators` be a list of |Operators| with common source and domain, and let
    `vectors` be a list of |VectorArrays| or vector-like |Operators| in the range
    of these operators. Given a |VectorArray| `domain` of vectors in the source of the
    operators, this algorithms determines a |VectorArray| `image` of range vectors
    such that the linear span of `image` contains:

    - `op.apply(U, mu=mu)` for all operators `op` in `operators`, for all possible |Parameters|
      `mu` and for all |VectorArrays| `U` contained in the linear span of `domain`,
    - `U` for all |VectorArrays| in `vectors`,
    - `v.as_vector(mu)` for all |Operators| in `vectors` and all possible |Parameters| `mu`.

    The algorithm will try to choose `image` as small as possible. However, no optimality
    is guaranteed.

    Parameters
    ----------
    operators
        See above.
    vectors
        See above.
    domain
        See above. If `None`, an empty `domain` |VectorArray| is assumed.
    extends
        For some operators, e.g. |EmpiricalInterpolatedOperator|, as well as for all
        elements of `vectors`, `image` is estimated independently from the choice of
        `domain`.  If `extends` is `True`, such operators are ignored. (This is useful
        in case these vectors have already been obtained by earlier calls to this
        function.)
    orthonormalize
        Compute an orthonormal basis for the linear span of `image` using the
        :func:`~pymor.algorithms.gram_schmidt.gram_schmidt` algorithm.
    product
        Inner product |Operator| w.r.t. which to orthonormalize.
    riesz_representatives
        If `True`, compute Riesz representatives of the vectors in `image` before
        orthonormalizing. (Useful for norm computation when the range of the
        `operators` is a dual space.)

    Returns
    -------
    The |VectorArray| `image`.

    Raises
    ------
    ImageCollectionError
        Is raised when for a given |Operator| no image estimate is possible.
    """
    assert operators or vectors
    domain_space = operators[0].source if operators else None
    image_space = operators[0].range if operators \
        else vectors[0].space if isinstance(vectors[0], VectorArrayInterface) \
        else vectors[0].range if vectors[0].range != NumpyVectorSpace(1) \
        else vectors[0].source
    assert all(
        isinstance(v, VectorArrayInterface) and v in image_space or
        v.source == NumpyVectorSpace(1) and v.range == image_space and v.linear or
        v.range == NumpyVectorSpace(1) and v.source == image_space and v.linear
        for v in vectors
    )
    assert all(op.source == domain_space and op.range == image_space for op in operators)
    assert domain is None or domain_space is None or domain in domain_space
    assert product is None or product.source == product.range == image_space

    def collect_operator_ranges(op, source, image):
        if isinstance(op, (LincombOperator, SelectionOperator)):
            for o in op.operators:
                collect_operator_ranges(o, source, image)
        elif isinstance(op, EmpiricalInterpolatedOperator):
            if hasattr(op, 'collateral_basis') and not extends:
                image.append(op.collateral_basis)
        elif isinstance(op, Concatenation):
            firstrange = op.first.range.empty()
            collect_operator_ranges(op.first, source, firstrange)
            collect_operator_ranges(op.second, firstrange, image)
        elif op.linear and not op.parametric:
            image.append(op.apply(source))
        else:
            raise ImageCollectionError(op)

    def collect_vector_ranges(op, image):
        if isinstance(op, (LincombOperator, SelectionOperator)):
            for o in op.operators:
                collect_vector_ranges(o, image)
        elif isinstance(op, AdjointOperator):
            if op.source not in image_space:
                raise ImageCollectionError(op)  # Not implemented
            operator = Concatenation(op.range_product, op.operator) if op.range_product else op.operator
            collect_operator_ranges(operator, NumpyVectorArray(np.ones(1)), image)
        elif op.linear and not op.parametric:
            image.append(op.as_vector())
        else:
            raise ImageCollectionError(op)

    image = image_space.empty()
    if not extends:
#.........这里部分代码省略.........

def estimate_image_hierarchical(operators=(), vectors=(), domain=None, extends=None,
                                orthonormalize=True, product=None, riesz_representatives=False):
    """Estimate the image of given operators for all mu.

    This is an extended version of :func:`estimate_image`, which calls
    :func:`estimate_image` individually for each vector of `domain`.

    As a result, the vectors in the returned `image` |VectorArray| will
    be ordered by the `domain` vector they correspond to (starting with
    vectors which correspond to the `functionals` and to the |Operators| for
    which the image is estimated independently from `domain`).

    This function also returns an `image_dims` list, such that the first
    `image_dims[i+1]` vectors of `image` correspond to the first `i`
    vectors of `domain` (the first `image_dims[0]` vectors correspond
    to `vectors` and to the |Operators| with fixed image estimate).

    Parameters
    ----------
    operators
        See :func:`estimate_image`.
    vectors
        See :func:`estimate_image`.
    domain
        See :func:`estimate_image`.
    extends
        When additional vectors have been appended to the `domain` |VectorArray|
        after :func:`estimate_image_hierarchical` has been called, and
        :func:`estimate_image_hierarchical` shall be called again for the extended
        `domain` array, `extends` can be set to `(image, image_dims)`, where
        `image`, `image_dims` are the return values of the last
        :func:`estimate_image_hierarchical` call. The old `domain` vectors will
        then be skipped during computation and `image`, `image_dims` will be
        modified in-place.
    orthonormalize
        See :func:`estimate_image`.
    product
        See :func:`estimate_image`.
    riesz_representatives
        See :func:`estimate_image`.

    Returns
    -------
    image
        See above.
    image_dims
        See above.

    Raises
    ------
    ImageCollectionError
        Is raised when for a given |Operator| no image estimate is possible.
    """
    assert operators or vectors
    domain_space = operators[0].source if operators else None
    image_space = operators[0].range if operators \
        else vectors[0].space if isinstance(vectors[0], VectorArrayInterface) \
        else vectors[0].range if vectors[0].range != NumpyVectorSpace(1) \
        else vectors[0].source
    assert all(
        isinstance(v, VectorArrayInterface) and v in image_space or
        v.source == NumpyVectorSpace(1) and v.range == image_space and v.linear or
        v.range == NumpyVectorSpace(1) and v.source == image_space and v.linear
        for v in vectors
    )
    assert all(op.source == domain_space and op.range == image_space for op in operators)
    assert domain is None or domain_space is None or domain in domain_space
    assert product is None or product.source == product.range == image_space
    assert extends is None or len(extends) == 2

    logger = getLogger('pymor.algorithms.image.estimate_image_hierarchical')

    if operators and domain is None:
        domain = domain_space.empty()

    if extends:
        image = extends[0]
        image_dims = extends[1]
        ind_range = list(range(len(image_dims) - 1, len(domain))) if operators else list(range(len(image_dims) - 1, 0))
    else:
        image = image_space.empty()
        image_dims = []
        ind_range = list(range(-1, len(domain))) if operators else [-1]

    for i in ind_range:
        logger.info('Estimating image for basis vector {} ...'.format(i))
        if i == -1:
            new_image = estimate_image(operators, vectors, None, extends=False,
                                       orthonormalize=False, product=product,
                                       riesz_representatives=riesz_representatives)
        else:
            new_image = estimate_image(operators, [], domain.copy(i), extends=True,
                                       orthonormalize=False, product=product,
                                       riesz_representatives=riesz_representatives)

        gram_schmidt_offset = len(image)
        image.append(new_image, remove_from_other=True)
        if orthonormalize:
            with logger.block('Orthonormalizing ...'):
                gram_schmidt(image, offset=gram_schmidt_offset, product=product, copy=False)
#.........这里部分代码省略.........

def estimate_image(operators=(), vectors=(),
                   domain=None, extends=False, orthonormalize=True, product=None,
                   riesz_representatives=False):
    """Estimate the image of given |Operators| for all mu.

    Let `operators` be a list of |Operators| with common source and range, and let
    `vectors` be a list of |VectorArrays| or vector-like |Operators| in the range
    of these operators. Given a |VectorArray| `domain` of vectors in the source of the
    operators, this algorithms determines a |VectorArray| `image` of range vectors
    such that the linear span of `image` contains:

    - `op.apply(U, mu=mu)` for all operators `op` in `operators`, for all possible |Parameters|
      `mu` and for all |VectorArrays| `U` contained in the linear span of `domain`,
    - `U` for all |VectorArrays| in `vectors`,
    - `v.as_range_array(mu)` for all |Operators| in `vectors` and all possible |Parameters| `mu`.

    The algorithm will try to choose `image` as small as possible. However, no optimality
    is guaranteed. The image estimation algorithm is specified by :class:`CollectOperatorRangeRules`
    and :class:`CollectVectorRangeRules`.

    Parameters
    ----------
    operators
        See above.
    vectors
        See above.
    domain
        See above. If `None`, an empty `domain` |VectorArray| is assumed.
    extends
        For some operators, e.g. |EmpiricalInterpolatedOperator|, as well as for all
        elements of `vectors`, `image` is estimated independently from the choice of
        `domain`.  If `extends` is `True`, such operators are ignored. (This is useful
        in case these vectors have already been obtained by earlier calls to this
        function.)
    orthonormalize
        Compute an orthonormal basis for the linear span of `image` using the
        :func:`~pymor.algorithms.gram_schmidt.gram_schmidt` algorithm.
    product
        Inner product |Operator| w.r.t. which to orthonormalize.
    riesz_representatives
        If `True`, compute Riesz representatives of the vectors in `image` before
        orthonormalizing (useful for dual norm computation when the range of the
        `operators` is a dual space).

    Returns
    -------
    The |VectorArray| `image`.

    Raises
    ------
    ImageCollectionError
        Is raised when for a given |Operator| no image estimate is possible.
    """
    assert operators or vectors
    domain_space = operators[0].source if operators else None
    image_space = operators[0].range if operators \
        else vectors[0].space if isinstance(vectors[0], VectorArrayInterface) \
        else vectors[0].range
    assert all(op.source == domain_space and op.range == image_space for op in operators)
    assert all(
        isinstance(v, VectorArrayInterface) and (
            v in image_space
        ) or
        isinstance(v, OperatorInterface) and (
            v.range == image_space and isinstance(v.source, NumpyVectorSpace) and v.linear
        )
        for v in vectors
    )
    assert domain is None or domain_space is None or domain in domain_space
    assert product is None or product.source == product.range == image_space

    image = image_space.empty()
    if not extends:
        rules = CollectVectorRangeRules(image)
        for v in vectors:
            try:
                rules.apply(v)
            except NoMatchingRuleError as e:
                raise ImageCollectionError(e.obj)

    if operators and domain is None:
        domain = domain_space.empty()
    for op in operators:
        rules = CollectOperatorRangeRules(domain, image, extends)
        try:
            rules.apply(op)
        except NoMatchingRuleError as e:
            raise ImageCollectionError(e.obj)

    if riesz_representatives and product:
        image = product.apply_inverse(image)

    if orthonormalize:
        gram_schmidt(image, product=product, copy=False)

    return image

def adaptive_rrf(A, source_product=None, range_product=None, tol=1e-4,
                 failure_tolerance=1e-15, num_testvecs=20, lambda_min=None, iscomplex=False):
    r"""Adaptive randomized range approximation of `A`.

    This is an implementation of Algorithm 1 in [BS18]_.

    Given the |Operator| `A`, the return value of this method is the |VectorArray|
    `B` with the property

    .. math::
        \Vert A - P_{span(B)} A \Vert \leq tol

    with a failure probability smaller than `failure_tolerance`, where the norm denotes the
    operator norm. The inner product of the range of `A` is given by `range_product` and
    the inner product of the source of `A` is given by `source_product`.

    Parameters
    ----------
    A
        The |Operator| A.
    source_product
        Inner product |Operator| of the source of A.
    range_product
        Inner product |Operator| of the range of A.
    tol
        Error tolerance for the algorithm.
    failure_tolerance
        Maximum failure probability.
    num_testvecs
        Number of test vectors.
    lambda_min
        The smallest eigenvalue of source_product.
        If `None`, the smallest eigenvalue is computed using scipy.
    iscomplex
        If `True`, the random vectors are chosen complex.

    Returns
    -------
    B
        |VectorArray| which contains the basis, whose span approximates the range of A.
    """

    assert source_product is None or isinstance(source_product, OperatorInterface)
    assert range_product is None or isinstance(range_product, OperatorInterface)
    assert isinstance(A, OperatorInterface)

    B = A.range.empty()

    R = A.source.random(num_testvecs, distribution='normal')
    if iscomplex:
        R += 1j*A.source.random(num_testvecs, distribution='normal')

    if source_product is None:
        lambda_min = 1
    elif lambda_min is None:
        def mv(v):
            return source_product.apply(source_product.source.from_numpy(v)).to_numpy()

        def mvinv(v):
            return source_product.apply_inverse(source_product.range.from_numpy(v)).to_numpy()
        L = LinearOperator((source_product.source.dim, source_product.range.dim), matvec=mv)
        Linv = LinearOperator((source_product.range.dim, source_product.source.dim), matvec=mvinv)
        lambda_min = eigsh(L, sigma=0, which="LM", return_eigenvectors=False, k=1, OPinv=Linv)[0]

    testfail = failure_tolerance / min(A.source.dim, A.range.dim)
    testlimit = np.sqrt(2. * lambda_min) * erfinv(testfail**(1. / num_testvecs)) * tol
    maxnorm = np.inf
    M = A.apply(R)

    while(maxnorm > testlimit):
        basis_length = len(B)
        v = A.source.random(distribution='normal')
        if iscomplex:
            v += 1j*A.source.random(distribution='normal')
        B.append(A.apply(v))
        gram_schmidt(B, range_product, atol=0, rtol=0, offset=basis_length, copy=False)
        M -= B.lincomb(B.inner(M, range_product).T)
        maxnorm = np.max(M.norm(range_product))

    return B