In [1]:

import os
import sys
from google.colab import drive
drive.mount('/content/drive')
sys.path.append('/content/drive/MyDrive/Projects/2023-HJ-Prox/src/')

from hj_prox import *
from functions import *
import numpy as np
import torch
import torch.nn as nn
import matplotlib.pyplot as plt
import time

torch.set_default_dtype(torch.float64)
torch.manual_seed(0)
device = 'cuda:0'

import os import sys from google.colab import drive drive.mount('/content/drive') sys.path.append('/content/drive/MyDrive/Projects/2023-HJ-Prox/src/') from hj_prox import * from functions import * import numpy as np import torch import torch.nn as nn import matplotlib.pyplot as plt import time torch.set_default_dtype(torch.float64) torch.manual_seed(0) device = 'cuda:0'

Drive already mounted at /content/drive; to attempt to forcibly remount, call drive.mount("/content/drive", force_remount=True).

Approximate Proximals and Moreau Envelopes¶

The Proximal of $f$ is given by \begin{equation} \text{prox}_{tf}(x) \triangleq \text{argmin}_{z\in\mathbb{R}^n} f(z) + \dfrac{1}{2t}\|z-x\|^2, \end{equation} and the Moreau envelope of $f$ is given by \begin{equation} u(x,t) \triangleq \inf_{z\in \mathbb{R}^n} f(z) + \dfrac{1}{2t}\|z-x\|^2 \end{equation}

We leverage the fact that the solution to the Moreau envelope above satisfies the Hamilton-Jacobi Equation \begin{equation} \begin{split} u_t^\delta + \frac{1}{2}\|Du^\delta \|^2 \ = \frac{\delta}{2} \Delta u^\delta \qquad &\text{ in } \mathbb{R}^n\times (0,T] \\ u = f \qquad &\text{ in } \mathbb{R}^n\times \{t = 0\} \end{split} \end{equation} when $\delta = 0$.

By adding a viscous term ($\delta > 0$), we are able to approximate the solution to the HJ equation using the Cole-Hopf formula to obtain \begin{equation} u^\delta(x,t) = - \delta \ln\Big(\Phi_t * \exp(-f/\delta)\Big)(x) = - \delta \ln \int_{\mathbb{R}^n} \Phi(x-y,t) \exp\left(\frac{-f(y)}{\delta}\right) dy \end{equation} where \begin{equation} \Phi(x,t) = \frac{1}{{(4\pi \delta t)}^{n/2}} \exp{\frac{-\|x\|^2}{4\delta t}}. \end{equation} This allows us to write the Moreau Envelope (and the proximal) explicitly as an expectation. In particular, we obtain \begin{equation} \text{prox}_{tf}(x) = \dfrac{\mathbb{E}_{y \sim \mathcal{N}_{x, \delta t I}} \left[y \exp\left(-\delta^{-1}f(y)\right) \right]}{\mathbb{E}_{y \sim \mathcal{N}_{x, \delta t I}} \left[\exp\left(-\delta^{-1}f(y)\right) \right]} \end{equation} and \begin{equation} u(x,t) \approx - \delta \ln \mathbb{E}_{y \sim \mathcal{N}_{x, \delta t I}} \left[y \exp\left(-\delta^{-1}f(y)\right) \right] \end{equation}

In [2]:

# plotting parameters
title_fontsize = 22
fontsize       = 20
fig1 = plt.figure()
my_blue = '#1f77b4'
my_orange = '#F97306'

# plotting parameters title_fontsize = 22 fontsize = 20 fig1 = plt.figure() my_blue = '#1f77b4' my_orange = '#F97306'

<Figure size 640x480 with 0 Axes>

Approximate Prox (Shrink) and Moreau Envelope for L1 Norm¶

Clean L1 Norm Proximal/Moreau Computation¶

In [3]:

f = l1_norm
analytic_prox = l1_norm_prox

# create a grid of x's for plotting
x = torch.linspace(-1,1,100).view(-1,1).to(device)
y_vals = torch.zeros(x.shape, device=device)
prox_true = torch.zeros(x.shape, device=device)
envelope_HJ = torch.zeros(x.shape, device=device)
envelope_true = torch.zeros(x.shape, device=device)
errs = torch.zeros(x.shape, device=device)

t = 2e-1
delta = 1e-2
alpha = 1.0
# y_vals = l1_norm(x) + 0.1*torch.randn(x.shape[0], device=device)
y_vals = f(x)
prox_true = analytic_prox(x, t=t)

