docs: def apply(self, a, return_Ya=False):
          r"""Apply the projection to an array.
  
          The computation is carried out without explicitly forming the
          matrix corresponding to the projection (which would be an array with
          ``shape==(N,N)``).
  
          See also :py:meth:`_apply`.
          """
          # is projection the zero operator?
          if self.V.shape[1] == 0:
              Pa = numpy.zeros(a.shape)
              if return_Ya:
                  return Pa, numpy.zeros((0, a.shape[1]))
              return Pa
          if return_Ya:
              x, Ya = self._apply(a, return_Ya=return_Ya)
          else:
              x = self._apply(a)
          for i in range(self.iterations-1):
              z = a - x
              w = self._apply(z)
              x = x + w
          if return_Ya:
              return x, Ya
          return x