I am trying to solve the below optimization problem using proximal gradient descent on a dataset in python:

$$f(\theta) = \arg\min_{\theta \in R^d}\frac{1}{m}\sum_{i=1}^m\Big [log(1+exp(x_i\theta))-y_ix_i\theta\Big]+\frac{\lambda_2}{2}||\theta||^2_2 + \lambda_1||\theta||_1$$

where $$x_i \in R^{1 \times d}$$ is sample $$i$$, $$\theta \in R^d$$ and $$y_i \in \{0,1\}$$ is the label for sample $$i$$ and $$\lambda>0$$.

The gradient of $$g(\theta)$$ being

$$\nabla g(\theta) = \frac{1}{m}\sum_{i=1}^m\Big[\frac{x_ie^{x\theta}}{1+e^{x_i\theta}} - x_iy_i\Big]+\theta \lambda_2$$

The dataset contains 784 columns and 2000 datapoints half of which i use for learning $$\theta$$ and the remaining for evaluating accuracy of the classifier. The $$\theta$$ learnt is used to predict labels given by $$\frac{1}{1+exp(-x\theta)}$$. I need to build the classifier using proximal operators, and plot the accuracy per iteration.

This is the code i have come up with but my accuracy is way off, it is expected to converge with >98% accuracy in one to 5 iterations max.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from sklearn.metrics import accuracy_score
np.random.seed(100)
A = data[:,0:784]
y = data[:,-1]
A_train = A[:1000,]
y_train = y[:1000,]
A_test = A[1000:,]
y_test = y[1000:,]
lam_1 = 10
lam_2 = 5
m = len(A_train)

def prox_op(theta,t):
a = np.maximum(theta - t, 0)
b = np.minimum(theta + t, 0)
return a+b

def func(theta):
func_sum = 0
for i in range(m):
func_sum += np.log(1+np.exp(A_train[i]*theta)) - y_train[i]*A_train[i]*theta
term2 = lam_2/2 * np.linalg.norm(theta)
term3 = lam_1 * np.linalg.norm([theta], ord = 1)
return func_sum/m + term2 + term3

for i in range(m):
grad_sum = (A_train[i] * np.exp(A_train[i] * theta)/(1 + np.exp(A_train[i] * theta))) - (A_train[i] * y_train[i])
return grad_sum + (theta * lam_2)

t = 0.01
iterations = 100
theta_new = np.random.uniform(-1,1,784)
theta_list = [theta_new]
func_list = [func(theta_new)]
acc_sc = []

for i in range(iterations-1):
theta_current = theta_list[-1]
theta_new_gd = theta_current - t * grad_g;
theta_new = prox_op(theta_new_gd,lam_1*t)
theta_list.append(theta_new)
func_list.append(func(theta_new))
yy_tr = 1/(1+np.exp(np.dot(-A_train, theta_new)))
acc_sc.append(accuracy_score(yy_tr>t*lam_1,y_train))

plt.plot(func_list)
plt.title('Proximal GD')
plt.xlabel('Iteration')
plt.ylabel('Function Value')
plt.show()



Could some one point where im going wrong please?

This is the plot of function value vs #iterations: plot

• I’d like to see a out of the objective function value versus iteration number. Mar 26, 2022 at 11:32
• @littleO Have added the plot. Mar 26, 2022 at 12:20
• Thanks. That plot looks strange though. There should be only one objective finding value at each iteration. Mar 26, 2022 at 12:36
• Is it because I'm using * and not the numpy.dot function? When i try to substitute the multiplications with numpy.dot, the function values go to infinity though. May be some of the multiplications need to be substituted but im not sure which. Mar 26, 2022 at 13:12

Using slightly different notations, denote $$\phi(\mathbf{w}) = g(\mathbf{w}) + \lambda_1 \|\mathbf{w}\|_1$$ where $$g(\mathbf{w}) =\frac{1}{N} \sum_{n=1}^N \left[ \log \left(1+e_n\right) - y_n \mathbf{x}_n^T \mathbf{w} \right] + \frac{\lambda_2}{2} \|\mathbf{w}\|^2_2$$ and the scalar $$e_n=\exp^{\mathbf{x}_n^T \mathbf{w}}$$.
The update requires the soft-thresholding operator $$\mathbf{w} \leftarrow S_{\lambda_1 t} \left[ \mathbf{w} -t \nabla g(\mathbf{w}) \right]$$ where $$\nabla g(\mathbf{w}) =\frac{1}{N} \sum_{n=1}^N \left[ \left( \frac{e_n}{1+e_n}-y_n \right) \mathbf{x}_n \right] + \lambda_2 \mathbf{w}$$.