def fit(self, X, y, X_val=None, y_val=None):
        """Train a network with the quasi-Newton method.
        
        Args:
            X (np.array of float): feature matrix for training
            y (np.array of float): target values for training
            X_val (np.array of float): feature matrix for validation
            y_val (np.array of float): target values for validation
        """
        y = y.reshape((len(y), 1))

        if sparse.issparse(X):
            X = X.tocsr()

        if X_val is not None:
            n_val = len(y_val)
            y_val = y_val.reshape((n_val, 1))

        # Set initial weights randomly.
        self.i = X.shape[1]
        self.l1 = self.l1 / self.i
        self.w = (np.random.rand((self.i + 2) * self.h + 1) - .5) * 1e-6
        self.w_opt = self.w
        self.n_opt = 0

        logger.info('training ...')
        n_obs = X.shape[0]
        batch = self.b
        n_epoch = self.n
        idx = range(n_obs)
        self.auc_opt = .5

        start = time.time()
        print('\tEPOCH TRAIN     VALID     BEST      TIME (m)')
        print('\t--------------------------------------------')

        # Before training
        p = self.predict_raw(X)
        auc = roc_auc_score(y, p)
        auc_val = auc
        if X_val is not None:
            p_val = self.predict_raw(X_val)
            auc_val = roc_auc_score(y_val, p_val)

        print('\t{:3d}:  {:.6f}  {:.6f}  {:.6f}  {:.2f}'.format(
              0, auc, auc_val, self.auc_opt,
              (time.time() - start) / SEC_PER_MIN))
     
        # Use 'while' instead of 'for' to increase n_epoch if the validation
        # error keeps improving at the end of n_epoch 
        epoch = 1
        while epoch <= n_epoch:
            # Shuffle inputs every epoch - it helps avoiding the local optimum
            # when batch < n_obs.
            np.random.shuffle(idx)

            # Find the optimal weights for batch input examples.
            # If batch == 1, it's the stochastic optimization, which is slow
            # but uses minimal memory.  If batch == n_obs, it's the batch
            # optimization, which is fast but uses maximum memory.
            # Otherwise, it's the mini-batch optimization, which balances the
            # speed and space trade-offs.
            for i in range(int(n_obs / batch) + 1):
                if (i + 1) * batch > n_obs:
                    sub_idx = idx[batch * i:n_obs]
                else:
                    sub_idx = idx[batch * i:batch * (i + 1)]

                x = X[sub_idx]
                neg_idx = [n_idx for n_idx, n_y in enumerate(y[sub_idx]) if n_y == 0.]
                pos_idx = [p_idx for p_idx, p_y in enumerate(y[sub_idx]) if p_y == 1.]
                x0 = x[neg_idx]
                x1 = x[pos_idx]
                # Update weights to minimize the cost function using the
                # quasi-Newton method (L-BFGS-B), where:
                #   func -- cost function
                #   jac -- jacobian (derivative of the cost function)
                #   maxiter -- number of iterations for L-BFGS-B
                ret = minimize(self.func,
                               self.w,
                               args=(x0, x1),
                               method='L-BFGS-B',
                               jac=self.fprime,
                               options={'maxiter': 5})
                self.w = ret.x

            p = self.predict_raw(X)
            auc = roc_auc_score(y, p)
            auc_val = auc

            if X_val is not None:
                p_val = self.predict_raw(X_val)
                auc_val = roc_auc_score(y_val, p_val)

                if auc_val > self.auc_opt:
                    self.auc_opt = auc_val
                    self.w_opt = self.w
                    self.n_opt = epoch

                    # If validation auc is still improving after n_epoch,
                    # try 10 more epochs
                    if epoch == n_epoch:
                        n_epoch += 5

            print('\t{:3d}:  {:.6f}  {:.6f}  {:.6f}  {:.2f}'.format(
                  epoch, auc, auc_val, self.auc_opt,
                  (time.time() - start) / SEC_PER_MIN))

            epoch += 1

        if X_val is not None:
            print('Optimal epoch is {0} ({1:.6f})'.format(self.n_opt,
                                                          self.auc_opt))
            self.w = self.w_opt

        logger.info('done training')
Train a network with the quasi-Newton method.
        
        Args:
            X (np.array of float): feature matrix for training
            y (np.array of float): target values for training
            X_val (np.array of float): feature matrix for validation
            y_val (np.array of float): target values for validation