prox_HJ = torch.zeros(prox_true.shape, device=device)
n_integral_samples = int(1e5)
for i in range(x.shape[0]):
  temp, ls_iters, temp_envelope = compute_prox(x[i].view(1,1), t, f, int_samples = n_integral_samples, delta = delta, alpha=alpha, device=device)
  prox_HJ[i] = temp
  envelope_HJ[i] = temp_envelope
  envelope_true[i] = f(prox_true[i].view(1,1)) + (1/(2*t))*torch.norm(prox_true[i] - x[i], p=2)**2

  errs[i] = f(temp) + (temp - x[i])/t


# PLOT
fig1 = plt.figure()
plt.style.use('seaborn-whitegrid')
ax = plt.axes()
ax.plot(x.cpu(), y_vals.cpu(), linewidth=3);
ax.plot(x.cpu(), envelope_true.cpu(), linewidth=3, color=my_orange)
ax.plot(x.cpu(), envelope_HJ.cpu(), '--g', linewidth=3)

# ax.set_xlabel("x-axis", fontsize=title_fontsize)
ax.legend(['$f$', 'u', '$u^\delta$'],fontsize=fontsize, loc=9)
ax.tick_params(labelsize=fontsize, which='both', direction='in')

save_str = 'l1_norm_envelope.pdf'
fig1.savefig(save_str, dpi=300 , bbox_inches="tight", pad_inches=0.0)

f = l1_norm analytic_prox = l1_norm_prox # create a grid of x's for plotting x = torch.linspace(-1,1,100).view(-1,1).to(device) y_vals = torch.zeros(x.shape, device=device) prox_true = torch.zeros(x.shape, device=device) envelope_HJ = torch.zeros(x.shape, device=device) envelope_true = torch.zeros(x.shape, device=device) errs = torch.zeros(x.shape, device=device) t = 2e-1 delta = 1e-2 alpha = 1.0 # y_vals = l1_norm(x) + 0.1*torch.randn(x.shape[0], device=device) y_vals = f(x) prox_true = analytic_prox(x, t=t) prox_HJ = torch.zeros(prox_true.shape, device=device) n_integral_samples = int(1e5) for i in range(x.shape[0]): temp, ls_iters, temp_envelope = compute_prox(x[i].view(1,1), t, f, int_samples = n_integral_samples, delta = delta, alpha=alpha, device=device) prox_HJ[i] = temp envelope_HJ[i] = temp_envelope envelope_true[i] = f(prox_true[i].view(1,1)) + (1/(2*t))*torch.norm(prox_true[i] - x[i], p=2)**2 errs[i] = f(temp) + (temp - x[i])/t # PLOT fig1 = plt.figure() plt.style.use('seaborn-whitegrid') ax = plt.axes() ax.plot(x.cpu(), y_vals.cpu(), linewidth=3); ax.plot(x.cpu(), envelope_true.cpu(), linewidth=3, color=my_orange) ax.plot(x.cpu(), envelope_HJ.cpu(), '--g', linewidth=3) # ax.set_xlabel("x-axis", fontsize=title_fontsize) ax.legend(['$f$', 'u', '$u^\delta$'],fontsize=fontsize, loc=9) ax.tick_params(labelsize=fontsize, which='both', direction='in') save_str = 'l1_norm_envelope.pdf' fig1.savefig(save_str, dpi=300 , bbox_inches="tight", pad_inches=0.0)

<ipython-input-3-31371be98f0c>:32: MatplotlibDeprecationWarning: The seaborn styles shipped by Matplotlib are deprecated since 3.6, as they no longer correspond to the styles shipped by seaborn. However, they will remain available as 'seaborn-v0_8-<style>'. Alternatively, directly use the seaborn API instead.
  plt.style.use('seaborn-whitegrid')

In [4]:

#Generate Files for Howard
filename = 'fig2a1_u.dat'
with open(filename, 'w') as csv_file:
  for idx, f_val in enumerate(envelope_true):
    csv_file.write('%0.5e %0.5e\n' % (x[idx], f_val))

filename = 'fig2a1_udelta.dat'
with open(filename, 'w') as csv_file:
  for idx, f_val in enumerate(envelope_HJ):
    csv_file.write('%0.5e %0.5e\n' % (x[idx], f_val))

