{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Block coordinate minimization (BCM) for sparse coding for a wave+spike model\n",
"\n",
"In-class group activity, EECS 598-006 W19\n",
"\n",
"Given a noisy signal $y = x + \\epsilon \\in R^N$,\n",
"this notebook performs *sparse coding* by solving\n",
"$$\n",
"\\hat{z} = \\arg \\min_z \\frac{1}{2} \\| y - D z \\|^2_2 + \\beta \\|z\\|_0\n",
"$$\n",
"to find sparse coefficients $\\hat{z}$\n",
"for a dictionary $D$.\n",
"After finding those coefficients,\n",
"the denoised estimate of the latent signal is\n",
"$ \\hat{x} = D \\hat{z}$.\n",
"\n",
"This notebook considers signals that are modeled\n",
"as the sum of\n",
"a smooth component\n",
"plus a few Kronecker impulses.\n",
"For this \"wave + spike\" model,\n",
"a reasonable dictionary is the\n",
"$N \\times 2N$ matrix\n",
"$D = [D_1 \\ D_2]$\n",
"where\n",
"$D_1$ is an (inverse) discrete cosine transform (DCT),\n",
"and\n",
"$D_2$ is the identity matrix.\n",
"\n",
"There is no analytical solution to the sparse coding problem above.\n",
"However, because $D_1$ and $D_2$ are both unitary matrices,\n",
"it is easy to perform two-block BCM\n",
"that alternates between updating\n",
"coefficient vectors $z_1$ and $z_2$,\n",
"where $z = (z_1,z_2)$.\n",
"This notebook applies BCM to a signal\n",
"formed by adding a (smooth) polynomial and a few spikes.\n",
"\n",
"2019-03-27 Jeff Fessler, University of Michigan"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Class version!!"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Packages needed\n",
"using MIRT: jim\n",
"using FFTW: idct\n",
"using LinearAlgebra: I, norm\n",
"using Plots\n",
"using LaTeXStrings\n",
"using Random: seed!"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"#?idct # uncomment if you want to know more about DCT"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# This 1D problem is small enough that we use the full DCT matrix and identity matrix.\n",
"# For a large-scale problem, we would make D1 and D2 a LinearMap instead.\n",
"N = 2^7\n",
"D1 = idct(Matrix(I,N,N), 1) # inverse DCT basis\n",
"D2 = Matrix(1.0I,N,N) # I\n",
"D = [D1 D2] # over-complete dictionary\n",
"jim(D')"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"