Guest course Dec. 7, 2022

Structure identification: from nonsmooth optimization to data science

Franck Iutzeler

The goal of this course/lab is to explore the notion of structure identification in nonsmooth optimization and its applications in data science.

1 - Proximity operator and stability¶

The proximity operator is defined for a function $g:\mathbb{R}^n\to\mathbb{R}\cup +\infty$ and a point $x\in\mathbb{R}^n$ as the mapping $$ \mathbf{prox}_{\gamma g}(x) = \arg\min_{u\in\mathbb{R}^n} \left\{ g(u) + \frac{1}{2\gamma}\| x - u \|^2 \right\} $$

1.a Stability in 1D¶

Let us first consider the absolute function $g = |\cdot|$ in $\mathbb{R}$ and plot the function $$ u \mapsto s_g(u; \gamma,x) = |u| + \frac{1}{2\gamma}\| x - u \|^2 $$ for different values of $\gamma,x$.

1. Show that $$ \mathbf{prox}_{\gamma g}(x) = \left\{ \begin{array}{ll} x-\gamma & \text{if } x\geq \gamma \\ x+\gamma & \text{if } x\leq -\gamma \\ 0 & \text{if } -\gamma \leq x \leq \gamma \end{array} \right. $$ Notice how the points of non-differentiabiity of $g$ trap the minimizers of the problem.

This is a typical example of how nonsmoothness can be used to promote certain patterns in the solutions of minimization problems.

1.b Stability in nD¶

Let us first consider the $\ell_1$ norm, ie. $g: x \mapsto \|x\|_1 = \sum_{i=1}^n |x_i|$ in $\mathbb{R}^n$ (where $x_i$ is the $i$-th component of vector $x$).

2. How can one compute the proximity operator of the $\ell_1$ norm from the one of the absolute value?

In the first part of the lab, we will consider the classical lasso problem: \begin{align} \tag{lasso} \min_{x\in\mathbb{R}^n} \frac{1}{2}\|Ax-b\|^2 + \lambda \|x\|_1 \end{align}

3. Observe the level sets of the lasso problem implemented below. Run the code several times to generate different sensing matrices $A$ and vector $b = A x_0 + \text{noise}$. Notice how the solution of the problem often seem to have a null first coordinate, as $x_0$.

2 - Stucture identification & Optimization¶

In this section, we will again consider the lasso problem $$ \min_{x\in\mathbb{R}^n} \frac{1}{2} \|Ax-b\|^2 + \lambda \|x\|_1 $$ using the data generated below and we will compare the solutions obtained by several optimization methods.

In Section 2.1 to 2.5, you will be asked to solve this problem by different manners. Make sure to fill at least one method from 2.1,2.2,2.3 and one from 2.4,2.5. You will be asked to compare the solution in Section 2.6.

2.1 General purpose minimization¶

4. Solve the problem with a general purpose solver such as scipy.optimize.minimize https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html#scipy.optimize.minimize. Store the result in a variable xopt_GP.

2.2 Modelling software¶

5. Solve the problem with a modelling software (eg. cvxpy https://www.cvxpy.org/) or a QP solver. Store the result in a variable xopt_CVXPY.

2.3 Dedicated library¶

6. Solve the problem with a dedicated library as sklearn https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Lasso.html. Store the result in a variable xopt_sklearn.

2.4 Proximal Gradient Method¶

Recall, the iterations of the proximal gradient are of the form $$ x_{k+1} = \mathbf{prox}_{\gamma g}\left( x_k - \gamma \nabla f(x_k) \right) $$ where $\gamma$ is a stepsize to choose.

7. Solve the problem with a proximal gradient method. Store the result in a variable xopt_PG.

2.5 Proximal Splitting methods: Douglas-Rachford, ADMM, etc.¶

Use any proximal splitting method you wish, eg. the Douglas-Rachford method which reads \begin{align*} x_k &= \mathbf{prox}_{\gamma g}(u_k)\\ z_k &= \mathbf{prox}_{\gamma f}(2x_k-u_k)\\ u_{k+1} &= u_k + z_k - x_k \end{align*}

8. Solve the problem with a splitting method such as Douglas-Rachford, ADMM, etc. Store the result in a variable xopt_Split. Several variables may converge to the same value, you can pick one arbitrarily and get back to it later.

2.6 Comparison¶

9. Compare the objective values for the points obtained above. Compare the sparsity pattern of these points and compare them to the one of $x_0$.

10. Compare the results obtain when the dimensions and the sparsity increases, eg.
n, m = 100, 200; prop_null = 0.6 or n, m = 400, 200; prop_null = 0.8

3 - a Practical example in Classification¶

Machine Learning as an Optimization problem¶

We have some data $\mathcal{D}$ consisting of $m$ examples; each example consisting of a feature vector $a_i\in\mathbb{R}^n$ and an observation $b_i\in \mathcal{O}$: $\mathcal{D} = \{[a_i,b_i]\}_{i=1..m}$. In this lab, we will consider the student performance dataset.

The goal of supervised learning is to construct a predictor for the observations when given feature vectors.

A popular approach is based on linear models which are based on finding a parameter $x$ such that the real number $\langle a_i , x \rangle$ is used to predict the value of the observation through a predictor function $g:\mathbb{R}\to \mathcal{O}$: $g(\langle a_i , x \rangle)$ is the predicted value from $a_i$.

In order to find such a parameter, we use the available data and a loss $\ell$ that penalizes the error made between the predicted $g(\langle a_i , x \rangle)$ and observed $b_i$ values. For each example $i$, the corresponding error function for a parameter $x$ is $f_i(x) = \ell( g(\langle a_i , x \rangle) ; b_i )$. Using the whole data, the parameter that minimizes the total error is the solution of the minimization problem $$ \min_{x\in\mathbb{R}^n} \frac{1}{m} \sum_{i=1}^m f_i(x) = \frac{1}{m} \sum_{i=1}^m \ell( g(\langle a_i , x \rangle) ; b_i ). $$

Binary Classification with Logisitic Regression¶

In our setup, the observations are binary: $\mathcal{O} = \{-1 , +1 \}$, and the Logistic loss is used to form the following optimization problem \begin{align*} \min_{x\in\mathbb{R}^n } f(x) := \frac{1}{m} \sum_{i=1}^m \log( 1+\exp(-b_i \langle a_i,x \rangle) ) + \frac{\lambda}{2} \|x\|_2^2. \end{align*} where the last term is added as a regularization (of type $\ell_2$, aka Tikhonov) to prevent overfitting.

Under some statistical hypotheses, $x^\star = \arg\min f(x)$ maximizes the likelihood of the labels knowing the features vector. Then, for a new point $d$ with features vector $a$, $$ p_1(a) = \mathbb{P}[a\in \text{ class } +1] = \frac{1}{1+\exp(-\langle a;x^\star \rangle)} $$ Thus, from $a$, if $p_1(a)$ is close to $1$, one can decide that $d$ belongs to class $1$; and the opposite decision if $p(a)$ is close to $0$.

L1-regularization¶

We will consider an $\ell_1$ regularization of this problem to promote sparsity of the iterates. A sparse final solution would select the most important features. The new function (below) is non-smooth but it has a smooth part, $f$, and a non-smooth part, $g$.

\begin{align*} \min_{x\in\mathbb{R}^n } F(x) := \underbrace{ \frac{1}{m} \sum_{i=1}^m \log( 1+\exp(-b_i \langle a_i,x \rangle) ) + \frac{\lambda_2}{2} \|x\|_2^2}_{f(x)} + \underbrace{\lambda_1 \|x\|_1 }_{g(x)}. \end{align*}

Features signification¶

The dataset is comprised of $27$ features described below and the goal is to predict if the student may pass its year or not. It is thus of importance to investigate which features are the most significant for the student success. We will see how the $\ell_1$ regularization can help to this goal.

3.1 Objective definition¶

The size of $x$ is given as $n$, and $L$ represents the Lipschitz constant of the gradient of $f$.

Oracles¶

Below are functions returning the gradient of $f$ and proximity operator of $g$.

Prediction Function¶

Below are two functions that compote the training accuracy (accuracy on the training set), and testing accuracy, given a parameter $x$.

3.2 Proximal Gradient algorithm¶

For minimizing a function $F:\mathbb{R}^n \to \mathbb{R}$ equal to $f+g$ where $f$ is differentiable and the $\mathbf{prox}$ of $g$ is known, given:

• the function to minimize F
• a 1st order oracle for $f$ f_grad
• a proximity operator for $g$ g_prox
• an initialization point x0
• the sought precision PREC
• a maximal number of iterations ITE_MAX
• a display boolean variable PRINT

this algorithm perform iterations of the form $$ x_{k+1} = \mathbf{prox}_{\gamma g}\left( x_k - \gamma \nabla f(x_k) \right) $$ where $\gamma$ is a stepsize to choose.

11. Fill the function below with the proximal gradient algorithm.

12. Run the algorithm and investigate the decrease of the algorithm.

13. Plot the support of the vector $x_k$ (i.e. one point for every non-null coordinate of $x_k$) versus the iteration number. What do you notice? Was it expected?

3.3 Regularization path.¶

We saw above that the algorithm selected some coordinates as the other get to zero. Considering our machine learning task, this translates into the algorithm selecting a subset of the features that will be used for the prediction step.

14. Change the parameter $\lambda_1$ of the problem (lam1) in the code above and investigate how it influences the number of selected features.

In order to quantify the influence of this feature selection, let us consider the regularization path that is the support of the final points obtained by our minimization method versus the value of $\lambda_1$.

15. For $\lambda_1 = 2^{-10},2^{-9}, .. , 2^{1}$, run the proximal gradient algorithm on the obtained problem and store the support of the final point, the prediction performance on the training set (prediction_train) and on the testing set (prediction_test). What is the best regularization choice and which are the most important features of the dataset?

3.4 Bonus questions¶

16. What happens if you replace the $\ell_1$ norm by the $\ell_\infty$ norm $\|x\|_\infty = \max_i |x_i|$?

17. What happens if you replace the proximal gradient by SAGA (see https://arxiv.org/pdf/1407.0202.pdf )