#Generate Files for Howard filename = 'fig2a1_u.dat' with open(filename, 'w') as csv_file: for idx, f_val in enumerate(envelope_true): csv_file.write('%0.5e %0.5e\n' % (x[idx], f_val)) filename = 'fig2a1_udelta.dat' with open(filename, 'w') as csv_file: for idx, f_val in enumerate(envelope_HJ): csv_file.write('%0.5e %0.5e\n' % (x[idx], f_val))

In [5]:

fig1 = plt.figure()
plt.style.use('seaborn-whitegrid')
ax = plt.axes()

ax.plot(x.cpu(), prox_true.cpu(), linewidth=3, color=my_orange)
ax.plot(x.cpu(), prox_HJ.cpu(), '--', linewidth=3, color='g')

# ax.set_xlabel("x-axis", fontsize=title_fontsize)
ax.legend(['True Prox', 'HJ Prox'],fontsize=fontsize, loc=2)
ax.tick_params(labelsize=fontsize, which='both', direction='in')
save_str = 'l1_norm_prox_comparison.pdf'
fig1.savefig(save_str, dpi=300 , bbox_inches="tight", pad_inches=0.0)

fig1 = plt.figure() plt.style.use('seaborn-whitegrid') ax = plt.axes() ax.plot(x.cpu(), prox_true.cpu(), linewidth=3, color=my_orange) ax.plot(x.cpu(), prox_HJ.cpu(), '--', linewidth=3, color='g') # ax.set_xlabel("x-axis", fontsize=title_fontsize) ax.legend(['True Prox', 'HJ Prox'],fontsize=fontsize, loc=2) ax.tick_params(labelsize=fontsize, which='both', direction='in') save_str = 'l1_norm_prox_comparison.pdf' fig1.savefig(save_str, dpi=300 , bbox_inches="tight", pad_inches=0.0)

<ipython-input-5-c80673add8e2>:2: MatplotlibDeprecationWarning: The seaborn styles shipped by Matplotlib are deprecated since 3.6, as they no longer correspond to the styles shipped by seaborn. However, they will remain available as 'seaborn-v0_8-<style>'. Alternatively, directly use the seaborn API instead.
  plt.style.use('seaborn-whitegrid')

Noisy L1 Norm Proximal/Moreau Computation¶

In [6]:

f = l1_norm_noisy
analytic_prox = l1_norm_prox

# create a grid of x's for plotting
x = torch.linspace(-1,1,100).view(-1,1).to(device)
y_vals = torch.zeros(x.shape, device=device)
y_vals = f(x)
prox_true = torch.zeros(x.shape, device=device)
envelope_HJ = torch.zeros(x.shape, device=device)
envelope_HJ2 = torch.zeros(x.shape, device=device)

t = 0.1
delta = 5e-2
delta2 = 1e-2
alpha = 1.0
prox_true = analytic_prox(x, t=t)

prox_HJ = torch.zeros(prox_true.shape, device=device)
prox_HJ2 = torch.zeros(prox_true.shape, device=device)
n_integral_samples = int(1e4)
for i in range(x.shape[0]):
  temp, ls_iters, temp_envelope = compute_prox(x[i].view(1,1), t, f, int_samples = n_integral_samples, delta = delta, alpha=alpha, device=device)
  prox_HJ[i] = temp
  envelope_HJ[i] = temp_envelope

for i in range(x.shape[0]):
  temp, ls_iters, temp_envelope = compute_prox(x[i].view(1,1), t, f, int_samples = n_integral_samples, delta = delta2, alpha=alpha, device=device)
  prox_HJ2[i] = temp
  envelope_HJ2[i] = temp_envelope

# ------------------------ PLOT without noiseless y ------------------------
fig1 = plt.figure()
plt.style.use('seaborn-whitegrid')
ax = plt.axes()
ax.plot(x.cpu(), y_vals.cpu(), linewidth=3);
ax.plot(x.cpu(), envelope_HJ.cpu(), linewidth=4, color='g')
ax.plot(x.cpu(), envelope_HJ2.cpu(), linewidth=4)
ax.plot(x.cpu(), envelope_l1_norm(x, t=t).cpu(), 'k', linewidth=3) ###### for save_str = 'eye_candy'

# ax.set_xlabel("x-axis", fontsize=title_fontsize)
ax.legend(['noisy $f$', '$u^{\delta_1}$', '$u^{\delta_2}$'],fontsize=fontsize, loc=9)
ax.tick_params(labelsize=fontsize, which='both', direction='in')

