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