I am currently studying this Nesterov's paper for project purposes, and I am trying to figure out how the smoothing and the minimization algorithm works


I have tried looking at the example provided on minimax strategies for matrix games,

$$ \Delta_n=\{x\in \mathbb{R}^n:\,x\ge0,\, \sum_{i=1}^{n}x^{(i)}=1\}\\ \min_{x\in \Delta_n}\,\left[{\langle c,x\rangle}_1+\max_{1\le j\le m}\left({\langle a_j,x\rangle}_1+b^{(j)}\right)\right] $$

and after smoothing of the objective function, I think the problem now looks like

$$ \min_{x\in \Delta_n}\,\left[{\langle c,x\rangle}_1+\mu\ln\left(\dfrac{1}{m}\sum_{j=1}^{m}e^{[{\langle a_j,x\rangle}+b^{(j)}]/\mu}\right)\right] $$

This example and the others in the paper are too complicated for me (I'm new to this, and confused). So I was hoping that someone who is familiar with the paper could give some guidance or provide much simpler examples with non-smooth objectives so that I can try smoothing it and then apply the algorithm above to get the $\epsilon$-solution.

  • 2
    $\begingroup$ Look, I've built my academic career on a foundation laid in large part by Nesterov and Nemirovskii. I have nothing but the highest respect for their work. But you are right, their work is complicated, and their papers are particularly difficult to read, even for those of us who have spent a lot of time trying. Do you really need to digest this paper, or do you simply need to solve a particular class of convex problems? If it is the latter let me suggest that you should be asking about that instead. $\endgroup$ Jun 26, 2015 at 17:05
  • 2
    $\begingroup$ Where precisely are you stuck in the paper ? What don't you understand ? What problem(s) do you have in mind ? As @MichaelGrant already noted, chances are, you can do without that paper. Without more detail it's difficult to help you. BTW, the underlying ideas in the paper are really simple, except that they're written in an "unclear / unnatural / non-standard" manner. Mind the quotes. This is characteristic of Nesterov. Anecdote: Nesterov "essentially" invented the so-called FISTA algorithm, but due to "unclear" writing, the credit went to other people (beck, teboulle ...) ... $\endgroup$
    – dohmatob
    Jun 27, 2015 at 9:21
  • 1
    $\begingroup$ You may find this somewhat simplification of the paper easier to read optimization-online.org/DB_FILE/2009/04/2286.pdf. They apply the method to compressed sensing problems, etc. $\endgroup$
    – dohmatob
    Jun 27, 2015 at 9:25
  • $\begingroup$ @MichaelGrant I'm actually planning to do my honours project on optimization. My supervisor suggested me to look at this paper and then find directions to work with. There're too many parts in the paper that I don't understand and hence I'm planning to start off with trying out some concrete examples before delving into the details. $\endgroup$
    – Sapphire
    Jun 27, 2015 at 15:32
  • 1
    $\begingroup$ Check out Vandenberghe's 236c notes, such as the slides on smoothing and the slides on fast proximal gradient algorithms. Much of the literature has been digested and presented in a short, clear way. At the end of each chapter is a list of references, including to papers by nesterov. $\endgroup$
    – littleO
    Jun 28, 2015 at 7:03

1 Answer 1


Details of Nesterov smoothing for computing Nash equilibria in matrix games

Let $A \in \mathbb{R}^{m\times n}$, $c \in \mathbb{R}^n$ and $b \in\mathbb{R}^m$, and consider a matrix game with payoff function $(x, u) \mapsto \langle Ax, u\rangle + \langle c, x\rangle + \langle b, u\rangle$. The Nash equilibrium problem is \begin{eqnarray} \underset{x \in \Delta_n}{\text{minimize }}\underset{u \in \Delta_m}{\text{maximize }}\langle Ax, u\rangle + \langle c, x\rangle + \langle b, u\rangle. \end{eqnarray}

From the "min" player's point of view, this problem can be re-written as \begin{eqnarray} \underset{x \in \Delta_n}{\text{minimize }}f(x), \end{eqnarray}

where the proper convex function $f: \mathbb{R}^n \rightarrow \mathbb{R}$ is defined by \begin{eqnarray} f(x) := \hat{f}(x) + \underset{u \in \Delta_m}{\text{max }}\langle Ax, u\rangle - \hat{\phi}(u),\hspace{.5em}\hat{f}(x) := \langle c, x\rangle,\hspace{.5em}\hat{\phi}(u) := \langle b, u\rangle. \end{eqnarray} Of course, the dual problem is \begin{eqnarray} \underset{u \in \Delta_m}{\text{maximize }}\phi(u), \end{eqnarray}

where the proper concave function $\phi: \mathbb{R}^m \rightarrow \mathbb{R}$ is defined by \begin{eqnarray} \phi(u) := -\hat{\phi}(u) + \underset{x \in \Delta_n}{\text{min }}\langle Ax, u\rangle + \hat{f}(x). \end{eqnarray} Now, for a positive scalar $\mu$, consider a smoothed version of $f$ \begin{eqnarray} \bar{f}_\mu(x) := \hat{f}(x) + \underbrace{\underset{u \in \Delta_m}{\text{max }}\langle Ax, u\rangle - \hat{\phi}(u) - \mu d_2(u)}_{f_\mu(x)}, \end{eqnarray} where $d_2$ is a prox-function (see Nesterov's paper for definition) for the simplex $\Delta_m$, with prox-center $\bar{u} \in \Delta_m$. For simplicity, take \begin{eqnarray} d_2(u) := \frac{1}{2}\|u - \bar{u}\|^2 \ge 0 = d_2(\bar{u}), \forall u \in \Delta_m, \end{eqnarray} where $\bar{u} := (1/m, 1/m, ..., 1/m)$ is the centroid of $\Delta_m$. Define a prox-function $d_1$ for $\Delta_n$ analogously. Note that $\sigma_1 = \sigma_2 = 1$ since the $d_j$'s are $1$-strongly convex. Also, noting that the $d_j$'s attain their maximum at the kinks of the respective simplexes, one computes \begin{eqnarray*} \begin{split} D_1 := \underset{x \in \Delta_n}{\text{max }}d_1(x) &= \frac{1}{2} \times \text{ squared distance between any kink and }\bar{x}\\ &= \frac{1}{2}\left((1-1/n)^2 + (n-1)(1/n)^2\right) = (1 - 1 / n) / 2 \end{split} \end{eqnarray*} and similarly $D_2 := \underset{u \in \Delta_m}{\text{max }}d_2(u) = (1 - 1 / m) / 2$.

Note that the functions $(f_\mu)_{\mu > 0}$ increase point-wise to $f$ in the limit $\mu \rightarrow 0^+$.

Now, let us prepare the ingredients necessary for implementing the algorithm.

Definition. Given a closed convex subset $C \subseteq \mathbb{R}^m$, and a point $x \in \mathbb{R}^m$, recall the definition of the orthogonal projection of $x$ unto $C$ \begin{eqnarray} proj_C(x) := \underset{z \in C}{\text{argmin }}\frac{1}{2}\|z-x\|^2. \end{eqnarray}

Geometrically, $proj_C(x)$ is the unique point of $C$ which is closest to $x$, and so $proj_C$ is a well-defined function from $\mathbb{R}^m$ to $C$. For example, if $C$ is the nonnegative $m$th-orthant $\mathbb{R}_+^m$, then $proj_C(x) \equiv (x)_+$, the component-wise maximum of $x$ and $0$.

Step 1: Call to the oracle, or computation of $f(x_k)$ and $\nabla f(x_k)$. For any $x \in \mathbb{R}^n$, define $v_\mu(x) := (Ax - b) / \mu$ and $u_\mu(x) := proj_{\Delta_m}(\bar{u} + v_\mu(x))$. One computes \begin{eqnarray*} \begin{split} f_\mu(x) := \underset{u \in \Delta_m}{\text{max }}\langle Ax, u\rangle - \hat{\phi}(u) - \mu d_2(u) &= -\underset{u \in \Delta_m}{\text{min }}\frac{1}{2}\mu\|u-\bar{u}\|^2 + \langle b - Ax, u\rangle\\ &= \frac{1}{2}\mu\|v_\mu(x)\|^2 - \underset{u \in \Delta_m}{\text{min }}\frac{1}{2}\mu\|u - (\bar{u} + v_\mu(x))\|^2, \end{split} \end{eqnarray*} by completing the square in $u$ and ignoring constant terms. One recognizes the minimization problem on the right-hand side as the orthogonal projection of the point $\bar{u} + v_\mu(x)$ unto the simplex $\Delta_m$. It is not difficult to show that $f_\mu$ is smooth with gradient $\nabla f_\mu: x \mapsto A^Tu_\mu(x)$ (one can invoke Danskin's theorem, for example) and $L_\mu: = \frac{1}{\mu\sigma_2}\|A\|^2$ is a Lipschitz constant for the latter. Thus, $\bar{f}_\mu$ is smooth and one computes \begin{eqnarray} \nabla \bar{f}_\mu(x_k) = \nabla \hat{f}(x_k) + \nabla f_\mu(x_k) = c + A^*u_\mu(x_k). \end{eqnarray}

Moreover, $L_\mu$ defined above is also a Lipschitz constant for $\nabla \bar{f}_\mu$.

Step 3: Computation of $z_k$. For any $s \in \mathbb{R}^n$, one has \begin{eqnarray*} \underset{x \in \Delta_n}{\text{argmin }}Ld_1(x) + \langle s, x\rangle = \underset{x \in \Delta_n}{\text{argmin }}\frac{1}{2}L_\mu\|x-\bar{x}\|^2 + \langle s, x\rangle = proj_{\Delta_n}\left(\bar{x} - \frac{1}{L_\mu}s\right), \end{eqnarray*} by completing the square in $x$ and ignoring constant terms. Thus, with $s = \sum_{i=0}^k\frac{i+1}{2}\nabla \bar{f}_\mu(x_i)$, we obtain the update rule \begin{eqnarray} z_k = proj_{\Delta_n}\left(\bar{x} - \frac{1}{L_\mu}\sum_{i=0}^k\frac{i+1}{2}\nabla \bar{f}_\mu(x_i)\right). \end{eqnarray}

Step 2: Computation of $y_k$. Similarly to the previous computation, one obtains the update rule \begin{eqnarray} y_k = T_{\Delta_n}(x_k) := \underset{y \in \Delta_n}{\text{argmin }}\langle \nabla \bar{f}_\mu(x), y\rangle + \frac{1}{2}L_\mu\|y-x\|^2 = proj_{\Delta_n}\left(x_k - \frac{1}{L_\mu}\nabla \bar{f}_\mu (x_k)\right). \end{eqnarray}

Stopping criterion. The algorithm is stopped once the primal-dual gap \begin{eqnarray} \begin{split} &f(x_k) - \phi(u_\mu(x_k))=\\ &\left(\langle c, x_k \rangle + \underset{u \in \Delta_m}{\text{max }}\langle Ax_k, u\rangle - \langle b, u\rangle\right) - \left(-\langle b, u_\mu(x_k)\rangle + \underset{z \in \Delta_n}{\text{min }}\langle Az, u_\mu(x_k)\rangle + \langle c, z\rangle\right)\\ &= \langle c, x_k \rangle + \underset{1 \le j \le m}{\text{max }}(Ax_k-b)_j + \langle b, u_\mu(x_k)\rangle - \underset{1 \le i \le n}{\text{min }}(A^Tu_\mu(x_k) + c)_i \end{split} \end{eqnarray} is below a tolerance $\epsilon > 0$ (say $\epsilon \sim 10^{-4}$). Note that in the above computation, we have used the fact that for any $r \in \mathbb{R}^m$, one has \begin{eqnarray} \begin{split} \underset{u \in \Delta_m}{\text{max }}\langle r, u\rangle &= \text{max of }\langle r, .\rangle\text{ on the boundary of }\Delta_m \\ &= \text{max of }\langle r, .\rangle\text{ on the line segment joining any two distinct kinks of } \Delta_m \\ &= \underset{1 \le i < j \le m}{\text{max }}\text{ }\underset{0 \le t \le 1}{\text{max }}tr_i + (1-t)r_j = \underset{1 \le i \le m}{\text{max }}r_i. \end{split} \end{eqnarray}

Algorithm parameters. In accord with Nesterov's recommendation, we may take \begin{eqnarray} \mu = \frac{\epsilon}{2D_2} = \frac{\epsilon}{2}, \text{ and }L_\mu=\frac{\|A\|^2}{\sigma_2 \mu} = \frac{\|A\|^2}{\mu}. \end{eqnarray}

Implementation The (quick-and-dirty) script below implements the algorithm presented above

Solving matrix games via Nesterov smoothing. Quick-and-dirty code!
# Author: Elvis DOHMATOB <[email protected]>

from math import sqrt
import numpy as np
from scipy import linalg

def proj_simplex(v, z=1.):
    """Projects v unto the simplex {x >= 0, x_0 + x_1 + ... x_n = z}.

    The method is John Duchi's O (n log n) Algorithm 1.
    # deterministic O(n log n)
    u = np.sort(v)[::-1]  # sort v in increasing order
    aux = (np.cumsum(u) - z) / np.arange(1., len(v) + 1.)
    return np.maximum(v - aux[np.nonzero(u > aux)[0][-1]], 0.)

def test_proj_simplex():
    v = np.array([1.1, 0., 0.])
    np.testing.assert_array_equal(proj_simplex(v), [1., 0., 0.])

    v = np.array([0., 0., 0.])
    np.testing.assert_array_equal(proj_simplex(v), [1. / 3, 1. / 3, 1. / 3])

def nesterov_ne(A, c=None, b=None, epsilon=1e-3, max_iter=np.inf,
                dynamic_mu=False, mu_init=1e-2, mu_factor=.9):
    """Computes an approx Nash equilibrium for matrix game via Nesterov
    smoothing [1].

    Formally, the problem is

        minimize maximize <Ax, u> + <c, x> + <b, u>,
        x in S_n u in S_m

    where S_k is the k-simplex.

    A: ndarray, shape (m, n)
        Payoff matrix.

    c: ndarray, shape (n,), optional (default None)
        <c, x> is added to the payoff function <Ax, u>.

    b: ndarray, shape (n,), optional (default None)
        <b, u> is added to the payoff function <Ax, u>

    epsilon: positive float, optional (default 1e-3)
        Tolerance on primal-dual gap.

    max_iter: int, optional (default np.inf)
        Maximum number of iterations to run. If no value is specified,
        then it is inferred using Nesterov's formulae.

    dynamic_mu: boolean, optional (default False)
        If true, then the smoothing parameter mu will be set dynamically.

    mu_init: positive float, optional (default 1e-2)
        Initial value for the smoothing parameter mu.

    mu_factor: positive float, optional (default .9)
        Factor by which the smoothing parameter mu is shrunk at each update.

    [1] Y. Nesterov, "Smooth minimization of non-smooth functions."
    # misc
    m, n = A.shape
    if not c is None: assert len(c) == n
    if not b is None: assert len(b) == m
    x_ = (1. / n) * np.ones(n)
    u_ = (1. / n) * np.ones(m)
    norm_A = linalg.norm(A, 2)
    D_1 = .5 * (1. - 1. / n)
    D_2 = .5 * (1. - 1. / m)

    # set parameters using formula 4.8 of [1]
    max_iter = min(int(np.floor(4 * norm_A * sqrt(D_1 * D_2) / epsilon)),
    mu_ = epsilon / (2. * D_2)
    x = x_.copy()
    grad_acc = np.zeros_like(x)
    values = []
    gaps = []
    mu = mu_init if dynamic_mu else mu_
    print "mu =", mu
    q = 2
    for k in range(max_iter):
        # misc
        L = norm_A ** 2 / mu  # there is an error in the L formula from [1]
        stepsize = 1. / L

        # make call to oracle
        aux = A.dot(x)
        if not b is None: aux -= b
        v = aux / mu
        v += u_
        u = proj_simplex(v)
        grad = A.T.dot(u)
        if not c is None: grad += c
        grad_acc += .5 * (k + 1.) * grad

        # callback
        value = u.dot(A.dot(x))
        gap = aux.max() - grad.min()
        if not c is None: gap += c.dot(x)
        if not b is None: gap += b.dot(u)
        assert gap + 1e-10 >= 0., "The world is a weird place!"
        print "Iter %03i/%03i: game value <Ax, u> = %g, primal-dual gap=%g" % (
            k + 1, max_iter, value, gap)

        # check convergence
        if gap < epsilon:
            print "Converged (primal-dual gap < %g)." % epsilon

        # y update
        y = proj_simplex(x - stepsize * grad)

        # z update
        z = proj_simplex(x_ - stepsize * grad_acc)

        # x update
        factor = 2. / (k + 3.)
        x = factor * z
        x += (1. - factor) * y

        # decrease mu ?
        if dynamic_mu and mu > mu_ and k > 0 and k % q == 0:
            # the idea is to decrease mu at iterations k = 2, 4, 8, 16, 32, ...
            q *= 2
            mu *= mu_factor
            print "Decreasing mu to %g" % mu

    return x, u, values, gaps

if __name__ == "__main__":
    import matplotlib.pyplot as plt
    rng = np.random.RandomState(42)
    for game in ["random", "ferguson"]:
        if game == "ferguson":
            A = np.array([[-2., 3.], [3, -4]])
            b = c = None
            title = "Ferguson game"
            A = rng.randn(500, 500)
            c = rng.randn(A.shape[1])
            b = rng.randn(A.shape[0])
            title = "Random %i x %i game" % A.shape
        x, u, values, gaps = nesterov_ne(A, c=c, b=b, dynamic_mu=True)
        print "x =", x
        print "u =", u
        plt.figure(figsize=(13.5, 7))
        ax = plt.subplot2grid((1, 2), (0, 0))
        ax.axhline(values[-1], linestyle="--")
        ax.set_ylabel("value of game")
        ax = plt.subplot2grid((1, 2), (0, 1))
        ax.set_ylabel("primal-dual gap")
        plt.savefig("%s_gaps.png" % title.replace(" ", "_"))
        plt.savefig("%s_values.png" % title.replace(" ", "_"))

Demo. The above script demos the algorithm on a matrix game with random payoff matrix, and on a matrix game from a book by Ferguson. When run, the above script should output

Iter 180/12324: game value <Ax, u> = 0.0856024, primal-dual gap=4.0856
Iter 181/12324: game value <Ax, u> = 0.0858599, primal-dual gap=4.08586
Iter 182/12324: game value <Ax, u> = 0.0854715, primal-dual gap=4.08547
Iter 183/12324: game value <Ax, u> = 0.0844478, primal-dual gap=4.08445
Iter 184/12324: game value <Ax, u> = 0.0827988, primal-dual gap=4.08372
Iter 185/12324: game value <Ax, u> = 0.0837844, primal-dual gap=1.87479
Iter 186/12324: game value <Ax, u> = 0.0835924, primal-dual gap=2.08359
Iter 187/12324: game value <Ax, u> = 0.0833272, primal-dual gap=1.13781
Iter 188/12324: game value <Ax, u> = 0.0833215, primal-dual gap=0.677185
Iter 189/12324: game value <Ax, u> = 0.0833295, primal-dual gap=0.154766
Iter 190/12324: game value <Ax, u> = 0.0833338, primal-dual gap=0.01287
Iter 191/12324: game value <Ax, u> = 0.0833336, primal-dual gap=0.00710599
Iter 192/12324: game value <Ax, u> = 0.0833333, primal-dual gap=0.000174806
Converged (primal-dual gap < 0.001).
x = [ 0.58331946  0.41668054]
u = [ 0.58335442  0.41664558]

enter image description here enter image description here

Note that $7/12 \approx 0.58330147$ and $8\frac{1}{3} \text{ cents} = 0.08\overset{.}{3}$ dollars.

Limitations and extensions. There are a number of conceptual problems with this smoothing scheme. For example, how do we select the smoothing parameter $\mu$? If it's too small, then the scheme becomes ill-conditioned (Lipschitz constants explode, etc.). If it's too large, then the approximation is very bad. A possible workaround is to start with a reasonably big value of $\mu$ and decrease it gradually across iterations (as was done in Algorithm 1 of this paper, for example). In fact, Nesterov himself has proposed the "Excessive Gap Technique" (google it) as a more principled way to go about smoothing. EGT has been used by Andrew Gilpin as an alternative way to go about solving imcomplete information two-person zero-sum sequential games (In such games, the player's strategy profiles are no longer simplexes but more general polyhera.), the player's like Texas Hold'em.

I hope this helps you get the ball rolling. Let me know if you have more issues.


  • $\begingroup$ Thanks for the detailed explanation. I'll take some time to digest it first. BTW, I tried running the code and I got x = [ 0.58368601 0.41631399] and u = [ 0. 1.] instead. I've yet to examine the code though. $\endgroup$
    – Sapphire
    Jun 28, 2015 at 4:54
  • 1
    $\begingroup$ I've update my answer to clarify the points to raise. Let me summarize here: The differentiability of $f_\mu$ and the expression of its gradient follows from Danskin's theorem (google it). Next, the computation of $D_j$ is straightforward from the geometry of the simplex (the maximum of $d_j$ is attained at the kinks of the simplex). I forgot a 1/2 factor in the expression for $D_j$. Fixed. Finally, for steps 2&3, simply complete the square to see where the projection comes in... $\endgroup$
    – dohmatob
    Jun 29, 2015 at 0:33
  • 1
    $\begingroup$ Last but not the least, grad_acc (in code) = $\sum_{i=0}^k\frac{i+1}{2}\nabla f(x_i)$ (in paper). Nothing mysterious here ;) $\endgroup$
    – dohmatob
    Jun 29, 2015 at 0:40
  • 1
    $\begingroup$ BTW, grad_acc means "accumulated gradient". $\endgroup$
    – dohmatob
    Jun 29, 2015 at 23:21
  • 1
    $\begingroup$ This is a fundamental result in matrix-game theory saying: mixed strategy u is a best reponse to an opponent iff each component pure strategy is also a best strategy. I strongly encourage you to prove it as an exercise. If you'd rather skip the exercise, then lookup proposition 3.1 of maths.lse.ac.uk/personal/stengel/TEXTE/agt-stengel.pdf $\endgroup$
    – dohmatob
    Jul 7, 2015 at 13:54

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .