def _conjugate_sinkhorn(self, x: np.ndarray, grad: np.ndarray, cost_matrix: np.ndarray) -> np.ndarray:
        """
        The conjugate sinkhorn_optimizer.

        :param x: Current adversarial examples.
        :param grad: The loss gradients.
        :param cost_matrix: A non-negative cost matrix.
        :return: Adversarial examples.
        """
        # Normalize inputs
        normalization = x.reshape(x.shape[0], -1).sum(-1).reshape(x.shape[0], 1, 1, 1)
        x = x.copy() / normalization

        # Dimension size for each example
        m = np.prod(x.shape[1:])

        # Initialize
        alpha = np.log(np.ones(x.shape) / m) + 0.5
        exp_alpha = np.exp(-alpha)

        beta = -self.regularization * grad
        beta = beta.astype(np.float64)
        exp_beta = np.exp(-beta)

        # Check for overflow
        if (exp_beta == np.inf).any():
            raise ValueError("Overflow error in `_conjugate_sinkhorn` for exponential beta.")

        cost_matrix_new = cost_matrix.copy() + 1
        cost_matrix_new = np.expand_dims(np.expand_dims(cost_matrix_new, 0), 0)

        i_nonzero = self._batch_dot(x, self._local_transport(cost_matrix_new, grad, self.kernel_size)) != 0
        i_nonzero_ = np.zeros(alpha.shape).astype(bool)
        i_nonzero_[:, :, :, :] = np.expand_dims(np.expand_dims(np.expand_dims(i_nonzero, -1), -1), -1)

        psi = np.ones(x.shape[0])

        var_k = np.expand_dims(np.expand_dims(np.expand_dims(psi, -1), -1), -1)
        var_k = np.exp(-var_k * cost_matrix - 1)

        convergence = np.array([-np.inf])

        for _ in range(self.conjugate_sinkhorn_max_iter):
            # Block coordinate descent iterates
            x[x == 0.0] = EPS_LOG  # Prevent divide by zero in np.log
            alpha[i_nonzero_] = (np.log(self._local_transport(var_k, exp_beta, self.kernel_size)) - np.log(x))[
                i_nonzero_
            ]
            exp_alpha = np.exp(-alpha)

            # Newton step
            var_g = -self.eps_step + self._batch_dot(
                exp_alpha, self._local_transport(cost_matrix * var_k, exp_beta, self.kernel_size)
            )

            var_h = -self._batch_dot(
                exp_alpha, self._local_transport(cost_matrix * cost_matrix * var_k, exp_beta, self.kernel_size)
            )

            delta = var_g / var_h

            # Ensure psi >= 0
            tmp = np.ones(delta.shape)
            neg = psi - tmp * delta < 0

            while neg.any() and np.min(tmp) > 1e-2:
                tmp[neg] /= 2
                neg = psi - tmp * delta < 0

            psi[i_nonzero] = np.maximum(psi - tmp * delta, 0)[i_nonzero]

            # Update K
            var_k = np.expand_dims(np.expand_dims(np.expand_dims(psi, -1), -1), -1)
            var_k = np.exp(-var_k * cost_matrix - 1)

            # Check for convergence
            next_convergence = self._conjugated_sinkhorn_evaluation(x, alpha, exp_alpha, exp_beta, psi, var_k)

            if (np.abs(convergence - next_convergence) <= 1e-4 + 1e-4 * np.abs(next_convergence)).all():
                break

            convergence = next_convergence

        result = exp_beta * self._local_transport(var_k, exp_alpha, self.kernel_size)
        result[~i_nonzero] = 0
        result *= normalization

        return result

        The conjugate sinkhorn_optimizer.

        :param x: Current adversarial examples.
        :param grad: The loss gradients.
        :param cost_matrix: A non-negative cost matrix.
        :return: Adversarial examples.
        
art/attacks/evasion/wasserstein.py