# save_str = 'l1_norm_envelope_noisy.pdf'
save_str = 'eye_candy.pdf'
fig1.savefig(save_str, dpi=300 , bbox_inches="tight", pad_inches=0.0)

f = l1_norm_noisy analytic_prox = l1_norm_prox # create a grid of x's for plotting x = torch.linspace(-1,1,100).view(-1,1).to(device) y_vals = torch.zeros(x.shape, device=device) y_vals = f(x) prox_true = torch.zeros(x.shape, device=device) envelope_HJ = torch.zeros(x.shape, device=device) envelope_HJ2 = torch.zeros(x.shape, device=device) t = 0.1 delta = 5e-2 delta2 = 1e-2 alpha = 1.0 prox_true = analytic_prox(x, t=t) prox_HJ = torch.zeros(prox_true.shape, device=device) prox_HJ2 = torch.zeros(prox_true.shape, device=device) n_integral_samples = int(1e4) for i in range(x.shape[0]): temp, ls_iters, temp_envelope = compute_prox(x[i].view(1,1), t, f, int_samples = n_integral_samples, delta = delta, alpha=alpha, device=device) prox_HJ[i] = temp envelope_HJ[i] = temp_envelope for i in range(x.shape[0]): temp, ls_iters, temp_envelope = compute_prox(x[i].view(1,1), t, f, int_samples = n_integral_samples, delta = delta2, alpha=alpha, device=device) prox_HJ2[i] = temp envelope_HJ2[i] = temp_envelope # ------------------------ PLOT without noiseless y ------------------------ fig1 = plt.figure() plt.style.use('seaborn-whitegrid') ax = plt.axes() ax.plot(x.cpu(), y_vals.cpu(), linewidth=3); ax.plot(x.cpu(), envelope_HJ.cpu(), linewidth=4, color='g') ax.plot(x.cpu(), envelope_HJ2.cpu(), linewidth=4) ax.plot(x.cpu(), envelope_l1_norm(x, t=t).cpu(), 'k', linewidth=3) ###### for save_str = 'eye_candy' # ax.set_xlabel("x-axis", fontsize=title_fontsize) ax.legend(['noisy $f$', '$u^{\delta_1}$', '$u^{\delta_2}$'],fontsize=fontsize, loc=9) ax.tick_params(labelsize=fontsize, which='both', direction='in') # save_str = 'l1_norm_envelope_noisy.pdf' save_str = 'eye_candy.pdf' fig1.savefig(save_str, dpi=300 , bbox_inches="tight", pad_inches=0.0)

<ipython-input-6-14052be998c3>:33: MatplotlibDeprecationWarning: The seaborn styles shipped by Matplotlib are deprecated since 3.6, as they no longer correspond to the styles shipped by seaborn. However, they will remain available as 'seaborn-v0_8-<style>'. Alternatively, directly use the seaborn API instead.
  plt.style.use('seaborn-whitegrid')

Approximate Prox and Moreau Envelope for Quadratic Function¶

In [7]:

A = torch.ones(1,1, device=device)
A = A.view(1,1)
b = torch.ones(1, device=device)
def f(x):
  return quadratic(x, A, b)

def analytic_prox(x, t=0.5):
  return quadratic_prox(x, A, b, t=t)

f = f

# create a grid of x's for plotting
x = torch.linspace(-2.0,0.5,100, device=device).view(-1,1)
y_vals = torch.zeros(x.shape, device=device)
envelope_HJ = torch.zeros(x.shape, device=device)
envelope_true = torch.zeros(x.shape, device=device)

t = 5e-1
delta = 5e-2
alpha = 1.0
y_vals = f(x)
errs = torch.zeros(x.shape, device=device)
prox_true = analytic_prox(x, t=t)
prox_HJ = torch.zeros(x.shape, device=device)

n_integral_samples = int(1e5)
for i in range(x.shape[0]):
  temp, ls_iters, temp_envelope = compute_prox(x[i].view(1,1), t, f, int_samples = n_integral_samples, delta = delta, alpha=alpha, device=device)
  prox_HJ[i] = temp
  errs[i] = f(temp) + (temp - x[i])/t
  envelope_HJ[i] = temp_envelope
  envelope_true[i] = f(prox_true[i].view(1,1)) + (1/(2*t))*torch.norm(prox_true[i] - x[i], p=2)**2

