def optimize_delta_bisection(self, data, skipstep=1, verbose=None):
        """Find optimal ridge penalty with bisection search.
        
        Parameters
        ----------
        data : array, shape (n_trials, n_channels, n_samples)
            Epoched data set. At least two trials are required.
        skipstep : int, optional
            Speed up calculation by skipping samples during cost function
            calculation.
            
        Returns
        -------
        self : :class:`VAR`
            The :class:`VAR` object to facilitate method chaining (see usage
            example).
        """
        data = atleast_3d(data)
        if data.shape[0] < 2:
            raise ValueError("At least two trials are required.")

        if verbose is None:
            verbose = config.getboolean('scot', 'verbose')

        maxsteps = 10
        maxdelta = 1e50

        a = -10
        b = 10

        trform = lambda x: np.sqrt(np.exp(x))

        msge = _get_msge_with_gradient_func(data.shape, self.p)

        ja, ka = msge(data, trform(a), self.xvschema, skipstep, self.p)
        jb, kb = msge(data, trform(b), self.xvschema, skipstep, self.p)

        # before starting the real bisection, assure the interval contains 0
        while np.sign(ka) == np.sign(kb):
            if verbose:
                print('Bisection initial interval (%f,%f) does not contain 0. '
                      'New interval: (%f,%f)' % (a, b, a * 2, b * 2))
            a *= 2
            b *= 2
            ja, ka = msge(data, trform(a), self.xvschema, skipstep, self.p)
            jb, kb = msge(data, trform(b), self.xvschema, skipstep, self.p)

            if trform(b) >= maxdelta:
                if verbose:
                    print('Bisection: could not find initial interval.')
                    print(' ********* Delta set to zero! ************ ')
                return 0

        nsteps = 0

        while nsteps < maxsteps:
            # point where the line between a and b crosses zero
            # this is not very stable!
            #c = a + (b-a) * np.abs(ka) / np.abs(kb-ka)
            c = (a + b) / 2
            j, k = msge(data, trform(c), self.xvschema, skipstep, self.p)
            if np.sign(k) == np.sign(ka):
                a, ka = c, k
            else:
                b, kb = c, k

            nsteps += 1
            tmp = trform([a, b, a + (b - a) * np.abs(ka) / np.abs(kb - ka)])
            if verbose:
                print('%d Bisection Interval: %f - %f, (projected: %f)' %
                      (nsteps, tmp[0], tmp[1], tmp[2]))

        self.delta = trform(a + (b - a) * np.abs(ka) / np.abs(kb - ka))
        if verbose:
            print('Final point: %f' % self.delta)
        return self