<jupyter_start><jupyter_text>A simple example, solved three ways1. CVXPY + MOSEK2. SD ADMM3. Coordinate descent<jupyter_code>%load_ext autoreload
%autoreload 2
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy import signal
from time import time
import seaborn as sns
import cvxpy as cvx
sns.set_style('darkgrid')
import sys
sys.path.append('..')
from osd import Problem
from osd.components import MeanSquareSmall, SmoothSecondDifference, SparseFirstDiffConvex, Boolean
from osd.signal_decomp_bcd import run_bcd
from osd.utilities import progress
from osd.signal_decomp_admm import calc_obj
TOL = 5e-6<jupyter_output><empty_output><jupyter_text>Data generation<jupyter_code>np.random.seed(42)
t = np.linspace(0, 1000, 200)
signal1 = np.sin(2 * np.pi * t * 1 / (500.))
signal2 = signal.square(2 * np.pi * t * 1 / (450.))
X_real = np.zeros((3, len(t)), dtype=float)
X_real[0] = 0.15 * np.random.randn(len(signal1))
X_real[1] = signal1
X_real[2] = signal2
y = np.sum(X_real, axis=0)
K, T = X_real.shape
plt.figure(figsize=(10, 6))
plt.plot(t, np.sum(X_real[1:], axis=0), label='true signal minus noise')
plt.plot(t, y, alpha=0.5, label='observed signal')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>CVXPY + MOSEK<jupyter_code>c1 = MeanSquareSmall(size=T)
c2 = SmoothSecondDifference(weight=1e3/T)
c3 = Boolean(scale=2, shift=-1)
components = [c1, c2, c3]
problem1 = Problem(y, components)
problem1.decompose(how='cvx')
opt_obj_val = problem1.objective_value
opt_obj_val
problem1.plot_decomposition(X_real=X_real);<jupyter_output>No decomposition available.<jupyter_text>SD ADMM<jupyter_code>problem2 = Problem(y, components)
problem2.decompose(how='admm', stopping_tolerance=1e-3)
problem2.objective_value
problem2.decompose(X_init=problem2.estimates, how='bcd', stopping_tolerance=1e-10)
problem2.objective_value
problem2.decompose(X_init=problem2.estimates, how='bcd', stopping_tolerance=1e-10)
problem2.objective_value
plt.figure()
plt.plot(np.r_[problem2.admm_result['obj_vals'], problem2.bcd_result['obj_vals']] - problem2.objective_value + 1e-12)
# plt.axvline(problem2.admm_result['it'], color='red', ls='--')
plt.title('objective value')
plt.yscale('log')
plt.show()
plt.figure()
plt.plot(np.r_[problem2.admm_result['optimality_residual'], problem2.bcd_result['optimality_residual']], label='residual')
plt.axvline(problem2.admm_result['it'], color='red', ls='--')
plt.yscale('log')
plt.legend()
plt.title('internal optimality residual')
plt.show()
problem2.plot_decomposition(X_real=X_real);<jupyter_output><empty_output><jupyter_text>Coordinate Descent<jupyter_code>problem3 = Problem(y, components)
problem3.decompose(how='bcd', stopping_tolerance=TOL)
problem3.objective_value
len(problem3.bcd_result['obj_vals'])
plt.figure()
ref = np.min(np.r_[problem3.bcd_result['obj_vals'], problem2.admm_result['obj_vals']])
plt.plot(problem3.bcd_result['obj_vals'] - ref, label='coordinate descent')
plt.plot(problem2.admm_result['obj_vals'] - ref, label='SD ADMM')
plt.title('suboptimality as compared to CVXPY')
plt.yscale('log')
plt.legend()
plt.show()
plt.figure()
plt.plot(problem3.bcd_result['optimality_residual'], label='coordinate descent')
plt.plot(problem2.admm_result['optimality_residual'], label='SD ADMM')
plt.yscale('log')
plt.title('internal optimality residual')
plt.legend()
plt.show()
plt.scatter(problem3.bcd_result['optimality_residual'], problem3.bcd_result['obj_vals'], 
            label='sd-bcd', marker='.')
plt.scatter(problem2.admm_result['optimality_residual'], problem2.admm_result['obj_vals'], 
            label='sd-admm', marker='.')
plt.xscale('log')
plt.yscale('log')
# plt.xlim(plt.ylim())
plt.xlabel('optimality residual')
plt.ylabel('subpotimality as compared to cvxpy')
# plt.gca().set_aspect('equal')
plt.legend()
plt.title('Comparison of algorithm optimality residual\nto actual difference between objective value and CVXPY value');
problem3.plot_decomposition(X_real=X_real);<jupyter_output><empty_output><jupyter_text>Hybrid<jupyter_code>problem2.objective_value
0.02509511305810197 - 0.02509508363937881<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14932236205227492, 'fraction_numerical': 0.03388189738625363, 'mean_word_length': 7.984782608695652, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '142981', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 1707, 'n_tokens_neox': 1562, 'n_words': 259, 'path': 'notebooks/Simple_nonconvex_example-3_ways.ipynb', 'repo_name': 'bmeyers/optimal-signal-demixing', 'size': '4132'}