# PLOT
fig1 = plt.figure()
plt.style.use('seaborn-whitegrid')
ax = plt.axes()
ax.plot(x.cpu(), y_vals.cpu(), linewidth=3);
ax.plot(x.cpu(), envelope_true.cpu(), linewidth=3, color=my_orange)
ax.plot(x.cpu(), envelope_HJ.cpu(), '--g', linewidth=3)

ax.legend(['$f$', 'u', '$u^\delta$'],fontsize=fontsize, loc=9)
ax.tick_params(labelsize=fontsize, which='both', direction='in')

save_str = 'quadratic_envelope.pdf'
fig1.savefig(save_str, dpi=300 , bbox_inches="tight", pad_inches=0.0)

A = torch.ones(1,1, device=device) A = A.view(1,1) b = torch.ones(1, device=device) def f(x): return quadratic(x, A, b) def analytic_prox(x, t=0.5): return quadratic_prox(x, A, b, t=t) f = f # create a grid of x's for plotting x = torch.linspace(-2.0,0.5,100, device=device).view(-1,1) y_vals = torch.zeros(x.shape, device=device) envelope_HJ = torch.zeros(x.shape, device=device) envelope_true = torch.zeros(x.shape, device=device) t = 5e-1 delta = 5e-2 alpha = 1.0 y_vals = f(x) errs = torch.zeros(x.shape, device=device) prox_true = analytic_prox(x, t=t) prox_HJ = torch.zeros(x.shape, device=device) n_integral_samples = int(1e5) for i in range(x.shape[0]): temp, ls_iters, temp_envelope = compute_prox(x[i].view(1,1), t, f, int_samples = n_integral_samples, delta = delta, alpha=alpha, device=device) prox_HJ[i] = temp errs[i] = f(temp) + (temp - x[i])/t envelope_HJ[i] = temp_envelope envelope_true[i] = f(prox_true[i].view(1,1)) + (1/(2*t))*torch.norm(prox_true[i] - x[i], p=2)**2 # PLOT fig1 = plt.figure() plt.style.use('seaborn-whitegrid') ax = plt.axes() ax.plot(x.cpu(), y_vals.cpu(), linewidth=3); ax.plot(x.cpu(), envelope_true.cpu(), linewidth=3, color=my_orange) ax.plot(x.cpu(), envelope_HJ.cpu(), '--g', linewidth=3) ax.legend(['$f$', 'u', '$u^\delta$'],fontsize=fontsize, loc=9) ax.tick_params(labelsize=fontsize, which='both', direction='in') save_str = 'quadratic_envelope.pdf' fig1.savefig(save_str, dpi=300 , bbox_inches="tight", pad_inches=0.0)

<ipython-input-7-e9418ee90c3a>:36: MatplotlibDeprecationWarning: The seaborn styles shipped by Matplotlib are deprecated since 3.6, as they no longer correspond to the styles shipped by seaborn. However, they will remain available as 'seaborn-v0_8-<style>'. Alternatively, directly use the seaborn API instead.
  plt.style.use('seaborn-whitegrid')

In [8]:

# Figs for Howard
filename = 'fig2a2_udelta.dat'
with open(filename, 'w') as csv_file:
  for idx, f_val in enumerate(envelope_HJ):
    csv_file.write('%0.5e %0.5e\n' % (x[idx], f_val))

# Figs for Howard filename = 'fig2a2_udelta.dat' with open(filename, 'w') as csv_file: for idx, f_val in enumerate(envelope_HJ): csv_file.write('%0.5e %0.5e\n' % (x[idx], f_val))

In [9]:

fig1 = plt.figure()
plt.style.use('seaborn-whitegrid')
ax = plt.axes()

ax.plot(x.cpu(), prox_true.cpu(), linewidth=3, color=my_orange)
ax.plot(x.cpu(), prox_HJ.cpu(), '--', linewidth=3, color='g')

# ax.set_xlabel("x-axis", fontsize=title_fontsize)
ax.legend(['True Prox', 'HJ Prox'],fontsize=fontsize, loc=2)
ax.tick_params(labelsize=fontsize, which='both', direction='in')
save_str = 'quadratic_prox_comparison.pdf'
fig1.savefig(save_str, dpi=300 , bbox_inches="tight", pad_inches=0.0)

