Skip to content
This repository was archived by the owner on Jun 11, 2026. It is now read-only.

PyTorch example for montonic fitting, includes interp1d #746

Open
awf wants to merge 2 commits into
base: master
Choose a base branch
from
Open

PyTorch example for montonic fitting, includes interp1d #746

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
26f64c6
PyTorch interp1d, and montonic fitting
awf May 3, 2021
b8b6ea5
Add weight decay, LBFGS, denser sampling
awf May 5, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
369 changes: 369 additions & 0 deletions users/awf/linterp_and_monotonic.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,369 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"#!conda install -y ipympl\n",
"\n",
"import torch\n",
"import numpy as np\n",
"\n",
"%matplotlib ipympl\n",
"\n",
"import matplotlib.pyplot as plt\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …",
"application/vnd.jupyter.widget-view+json": {
"version_major": 2,
"version_minor": 0,
"model_id": "bb8d6f0e7c5b447c81fca376986d44af"
}
},
"metadata": {}
},
{
"output_type": "execute_result",
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7f6831641220>]"
]
},
"metadata": {},
"execution_count": 2
}
],
"source": [
"\n",
"x_dense = torch.arange(0,3.01,0.1)\n",
"f = lambda x: x ** 2 - 0.2 * x**3\n",
"plt.plot(x_dense, f(x_dense))\n",
"\n",
"n_samples = 10\n",
"samples_x = torch.rand(n_samples) * 2.8 + 0.1\n",
"\n",
"# Sample, but lower\n",
"samples_y = f(samples_x)\n",
"\n",
"plt.plot(samples_x, samples_y, '.')\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …",
"application/vnd.jupyter.widget-view+json": {
"version_major": 2,
"version_minor": 0,
"model_id": "35375136dd694f1db1dff11bffd5f974"
}
},
"metadata": {}
}
],
"source": [
"# Grab from \n",
"def interp1d(x, y, xnew):\n",
" \"\"\"\n",
" PyTorch interp1d, as in scipy.interpolate(x,y,assume_sorted=True)(xnew)\n",
" \n",
" Limited to 1D-1D for now, should not be hard to enhance if needed\n",
" \n",
" github.com/awf\n",
" \"\"\"\n",
" # For each point in x, we want to find index of the knot above it, and hence that below it\n",
" # so knots_x[ind-1] <= x < knots_x[ind]\n",
" inds = torch.bucketize(xnew, x, right=True)\n",
"\n",
" # Call those points xlo, xhi\n",
" xlo = x[inds-1]\n",
" xhi = x[inds]\n",
" ylo = y[inds-1]\n",
" yhi = y[inds]\n",
" \n",
" dx = xhi - xlo\n",
" dy = yhi - ylo\n",
"\n",
" # Then t = (xnew - xlo)/dx\n",
" # ynew = ylo + t * dy\n",
" t = (xnew - xlo) / dx\n",
" return ylo + t * dy\n",
"\n",
"def test_interp1d():\n",
" # scipy example from https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html\n",
" import scipy.interpolate\n",
" import numpy as np\n",
" x = np.arange(0, 10)\n",
" y = np.exp(-x/3.0)\n",
" f = scipy.interpolate.interp1d(x, y)\n",
" xnew = np.arange(0, 9, 0.1)\n",
" ynew = f(xnew) # use interpolation function returned by `interp1d`\n",
"\n",
" # Now do it with torch\n",
" tx = torch.arange(0, 10)\n",
" ty = torch.exp(-tx/3.0)\n",
" txnew = torch.arange(0, 9, 0.1)\n",
" tynew = interp1d(tx,ty,txnew)\n",
" plt.figure()\n",
" plt.plot(x, y, 'ro', xnew, ynew, 'r-',\n",
" tx, ty, 'kx', txnew, tynew, 'k:')\n",
" plt.show()\n",
" \n",
"test_interp1d()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …",
"application/vnd.jupyter.widget-view+json": {
"version_major": 2,
"version_minor": 0,
"model_id": "d7833134ba9a4b5f82281af887dc1700"
}
},
"metadata": {}
},
{
"output_type": "execute_result",
"data": {
"text/plain": [
"SGD (\n",
"Parameter Group 0\n",
" dampening: 0\n",
" lr: 0.01\n",
" momentum: 0.9\n",
" nesterov: False\n",
" weight_decay: 0\n",
")"
]
},
"metadata": {},
"execution_count": 4
}
],
"source": [
"\n",
"\n",
"# Fitting. Define approximator.\n",
"knots_x = torch.linspace(0.0,3.01,6)\n",
"\n",
"# Model is piecewise linear interpolation between knots (\"linear spline\")\n",
"model = lambda weights, x: interp1d(knots_x, weights, x)\n",
"\n",
"# Define losses\n",
"\n",
"# 1. Data term: function should agree with samples\n",
"loss_data = lambda weights: (model(weights, samples_x) - samples_y).abs().sum()\n",
"\n",
"# 2. Monotonicity: pay a penalty wherever weights[i+1] < weights[i]\n",
"loss_monotonic = lambda weights: torch.relu(weights[:-1] - weights[1:]).sum()\n",
"\n",
"loss = lambda weights: loss_data(weights) + loss_monotonic(weights)\n",
"\n",
"\n",
"# Initial estimate of parameters\n",
"torch.manual_seed(41)\n",
"weights = torch.cumsum(torch.randn(knots_x.shape),0)\n",
"\n",
"weights.requires_grad = True\n",
"\n",
"fig = plt.figure()\n",
"with torch.no_grad():\n",
" samples_pred_y = model(weights, samples_x)\n",
" plt.plot(samples_x, samples_y, '.', samples_x, samples_pred_y, 'x')\n",
" y_dense = model(weights, x_dense)\n",
" plt.plot(x_dense, y_dense, '-', knots_x, weights, 'o')\n",
"\n",
"\n",
"# Initialize optimizer\n",
"\n",
"#optimizer = torch.optim.Adam((weights,), lr=0.01, betas=(0.9, 0.999), eps=1e-08, weight_decay=0, amsgrad=False)\n",
"optimizer = torch.optim.SGD((weights,), lr=0.01,momentum=0.9)\n",
"iter_count = 0\n",
"optimizer"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"for k in range(150):\n",
" # Run one iteration\n",
" iter_count += 1\n",
" optimizer.zero_grad()\n",
" loss_value = loss(weights)\n",
" loss_value.backward()\n",
" optimizer.step()\n",
"\n",
" with torch.no_grad():\n",
" fig.clear()\n",
" plt.plot(samples_x, samples_y, 'r.')\n",
" plt.plot(x_dense, model(weights, x_dense), 'b-', knots_x, weights, 'bo')\n",
" plt.title(f'iter {iter_count} loss {loss_value:.2f} {optimizer.__class__}')\n",
" plt.show()\n",
" fig.canvas.draw()\n"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …",
"application/vnd.jupyter.widget-view+json": {
"version_major": 2,
"version_minor": 0,
"model_id": "7fa7b95ef085468f944fa70fbc214a1c"
}
},
"metadata": {}
}
],
"source": [
"# One-sided Fitting. Same model as above, but redefine here anyway\n",
"n_knots = 16\n",
"\n",
"# Fitting. Define approximator.\n",
"knots_x = torch.linspace(0.0,3.01,n_knots)\n",
"\n",
"# Model is piecewise linear interpolation between knots (\"linear spline\")\n",
"model = lambda weights, x: interp1d(knots_x, weights, x)\n",
"\n",
"# Define losses\n",
"\n",
"# 1. Data term: function should not exceed samples: relu(data - model)\n",
"p = 2 # p-norm loss on one-sided term. In principle, we should pay infinity for violating the bound, but quadratic should be enough\n",
"loss_data = lambda weights: (torch.relu(samples_y - model(weights, samples_x)) ** p).sum()\n",
"\n",
"# 2. Monotonicity: pay a penalty wherever weights[i] > weights[i+1]\n",
"loss_monotonic = lambda weights: torch.relu(weights[:-1] - weights[1:]).sum()\n",
"\n",
"weight_decay = 0.1\n",
"\n",
"# 3. Weight decay: \"gravity\" to pull model down\n",
"loss_weight_decay = lambda weights: 0.1 * weights.square().sum()\n",
"\n",
"loss = lambda weights: loss_data(weights) + loss_monotonic(weights) + loss_weight_decay(weights)\n",
"\n",
"# Initial estimate of parameters\n",
"torch.manual_seed(41)\n",
"weights = torch.cumsum(torch.randn(knots_x.shape),0)\n",
"\n",
"# Parameters => requires gradient\n",
"weights.requires_grad = True\n",
"\n",
"# Plot the samples, their model predictions, and the model parameters (i.e. spline knots)\n",
"fig = plt.figure()\n",
"with torch.no_grad():\n",
" samples_pred_y = model(weights, samples_x)\n",
" plt.plot(samples_x, samples_y, '.', samples_x, samples_pred_y, 'x')\n",
" y_dense = model(weights, x_dense)\n",
" plt.plot(x_dense, y_dense, '-', knots_x, weights, 'o')\n",
"fig.canvas.draw()\n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"\n",
"# Initialize optimizer: SGD/LBFGS are the only ones with decent behaviour\n",
"optimizer = torch.optim.LBFGS\n",
"max_iter = 150\n",
"if optimizer == torch.optim.Adam:\n",
" optimizer = optimizer((weights,), lr=0.01, betas=(0.9, 0.999), eps=1e-08, weight_decay=0, amsgrad=False)\n",
"if optimizer == torch.optim.Adagrad:\n",
" optimizer = optimizer((weights,), lr=1.0, lr_decay=0, weight_decay=0, initial_accumulator_value=0, eps=1e-10)\n",
"if optimizer == torch.optim.SGD:\n",
" optimizer = optimizer((weights,), lr=0.01,momentum=0.9)\n",
"if optimizer == torch.optim.LBFGS:\n",
" max_iter = 150\n",
" optimizer = optimizer((weights,), lr=.1, max_iter=20, max_eval=None, \n",
" tolerance_grad=1e-07, tolerance_change=1e-09, \n",
" history_size=4, line_search_fn=None)\n",
"iter_count = 0\n",
"\n",
"weights_best = weights\n",
"loss_best = 1e20\n",
"best_iter = iter_count\n",
"for k in range(max_iter):\n",
" # Run one iteration\n",
" iter_count += 1\n",
" def closure():\n",
" optimizer.zero_grad()\n",
" loss_value = loss(weights)\n",
" loss_value.backward()\n",
" return loss_value\n",
" \n",
" optimizer.step(closure)\n",
"\n",
" with torch.no_grad():\n",
" loss_value = loss(weights)\n",
" if loss_value < loss_best:\n",
" weights_best = weights.clone().detach()\n",
" loss_best = loss_value\n",
" best_iter = iter_count\n",
" fig.clear()\n",
" plt.plot(samples_x, samples_y, 'r.', samples_x, model(weights_best, samples_x), 'rx')\n",
" plt.plot(x_dense, model(weights_best, x_dense), 'b-', knots_x, weights_best, 'bo')\n",
" plt.title(f'iter {iter_count}, loss {loss_best:.3f} at iter {best_iter}, weight_decay={weight_decay}')\n",
" plt.show()\n",
" fig.canvas.draw()\n"
]
}
],
"metadata": {
"kernelspec": {
"name": "python388jvsc74a57bd021ad436eab6efc5a79d7ba41ee8a80a7e2ac2e5a53e8d2a30f44497713f6797f",
"display_name": "Python 3.8.8 64-bit ('ptnightly': conda)"
},
"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.8.8"
},
"metadata": {
"interpreter": {
"hash": "21ad436eab6efc5a79d7ba41ee8a80a7e2ac2e5a53e8d2a30f44497713f6797f"
}
}
},
"nbformat": 4,
"nbformat_minor": 4
}