Misc finishing work

cvxgrp · Sep 15, 2024 · d416981 · d416981
1 parent 2d1891c
commit d416981
Show file tree

Hide file tree

Showing 4 changed files with 219 additions and 9 deletions.
diff --git a/examples/cvxpy-example.ipynb b/examples/cvxpy-example.ipynb
@@ -0,0 +1,208 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "55d17ccd-1abc-4c82-a62d-2530d2937c0a",
+   "metadata": {},
+   "outputs": [
+    {
+     "ename": "ModuleNotFoundError",
+     "evalue": "No module named 'randalo'",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
+      "Cell \u001b[0;32mIn[2], line 10\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mcvxpy\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mcp\u001b[39;00m\n\u001b[1;32m      6\u001b[0m \u001b[38;5;66;03m#import torch\u001b[39;00m\n\u001b[1;32m      7\u001b[0m \u001b[38;5;66;03m#from torch.nn import functional as F\u001b[39;00m\n\u001b[1;32m      8\u001b[0m \u001b[38;5;66;03m#from tqdm import tqdm\u001b[39;00m\n\u001b[0;32m---> 10\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mrandalo\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m RandALO\n\u001b[1;32m     12\u001b[0m np\u001b[38;5;241m.\u001b[39mrandom\u001b[38;5;241m.\u001b[39mseed(\u001b[38;5;241m0\u001b[39m)\n\u001b[1;32m     13\u001b[0m \u001b[38;5;66;03m#torch.manual_seed(0)\u001b[39;00m\n",
+      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'randalo'"
+     ]
+    }
+   ],
+   "source": [
+    "import time\n",
+    "\n",
+    "from matplotlib import pyplot as plt\n",
+    "import numpy as np\n",
+    "import cvxpy as cp\n",
+    "#import torch\n",
+    "#from torch.nn import functional as F\n",
+    "#from tqdm import tqdm\n",
+    "\n",
+    "import randalo\n",
+    "\n",
+    "np.random.seed(0)\n",
+    "#torch.manual_seed(0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a2d8471f-feba-4e26-b700-c9b7518bf18d",
+   "metadata": {},
+   "source": [
+    "# Using RandALO with CVXPY\n",
+    "\n",
+    "Since the CVXPY modeling language is significantly more general than the settings where RandALO is appropriate and CVXPYlayers (the differentiation library for CVXPY scales poorly to large problems), we provide a simple modeling language to describe a loss function and a regularizer. We also provide helper methods to transform these loss and regularizers into a CVXPY problem and a linear operator representing the Jacobian."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0761a852-bf1f-4eac-afeb-598e33905015",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "n, p = 100, 30\n",
+    "X = np.random.randn(n, p)\n",
+    "beta = np.zeros(p)\n",
+    "p[0:3] = 1.0\n",
+    "p[9:12] = 2.0\n",
+    "y = X @ beta + 0.1 * np.random.randn(n)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "74eba39c-c246-42b9-bc68-87bfc28df0b0",
+   "metadata": {},
+   "source": [
+    "## Forming the model\n",
+    "\n",
+    "We provide two loss functions: `randalo.LogisticLoss` and `randalo.MLELoss`. Custom loss functions can be written by subclassing `randalo.Loss`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "76bc4a11-8983-49bf-b652-7366476f7e38",
+   "metadata": {},
+   "outputs": [
+    {
+     "ename": "NameError",
+     "evalue": "name 'randalo' is not defined",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
+      "Cell \u001b[0;32mIn[3], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m loss \u001b[38;5;241m=\u001b[39m \u001b[43mrandalo\u001b[49m\u001b[38;5;241m.\u001b[39mMLELoss()\n",
+      "\u001b[0;31mNameError\u001b[0m: name 'randalo' is not defined"
+     ]
+    }
+   ],
+   "source": [
+    "loss = randalo.MLELoss()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "558eb322-9cfd-4745-8cda-4d0fe42d5bb2",
+   "metadata": {},
+   "source": [
+    "Additionally, we support a number of regularizers: `randalo.L2Regularizer`, `randalo.L1Regularizer`, and `SquareRegularizer`. Each of the default regularizers takes an optional matrix (or indices of the variable it should operate on). Custom regularizers may be written by subclassing `Regularizer`. Hyperparameters can be implemented with `randalo.HyperParameter`.\n",
+    "\n",
+    "In this notebook we implement group lasso."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "4512444a-0033-45ec-ba34-c6bada1701c1",
+   "metadata": {},
+   "outputs": [
+    {
+     "ename": "NameError",
+     "evalue": "name 'randalo' is not defined",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
+      "Cell \u001b[0;32mIn[4], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m alpha \u001b[38;5;241m=\u001b[39m \u001b[43mrandalo\u001b[49m\u001b[38;5;241m.\u001b[39mHyperParameter()\n\u001b[1;32m      3\u001b[0m regularizer \u001b[38;5;241m=\u001b[39m alpha \u001b[38;5;241m*\u001b[39m \u001b[38;5;28msum\u001b[39m(randalo\u001b[38;5;241m.\u001b[39mL2Regularizer(\u001b[38;5;28mlist\u001b[39m(\u001b[38;5;28mrange\u001b[39m(i, i \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m3\u001b[39m)) \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(\u001b[38;5;241m0\u001b[39m, p, \u001b[38;5;241m3\u001b[39m)))\n",
+      "\u001b[0;31mNameError\u001b[0m: name 'randalo' is not defined"
+     ]
+    }
+   ],
+   "source": [
+    "alpha = randalo.HyperParameter()\n",
+    "\n",
+    "regularizer = alpha * sum(randalo.L2Regularizer(list(range(i, i + 3)) for i in range(0, p, 3)))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "e578a365-4ea8-448d-afbc-527ab3baf703",
+   "metadata": {},
+   "outputs": [
+    {
+     "ename": "NameError",
+     "evalue": "name 'randalo' is not defined",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
+      "Cell \u001b[0;32mIn[5], line 3\u001b[0m\n\u001b[1;32m      1\u001b[0m b \u001b[38;5;241m=\u001b[39m cp\u001b[38;5;241m.\u001b[39mVariable(p)\n\u001b[0;32m----> 3\u001b[0m prob, J \u001b[38;5;241m=\u001b[39m \u001b[43mrandalo\u001b[49m\u001b[38;5;241m.\u001b[39mgen_cvxpy_jacobian(loss, regularizer, X, b, y) \n",
+      "\u001b[0;31mNameError\u001b[0m: name 'randalo' is not defined"
+     ]
+    }
+   ],
+   "source": [
+    "b = cp.Variable(p)\n",
+    "\n",
+    "prob, J = randalo.gen_cvxpy_jacobian(loss, regularizer, X, b, y)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "77f085f3-d364-483d-8538-7d8c3594eb39",
+   "metadata": {},
+   "outputs": [
+    {
+     "ename": "NameError",
+     "evalue": "name 'alpha' is not defined",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
+      "Cell \u001b[0;32mIn[6], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[43malpha\u001b[49m\u001b[38;5;241m.\u001b[39mvalue \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m0.1\u001b[39m\n\u001b[1;32m      2\u001b[0m prob\u001b[38;5;241m.\u001b[39msolve()\n\u001b[1;32m      4\u001b[0m alo \u001b[38;5;241m=\u001b[39m randalo\u001b[38;5;241m.\u001b[39mRandALO(loss, J, y, b\u001b[38;5;241m.\u001b[39mvalue)\n",
+      "\u001b[0;31mNameError\u001b[0m: name 'alpha' is not defined"
+     ]
+    }
+   ],
+   "source": [
+    "alpha.value = 0.1\n",
+    "prob.solve()\n",
+    "\n",
+    "alo = randalo.RandALO(loss, J, y, b.value)\n",
+    "alo.evaluate(loss)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ccd15851-7844-4455-a58c-d803152d0cad",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "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.12.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/examples/scikit-learn.ipynb b/examples/scikit-learn.ipynb
@@ -366,7 +366,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "alo",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -380,9 +380,9 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.5"
+   "version": "3.12.6"
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
 }
diff --git a/randalo/__init__.py b/randalo/__init__.py
@@ -5,5 +5,5 @@
 from . import utils
 
 from .randalo import RandALO
-from .modeling_layer import SquareRegularizer, L1Regularizer, LogisticLoss, MSELoss
+from .modeling_layer import HyperParameter, Regularizer, SquareRegularizer, L1Regularizer, L2Regularizer, Loss, LogisticLoss, MSELoss
 from .reductions import Jacobian, gen_cvxpy_jacobian
diff --git a/randalo/modeling_layer.py b/randalo/modeling_layer.py
@@ -109,17 +109,18 @@ def to_cvxpy(self, variable):
 
     def get_constraint_hessian_mask(self, beta_hat, epsilon=1e-6):
         scale = self._scale()
-        mask = None
+        if scale == 0.0:
+            return None, None, None
         if self.linear is None:
             return None, torch.diag(
-                    2 * scale * torch.ones_like(beta_hat, dtype=beta_hat.dtype)), mask
+                    2 * scale * torch.ones_like(beta_hat, dtype=beta_hat.dtype)), None
         elif isinstance(linear, list):
-            diag = torch.zeros_like(mask, dtype=beta_hat.dtype)
+            diag = torch.zeros_like(beta_hat)
             diag[linear] = scale
-            return None, torch.diag(diag), mask
+            return None, torch.diag(diag), None
         else:
             A = utils.to_tensor(linear)
-            return None, torch.diag(scale * (A.mT @ A)), mask
+            return None, torch.diag(scale * (A.mT @ A)), None
 
 
 class L1Regularizer(Regularizer):
@@ -166,6 +167,7 @@ def get_constraint_hessian_mask(self, beta_hat, epsilon=1e-6):
             diag = torch.zero_like(beta_hat)
             diag[linear] = 1.0
             return None, self._scale() * (torch.diag(diag) - tilde_b @ tilde_b.T), None
+
         else:
             Lb = linear @ beta_hat
             norm = torch.linalg.norm(Lb)