fig1 = plt.figure() plt.style.use('seaborn-whitegrid') ax = plt.axes() ax.plot(x.cpu(), prox_true.cpu(), linewidth=3, color=my_orange) ax.plot(x.cpu(), prox_HJ.cpu(), '--', linewidth=3, color='g') # ax.set_xlabel("x-axis", fontsize=title_fontsize) ax.legend(['True Prox', 'HJ Prox'],fontsize=fontsize, loc=2) ax.tick_params(labelsize=fontsize, which='both', direction='in') save_str = 'quadratic_prox_comparison.pdf' fig1.savefig(save_str, dpi=300 , bbox_inches="tight", pad_inches=0.0)

<ipython-input-9-f213e4ce4ac9>:2: MatplotlibDeprecationWarning: The seaborn styles shipped by Matplotlib are deprecated since 3.6, as they no longer correspond to the styles shipped by seaborn. However, they will remain available as 'seaborn-v0_8-<style>'. Alternatively, directly use the seaborn API instead.
  plt.style.use('seaborn-whitegrid')

In [10]:

# Figs for Howard
filename = 'fig2a2_hjprox.dat'
with open(filename, 'w') as csv_file:
  for idx, f_val in enumerate(prox_HJ):
    csv_file.write('%0.5e %0.5e\n' % (x[idx], f_val))

# Figs for Howard filename = 'fig2a2_hjprox.dat' with open(filename, 'w') as csv_file: for idx, f_val in enumerate(prox_HJ): csv_file.write('%0.5e %0.5e\n' % (x[idx], f_val))

Noisy Quadratic Moreau and Prox¶

In [11]:

def f(x):
  return quadratic_noisy(x, A, b)

y_vals = torch.zeros(x.shape, device=device)
envelope_HJ = torch.zeros(x.shape, device=device)
envelope_true = torch.zeros(x.shape, device=device)

t = 5e-1
delta = 5e-2
delta2 = 1e-2
alpha = 1.0
# y_vals = l1_norm(x) + 0.1*torch.randn(x.shape[0], device=device)
y_vals = f(x)
prox_true = analytic_prox(x, t=t)

prox_HJ = torch.zeros(prox_true.shape, device=device)
prox_HJ2 = torch.zeros(prox_true.shape, device=device)
n_integral_samples = int(1e4)
for i in range(x.shape[0]):
  temp, ls_iters, temp_envelope = compute_prox(x[i].view(1,1), t, f, int_samples = n_integral_samples, delta = delta, alpha=alpha, device=device)
  prox_HJ[i] = temp
  envelope_HJ[i] = temp_envelope

for i in range(x.shape[0]):
  temp, ls_iters, temp_envelope = compute_prox(x[i].view(1,1), t, f, int_samples = n_integral_samples, delta = delta2, alpha=alpha, device=device)
  prox_HJ2[i] = temp
  envelope_HJ2[i] = temp_envelope

# PLOT
fig1 = plt.figure()
plt.style.use('seaborn-whitegrid')
ax = plt.axes()
ax.plot(x.cpu(), y_vals.cpu(), linewidth=3);
ax.plot(x.cpu(), envelope_HJ.cpu(), linewidth=4, color='g')
ax.plot(x.cpu(), envelope_HJ2.cpu(), linewidth=4)

# ax.set_xlabel("x-axis", fontsize=title_fontsize)
ax.legend(['noisy $f$', '$u^{\delta_1}$', '$u^{\delta_2}$'],fontsize=fontsize, loc=9)
ax.tick_params(labelsize=fontsize, which='both', direction='in')

save_str = 'quadratic_envelope_noisy.pdf'
fig1.savefig(save_str, dpi=300 , bbox_inches="tight", pad_inches=0.0)

