Generalized LASSO¶
Solve the total variation denoising problem based on Stephen Boyd, et al. “Distributed optimization and statistical learning via the alternating direction method of multipliers”. Now Publishers Inc, 2011.
Problem Description¶
\[\min \frac{1}{2} ||Ax-b||_2^2+\lambda||Fx||_1\]
Here we consider special case when \(A\) is an indentity matrix and \(F\) is the difference matrix, i.e.,
\[\min_x \frac{1}{2}||x-b||_2^2+\lambda\sum_{i=1}^{n-1}|x_{i+1}-x_i|\]
Modules Importing¶
Import all necessary modules and add PyGRANSO src folder to system path.
[1]:
import time
import torch
from pygranso.pygranso import pygranso
from pygranso.pygransoStruct import pygransoStruct
Data Initialization¶
Specify torch device, and generate data
[2]:
device = torch.device('cpu')
n = 80
eta = 0.5 # parameter for penalty term
torch.manual_seed(1)
b = torch.rand(n,1)
pos_one = torch.ones(n-1)
neg_one = -torch.ones(n-1)
F = torch.zeros(n-1,n)
F[:,0:n-1] += torch.diag(neg_one,0)
F[:,1:n] += torch.diag(pos_one,0)
# All the user-provided data (vector/matrix/tensor) must be in torch tensor format.
# As PyTorch tensor is single precision by default, one must explicitly set `dtype=torch.double`.
# Also, please make sure the device of provided torch tensor is the same as opts.torch_device.
F = F.to(device=device, dtype=torch.double) # double precision requireed in torch operations
b = b.to(device=device, dtype=torch.double)
Function Set-Up¶
Encode the optimization variables, and objective and constraint functions.
Note: please strictly follow the format of comb_fn, which will be used in the PyGRANSO main algortihm.
[3]:
# variables and corresponding dimensions.
var_in = {"x": [n,1]}
def user_fn(X_struct,F,b):
x = X_struct.x
# objective function
f = (x-b).T @ (x-b) + eta * torch.norm( F@x, p = 1)
# inequality constraint, matrix form
ci = None
# equality constraint
ce = None
return [f,ci,ce]
comb_fn = lambda X_struct : user_fn(X_struct,F,b)
User Options¶
Specify user-defined options for PyGRANSO
[4]:
opts = pygransoStruct()
opts.torch_device = device
opts.x0 = torch.ones((n,1)).to(device=device, dtype=torch.double)
opts.print_level = 1
opts.print_frequency = 10
Main Algorithm¶
[5]:
start = time.time()
soln = pygranso(var_spec = var_in,combined_fn = comb_fn,user_opts = opts)
end = time.time()
print("Total Wall Time: {}s".format(end - start))
╔═════ QP SOLVER NOTICE ════════════════════════════════════════════════════════════════════════╗
║ PyGRANSO requires a quadratic program (QP) solver that has a quadprog-compatible interface, ║
║ the default is osqp. Users may provide their own wrapper for the QP solver. ║
║ To disable this notice, set opts.quadprog_info_msg = False ║
╚═══════════════════════════════════════════════════════════════════════════════════════════════╝
══════════════════════════════════════════════════════════════════════════════════════════════╗
PyGRANSO: A PyTorch-enabled port of GRANSO with auto-differentiation ║
Version 1.2.0 ║
Licensed under the AGPLv3, Copyright (C) 2021-2022 Tim Mitchell and Buyun Liang ║
══════════════════════════════════════════════════════════════════════════════════════════════╣
Problem specifications: ║
# of variables : 80 ║
# of inequality constraints : 0 ║
# of equality constraints : 0 ║
═════╦════════════╦════════════════╦═════════════╦═══════════════════════╦════════════════════╣
║ Penalty Fn ║ ║ Violation ║ <--- Line Search ---> ║ <- Stationarity -> ║
Iter ║ Mu │ Value ║ Objective ║ Ineq │ Eq ║ SD │ Evals │ t ║ Grads │ Value ║
═════╬════════════╬════════════════╬═════════════╬═══════════════════════╬════════════════════╣
0 ║ - │ - ║ 23.9720755177 ║ - │ - ║ - │ 1 │ 0.000000 ║ 1 │ 9.792257 ║
10 ║ - │ - ║ 6.11242762691 ║ - │ - ║ QN │ 4 │ 0.125000 ║ 1 │ 1.886202 ║
20 ║ - │ - ║ 5.55928361973 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 1 │ 0.992414 ║
30 ║ - │ - ║ 5.31815576231 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 1 │ 1.091197 ║
40 ║ - │ - ║ 5.20063323032 ║ - │ - ║ QN │ 6 │ 0.031250 ║ 1 │ 0.859940 ║
50 ║ - │ - ║ 5.15831247662 ║ - │ - ║ QN │ 6 │ 0.031250 ║ 1 │ 0.953927 ║
60 ║ - │ - ║ 5.13076954691 ║ - │ - ║ QN │ 6 │ 0.031250 ║ 1 │ 0.287711 ║
70 ║ - │ - ║ 5.11261199123 ║ - │ - ║ QN │ 6 │ 0.031250 ║ 1 │ 0.315218 ║
80 ║ - │ - ║ 5.09731135000 ║ - │ - ║ QN │ 3 │ 0.250000 ║ 1 │ 0.135741 ║
90 ║ - │ - ║ 5.08764960138 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 1 │ 0.059596 ║
100 ║ - │ - ║ 5.08307079504 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 1 │ 0.021728 ║
110 ║ - │ - ║ 5.08033270693 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 1 │ 0.019296 ║
120 ║ - │ - ║ 5.07942940005 ║ - │ - ║ QN │ 6 │ 0.031250 ║ 1 │ 0.003682 ║
130 ║ - │ - ║ 5.07875782846 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 1 │ 0.008897 ║
140 ║ - │ - ║ 5.07840647141 ║ - │ - ║ QN │ 6 │ 0.031250 ║ 1 │ 0.002367 ║
150 ║ - │ - ║ 5.07822900856 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 1 │ 0.003504 ║
160 ║ - │ - ║ 5.07812046646 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 2 │ 6.63e-04 ║
170 ║ - │ - ║ 5.07805479861 ║ - │ - ║ QN │ 7 │ 0.015625 ║ 2 │ 0.001014 ║
180 ║ - │ - ║ 5.07801875464 ║ - │ - ║ QN │ 7 │ 0.015625 ║ 3 │ 5.21e-04 ║
190 ║ - │ - ║ 5.07799560661 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 4 │ 3.36e-04 ║
═════╬════════════╬════════════════╬═════════════╬═══════════════════════╬════════════════════╣
║ Penalty Fn ║ ║ Violation ║ <--- Line Search ---> ║ <- Stationarity -> ║
Iter ║ Mu │ Value ║ Objective ║ Ineq │ Eq ║ SD │ Evals │ t ║ Grads │ Value ║
═════╬════════════╬════════════════╬═════════════╬═══════════════════════╬════════════════════╣
200 ║ - │ - ║ 5.07798294167 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 3 │ 3.07e-04 ║
210 ║ - │ - ║ 5.07797328047 ║ - │ - ║ QN │ 4 │ 0.125000 ║ 2 │ 2.42e-04 ║
220 ║ - │ - ║ 5.07796701771 ║ - │ - ║ QN │ 4 │ 0.125000 ║ 2 │ 2.53e-04 ║
230 ║ - │ - ║ 5.07796178749 ║ - │ - ║ QN │ 3 │ 0.250000 ║ 2 │ 2.02e-04 ║
240 ║ - │ - ║ 5.07795849453 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 3 │ 4.33e-04 ║
250 ║ - │ - ║ 5.07795727020 ║ - │ - ║ QN │ 4 │ 0.125000 ║ 3 │ 1.02e-04 ║
260 ║ - │ - ║ 5.07795615914 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 6 │ 8.15e-05 ║
270 ║ - │ - ║ 5.07795585543 ║ - │ - ║ QN │ 4 │ 0.125000 ║ 12 │ 5.60e-05 ║
280 ║ - │ - ║ 5.07795568966 ║ - │ - ║ QN │ 6 │ 0.031250 ║ 15 │ 3.50e-05 ║
290 ║ - │ - ║ 5.07795558861 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 22 │ 1.64e-05 ║
300 ║ - │ - ║ 5.07795552833 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 30 │ 1.32e-05 ║
310 ║ - │ - ║ 5.07795548293 ║ - │ - ║ QN │ 7 │ 0.015625 ║ 36 │ 1.30e-05 ║
320 ║ - │ - ║ 5.07795545267 ║ - │ - ║ QN │ 4 │ 0.125000 ║ 41 │ 1.59e-05 ║
330 ║ - │ - ║ 5.07795539919 ║ - │ - ║ QN │ 3 │ 0.250000 ║ 27 │ 2.61e-05 ║
340 ║ - │ - ║ 5.07795531984 ║ - │ - ║ QN │ 2 │ 0.500000 ║ 7 │ 7.33e-06 ║
350 ║ - │ - ║ 5.07795529483 ║ - │ - ║ QN │ 3 │ 0.250000 ║ 16 │ 1.38e-05 ║
360 ║ - │ - ║ 5.07795527117 ║ - │ - ║ QN │ 3 │ 0.250000 ║ 17 │ 2.24e-05 ║
370 ║ - │ - ║ 5.07795526252 ║ - │ - ║ QN │ 4 │ 0.125000 ║ 19 │ 1.86e-05 ║
380 ║ - │ - ║ 5.07795525902 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 27 │ 1.40e-05 ║
390 ║ - │ - ║ 5.07795525772 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 35 │ 1.04e-05 ║
═════╬════════════╬════════════════╬═════════════╬═══════════════════════╬════════════════════╣
║ Penalty Fn ║ ║ Violation ║ <--- Line Search ---> ║ <- Stationarity -> ║
Iter ║ Mu │ Value ║ Objective ║ Ineq │ Eq ║ SD │ Evals │ t ║ Grads │ Value ║
═════╬════════════╬════════════════╬═════════════╬═══════════════════════╬════════════════════╣
400 ║ - │ - ║ 5.07795525717 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 45 │ 6.25e-06 ║
410 ║ - │ - ║ 5.07795525684 ║ - │ - ║ QN │ 3 │ 0.250000 ║ 54 │ 3.75e-06 ║
420 ║ - │ - ║ 5.07795525661 ║ - │ - ║ QN │ 8 │ 0.007812 ║ 64 │ 3.88e-06 ║
430 ║ - │ - ║ 5.07795525647 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 74 │ 3.47e-06 ║
440 ║ - │ - ║ 5.07795525636 ║ - │ - ║ QN │ 4 │ 0.125000 ║ 84 │ 2.74e-06 ║
450 ║ - │ - ║ 5.07795525628 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 90 │ 2.55e-06 ║
460 ║ - │ - ║ 5.07795525623 ║ - │ - ║ QN │ 6 │ 0.031250 ║ 90 │ 1.50e-06 ║
470 ║ - │ - ║ 5.07795525620 ║ - │ - ║ QN │ 4 │ 0.125000 ║ 90 │ 5.22e-07 ║
480 ║ - │ - ║ 5.07795525619 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 90 │ 3.24e-07 ║
490 ║ - │ - ║ 5.07795525618 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 90 │ 2.32e-07 ║
500 ║ - │ - ║ 5.07795525618 ║ - │ - ║ QN │ 6 │ 0.031250 ║ 90 │ 1.52e-07 ║
510 ║ - │ - ║ 5.07795525617 ║ - │ - ║ QN │ 6 │ 0.031250 ║ 90 │ 9.47e-08 ║
520 ║ - │ - ║ 5.07795525617 ║ - │ - ║ QN │ 5 │ 0.062500 ║ 90 │ 6.51e-08 ║
530 ║ - │ - ║ 5.07795525617 ║ - │ - ║ QN │ 6 │ 0.031250 ║ 90 │ 4.95e-08 ║
540 ║ - │ - ║ 5.07795525617 ║ - │ - ║ QN │ 6 │ 0.031250 ║ 90 │ 4.09e-08 ║
550 ║ - │ - ║ 5.07795525617 ║ - │ - ║ QN │ 4 │ 0.125000 ║ 90 │ 3.47e-08 ║
560 ║ - │ - ║ 5.07795525617 ║ - │ - ║ QN │ 4 │ 0.125000 ║ 90 │ 2.82e-08 ║
570 ║ - │ - ║ 5.07795525617 ║ - │ - ║ QN │ 4 │ 0.125000 ║ 90 │ 1.35e-08 ║
═════╩════════════╩════════════════╩═════════════╩═══════════════════════╩════════════════════╣
F = final iterate, B = Best (to tolerance), MF = Most Feasible ║
Optimization results: ║
═════╦════════════╦════════════════╦═════════════╦═══════════════════════╦════════════════════╣
F ║ │ ║ 5.07795525617 ║ - │ - ║ │ │ ║ │ ║
B ║ │ ║ 5.07795525617 ║ - │ - ║ │ │ ║ │ ║
═════╩════════════╩════════════════╩═════════════╩═══════════════════════╩════════════════════╣
Iterations: 578 ║
Function evaluations: 2802 ║
PyGRANSO termination code: 0 --- converged to stationarity tolerance. ║
══════════════════════════════════════════════════════════════════════════════════════════════╝
Total Wall Time: 7.897493124008179s