def f(x): return quadratic_noisy(x, A, b) y_vals = torch.zeros(x.shape, device=device) envelope_HJ = torch.zeros(x.shape, device=device) envelope_true = torch.zeros(x.shape, device=device) t = 5e-1 delta = 5e-2 delta2 = 1e-2 alpha = 1.0 # y_vals = l1_norm(x) + 0.1*torch.randn(x.shape[0], device=device) y_vals = f(x) prox_true = analytic_prox(x, t=t) prox_HJ = torch.zeros(prox_true.shape, device=device) prox_HJ2 = torch.zeros(prox_true.shape, device=device) n_integral_samples = int(1e4) for i in range(x.shape[0]): temp, ls_iters, temp_envelope = compute_prox(x[i].view(1,1), t, f, int_samples = n_integral_samples, delta = delta, alpha=alpha, device=device) prox_HJ[i] = temp envelope_HJ[i] = temp_envelope for i in range(x.shape[0]): temp, ls_iters, temp_envelope = compute_prox(x[i].view(1,1), t, f, int_samples = n_integral_samples, delta = delta2, alpha=alpha, device=device) prox_HJ2[i] = temp envelope_HJ2[i] = temp_envelope # PLOT fig1 = plt.figure() plt.style.use('seaborn-whitegrid') ax = plt.axes() ax.plot(x.cpu(), y_vals.cpu(), linewidth=3); ax.plot(x.cpu(), envelope_HJ.cpu(), linewidth=4, color='g') ax.plot(x.cpu(), envelope_HJ2.cpu(), linewidth=4) # ax.set_xlabel("x-axis", fontsize=title_fontsize) ax.legend(['noisy $f$', '$u^{\delta_1}$', '$u^{\delta_2}$'],fontsize=fontsize, loc=9) ax.tick_params(labelsize=fontsize, which='both', direction='in') save_str = 'quadratic_envelope_noisy.pdf' fig1.savefig(save_str, dpi=300 , bbox_inches="tight", pad_inches=0.0)

<ipython-input-11-754e1ba4d0d8>:31: MatplotlibDeprecationWarning: The seaborn styles shipped by Matplotlib are deprecated since 3.6, as they no longer correspond to the styles shipped by seaborn. However, they will remain available as 'seaborn-v0_8-<style>'. Alternatively, directly use the seaborn API instead.
  plt.style.use('seaborn-whitegrid')

In [12]:

# Figs for Howard
filename = 'fig2d2_noisyf.dat'
with open(filename, 'w') as csv_file:
  for idx, f_val in enumerate(y_vals):
    csv_file.write('%0.5e %0.5e\n' % (x[idx], f_val))

# Figs for Howard
filename = 'fig2d2_udelta1.dat'
with open(filename, 'w') as csv_file:
  for idx, f_val in enumerate(envelope_HJ):
    csv_file.write('%0.5e %0.5e\n' % (x[idx], f_val))

# Figs for Howard
filename = 'fig2d2_udelta2.dat'
with open(filename, 'w') as csv_file:
  for idx, f_val in enumerate(envelope_HJ2):
    csv_file.write('%0.5e %0.5e\n' % (x[idx], f_val))

# Figs for Howard filename = 'fig2d2_noisyf.dat' with open(filename, 'w') as csv_file: for idx, f_val in enumerate(y_vals): csv_file.write('%0.5e %0.5e\n' % (x[idx], f_val)) # Figs for Howard filename = 'fig2d2_udelta1.dat' with open(filename, 'w') as csv_file: for idx, f_val in enumerate(envelope_HJ): csv_file.write('%0.5e %0.5e\n' % (x[idx], f_val)) # Figs for Howard filename = 'fig2d2_udelta2.dat' with open(filename, 'w') as csv_file: for idx, f_val in enumerate(envelope_HJ2): csv_file.write('%0.5e %0.5e\n' % (x[idx], f_val))

In [13]:

fig1 = plt.figure()
plt.style.use('seaborn-whitegrid')
ax = plt.axes()

# ax.plot(x.cpu(), prox_true.cpu(), linewidth=4)
ax.plot(x.cpu(), prox_HJ.cpu(), linewidth=4, color='g')

# ax.set_xlabel("x-axis", fontsize=title_fontsize)
ax.legend(['HJ Prox from noisy \n samples'],fontsize=fontsize+2, loc=2)
# ax.tick_params(labelsize=fontsize, which='both', direction='in')
save_str = 'quadratic_prox_noisy.pdf'
fig1.savefig(save_str, dpi=300 , bbox_inches="tight", pad_inches=0.0)

fig1 = plt.figure() plt.style.use('seaborn-whitegrid') ax = plt.axes() # ax.plot(x.cpu(), prox_true.cpu(), linewidth=4) ax.plot(x.cpu(), prox_HJ.cpu(), linewidth=4, color='g') # ax.set_xlabel("x-axis", fontsize=title_fontsize) ax.legend(['HJ Prox from noisy \n samples'],fontsize=fontsize+2, loc=2) # ax.tick_params(labelsize=fontsize, which='both', direction='in') save_str = 'quadratic_prox_noisy.pdf' fig1.savefig(save_str, dpi=300 , bbox_inches="tight", pad_inches=0.0)

<ipython-input-13-dfba4c3aef84>:2: MatplotlibDeprecationWarning: The seaborn styles shipped by Matplotlib are deprecated since 3.6, as they no longer correspond to the styles shipped by seaborn. However, they will remain available as 'seaborn-v0_8-<style>'. Alternatively, directly use the seaborn API instead.
  plt.style.use('seaborn-whitegrid')

Approximate Prox and Moreau Envelope for Log Barrier¶

In [15]:

f = log_barrier
analytic_prox = log_barrier_prox

# create a grid of x's for plotting
x = torch.linspace(2,10,100).view(-1,1).to(device)
y_vals = torch.zeros(x.shape, device=device)
prox_true = torch.zeros(x.shape, device=device)
envelope_HJ = torch.zeros(x.shape, device=device)
envelope_true = torch.zeros(x.shape, device=device)
errs = torch.zeros(x.shape, device=device)

t = 2.0
delta = 1e-1
alpha = 1.0
# y_vals = l1_norm(x) + 0.1*torch.randn(x.shape[0], device=device)
y_vals = f(x)
prox_true = analytic_prox(x, t=t)

prox_HJ = torch.zeros(prox_true.shape, device=device)
n_integral_samples = int(1e4)
for i in range(x.shape[0]):
  temp, ls_iters, temp_envelope = compute_prox(x[i].view(1,1), t, f, int_samples = n_integral_samples, delta = delta, alpha=alpha, device=device)
  prox_HJ[i] = temp
  envelope_HJ[i] = temp_envelope
  envelope_true[i] = f(prox_true[i].view(1,1)) + (1/(2*t))*torch.norm(prox_true[i] - x[i], p=2)**2

  errs[i] = f(temp) + (temp - x[i])/t


# PLOT
fig1 = plt.figure()
plt.style.use('seaborn-whitegrid')
ax = plt.axes()
ax.plot(x.cpu(), y_vals.cpu(), linewidth=3);
ax.plot(x.cpu(), envelope_true.cpu(), linewidth=3, color=my_orange)
ax.plot(x.cpu(), envelope_HJ.cpu(), '--g', linewidth=3)

# ax.set_xlabel("x-axis", fontsize=title_fontsize)
ax.legend(['$f$', 'u', '$u^\delta$'],fontsize=fontsize, loc=9)
ax.tick_params(labelsize=fontsize, which='both', direction='in')

save_str = 'log_barrier_envelope.pdf'
fig1.savefig(save_str, dpi=300 , bbox_inches="tight", pad_inches=0.0)

f = log_barrier analytic_prox = log_barrier_prox # create a grid of x's for plotting x = torch.linspace(2,10,100).view(-1,1).to(device) y_vals = torch.zeros(x.shape, device=device) prox_true = torch.zeros(x.shape, device=device) envelope_HJ = torch.zeros(x.shape, device=device) envelope_true = torch.zeros(x.shape, device=device) errs = torch.zeros(x.shape, device=device) t = 2.0 delta = 1e-1 alpha = 1.0 # y_vals = l1_norm(x) + 0.1*torch.randn(x.shape[0], device=device) y_vals = f(x) prox_true = analytic_prox(x, t=t) prox_HJ = torch.zeros(prox_true.shape, device=device) n_integral_samples = int(1e4) for i in range(x.shape[0]): temp, ls_iters, temp_envelope = compute_prox(x[i].view(1,1), t, f, int_samples = n_integral_samples, delta = delta, alpha=alpha, device=device) prox_HJ[i] = temp envelope_HJ[i] = temp_envelope envelope_true[i] = f(prox_true[i].view(1,1)) + (1/(2*t))*torch.norm(prox_true[i] - x[i], p=2)**2 errs[i] = f(temp) + (temp - x[i])/t # PLOT fig1 = plt.figure() plt.style.use('seaborn-whitegrid') ax = plt.axes() ax.plot(x.cpu(), y_vals.cpu(), linewidth=3); ax.plot(x.cpu(), envelope_true.cpu(), linewidth=3, color=my_orange) ax.plot(x.cpu(), envelope_HJ.cpu(), '--g', linewidth=3) # ax.set_xlabel("x-axis", fontsize=title_fontsize) ax.legend(['$f$', 'u', '$u^\delta$'],fontsize=fontsize, loc=9) ax.tick_params(labelsize=fontsize, whic