sparse deconvolution priors interactive

%23%20%2F%2F%2F%20script%0A%23%20requires-python%20%3D%20%22%3E%3D3.11%22%0A%23%20dependencies%20%3D%20%5B%22marimo%3E%3D0.23.9%22%2C%20%22matplotlib%3D%3D3.10.8%22%2C%20%22numpy%3D%3D2.4.3%22%5D%0A%23%20%0A%23%20%5Btool.marimo.export%5D%0A%23%20lock_kind%20%3D%20%22resolved%22%0A%23%20%2F%2F%2F%0A%0Aimport%20marimo%0A%0A__generated_with%20%3D%20%220.25.1%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22)%0A%0A%0A%40app.cell%0Adef%20imports()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20import%20numpy%20as%20np%0A%0A%20%20%20%20plt.style.use(%22seaborn-v0_8-whitegrid%22)%0A%20%20%20%20return%20mo%2C%20np%2C%20plt%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20intro(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Sparse%20deconvolution%3A%20Gaussian%20versus%20Laplacian%20priors%0A%0A%20%20%20%20An%20inverse%20problem%20runs%20a%20known%20model%20backwards.%20Here%20a%20hidden%20sparse%20signal%20is%20blurred%20and%0A%20%20%20%20contaminated%20by%20noise%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cmathbf%20y%20%3D%20%5Cmathbf%20H%5Cmathbf%20x%20%2B%20%5Cmathbf%20n.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20Direct%20inversion%20magnifies%20noise%20wherever%20the%20blur%20has%20almost%20erased%20a%20frequency.%20A%20prior%0A%20%20%20%20stabilizes%20those%20weak%20directions%E2%80%94but%20its%20**shape**%20decides%20what%20kind%20of%20answer%20is%20preferred.%0A%20%20%20%20This%20notebook%20compares%3A%0A%0A%20%20%20%20-%20a%20Gaussian%20prior%2C%20giving%20quadratic%20Tikhonov%20regularization%20and%20a%20closed-form%20solution%3B%0A%20%20%20%20-%20a%20Laplacian%20prior%2C%20giving%20an%20%24%5Cell_1%24%20penalty%20and%20a%20sparse%20MAP%20estimate%20found%20numerically.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20inverse_helpers(np)%3A%0A%20%20%20%20def%20circular_blur_spectrum(length%2C%20width)%3A%0A%20%20%20%20%20%20%20%20indices%20%3D%20np.arange(length)%0A%20%20%20%20%20%20%20%20distance%20%3D%20np.minimum(indices%2C%20length%20-%20indices)%0A%20%20%20%20%20%20%20%20kernel%20%3D%20np.exp(-0.5%20*%20(distance%20%2F%20width)%20**%202)%0A%20%20%20%20%20%20%20%20kernel%20%3D%20kernel%20%2F%20kernel.sum()%0A%20%20%20%20%20%20%20%20return%20kernel%2C%20np.fft.fft(kernel)%0A%0A%0A%20%20%20%20def%20apply_blur(signal%2C%20spectrum)%3A%0A%20%20%20%20%20%20%20%20return%20np.fft.ifft(np.fft.fft(signal)%20*%20spectrum).real%0A%0A%0A%20%20%20%20def%20apply_adjoint(signal%2C%20spectrum)%3A%0A%20%20%20%20%20%20%20%20return%20np.fft.ifft(np.fft.fft(signal)%20*%20np.conj(spectrum)).real%0A%0A%0A%20%20%20%20def%20soft_threshold(values%2C%20threshold)%3A%0A%20%20%20%20%20%20%20%20return%20np.sign(values)%20*%20np.maximum(np.abs(values)%20-%20threshold%2C%200.0)%0A%0A%0A%20%20%20%20def%20ridge_deconvolution(observation%2C%20spectrum%2C%20regularization)%3A%0A%20%20%20%20%20%20%20%20numerator%20%3D%20np.conj(spectrum)%20*%20np.fft.fft(observation)%0A%20%20%20%20%20%20%20%20denominator%20%3D%20np.abs(spectrum)%20**%202%20%2B%20regularization%0A%20%20%20%20%20%20%20%20return%20np.fft.ifft(numerator%20%2F%20denominator).real%0A%0A%0A%20%20%20%20def%20naive_deconvolution(observation%2C%20spectrum%2C%20floor%3D1e-6)%3A%0A%20%20%20%20%20%20%20%20inverse%20%3D%20np.zeros_like(spectrum%2C%20dtype%3Dcomplex)%0A%20%20%20%20%20%20%20%20stable%20%3D%20np.abs(spectrum)%20%3E%3D%20floor%0A%20%20%20%20%20%20%20%20inverse%5Bstable%5D%20%3D%201.0%20%2F%20spectrum%5Bstable%5D%0A%20%20%20%20%20%20%20%20return%20np.fft.ifft(np.fft.fft(observation)%20*%20inverse).real%0A%0A%0A%20%20%20%20def%20fista_l1_deconvolution(observation%2C%20spectrum%2C%20regularization%2C%20iterations)%3A%0A%20%20%20%20%20%20%20%20estimate%20%3D%20np.zeros_like(observation)%0A%20%20%20%20%20%20%20%20momentum%20%3D%20estimate.copy()%0A%20%20%20%20%20%20%20%20acceleration%20%3D%201.0%0A%20%20%20%20%20%20%20%20lipschitz%20%3D%20max(float(np.max(np.abs(spectrum)%20**%202))%2C%201e-12)%0A%20%20%20%20%20%20%20%20objective%20%3D%20np.empty(iterations)%0A%20%20%20%20%20%20%20%20for%20iteration%20in%20range(iterations)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20gradient%20%3D%20apply_adjoint(apply_blur(momentum%2C%20spectrum)%20-%20observation%2C%20spectrum)%0A%20%20%20%20%20%20%20%20%20%20%20%20next_estimate%20%3D%20soft_threshold(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20momentum%20-%20gradient%20%2F%20lipschitz%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20regularization%20%2F%20lipschitz%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20%20%20%20%20next_acceleration%20%3D%200.5%20*%20(1.0%20%2B%20np.sqrt(1.0%20%2B%204.0%20*%20acceleration**2))%0A%20%20%20%20%20%20%20%20%20%20%20%20momentum%20%3D%20next_estimate%20%2B%20(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20(acceleration%20-%201.0)%20%2F%20next_acceleration%0A%20%20%20%20%20%20%20%20%20%20%20%20)%20*%20(next_estimate%20-%20estimate)%0A%20%20%20%20%20%20%20%20%20%20%20%20estimate%20%3D%20next_estimate%0A%20%20%20%20%20%20%20%20%20%20%20%20acceleration%20%3D%20next_acceleration%0A%20%20%20%20%20%20%20%20%20%20%20%20residual%20%3D%20apply_blur(estimate%2C%20spectrum)%20-%20observation%0A%20%20%20%20%20%20%20%20%20%20%20%20objective%5Biteration%5D%20%3D%20(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%200.5%20*%20float(residual%20%40%20residual)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%2B%20regularization%20*%20float(np.abs(estimate).sum())%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20return%20estimate%2C%20objective%0A%0A%0A%20%20%20%20def%20signal_rmse(reference%2C%20estimate)%3A%0A%20%20%20%20%20%20%20%20return%20float(np.sqrt(np.mean((reference%20-%20estimate)%20**%202)))%0A%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20apply_blur%2C%0A%20%20%20%20%20%20%20%20circular_blur_spectrum%2C%0A%20%20%20%20%20%20%20%20fista_l1_deconvolution%2C%0A%20%20%20%20%20%20%20%20naive_deconvolution%2C%0A%20%20%20%20%20%20%20%20ridge_deconvolution%2C%0A%20%20%20%20%20%20%20%20signal_rmse%2C%0A%20%20%20%20%20%20%20%20soft_threshold%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20forward_explanation(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%201%20%C2%B7%20Build%20the%20forward%20problem%0A%0A%20%20%20%20The%20hidden%20vector%20is%20sparse%20in%20the%20sample%20domain.%20Circular%20convolution%20by%20a%20Gaussian%20kernel%0A%20%20%20%20multiplies%20its%20discrete%20Fourier%20transform%20by%20the%20blur%20response%20%24H%5Bk%5D%24.%20When%20%24%7CH%5Bk%5D%7C%24%20is%20tiny%2C%0A%20%20%20%20direct%20division%20by%20%24H%5Bk%5D%24%20amplifies%20even%20a%20small%20noise%20component%20enormously.%0A%0A%20%20%20%20Adjust%20the%20number%20of%20active%20coefficients%2C%20blur%20width%2C%20and%20noise%20level.%20The%20same%20realization%20is%0A%20%20%20%20then%20used%20by%20both%20priors%20below.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20forward_controls(mo)%3A%0A%20%20%20%20signal_sign_control%20%3D%20mo.ui.dropdown(%0A%20%20%20%20%20%20%20%20options%3D%5B%22Signed%20spikes%22%2C%20%22Positive%20spikes%22%5D%2C%0A%20%20%20%20%20%20%20%20value%3D%22Signed%20spikes%22%2C%0A%20%20%20%20%20%20%20%20label%3D%22Sparse%20signal%22%2C%0A%20%20%20%20%20%20%20%20full_width%3DTrue%2C%0A%20%20%20%20)%0A%20%20%20%20spike_count_control%20%3D%20mo.ui.slider(%0A%20%20%20%20%20%20%20%203%2C%0A%20%20%20%20%20%20%20%2018%2C%0A%20%20%20%20%20%20%20%20value%3D9%2C%0A%20%20%20%20%20%20%20%20step%3D1%2C%0A%20%20%20%20%20%20%20%20label%3D%22Number%20of%20nonzero%20coefficients%22%2C%0A%20%20%20%20%20%20%20%20show_value%3DTrue%2C%0A%20%20%20%20%20%20%20%20full_width%3DTrue%2C%0A%20%20%20%20)%0A%20%20%20%20blur_width_control%20%3D%20mo.ui.slider(%0A%20%20%20%20%20%20%20%200.8%2C%0A%20%20%20%20%20%20%20%206.0%2C%0A%20%20%20%20%20%20%20%20value%3D3.0%2C%0A%20%20%20%20%20%20%20%20step%3D0.2%2C%0A%20%20%20%20%20%20%20%20label%3D%22Gaussian%20blur%20width%22%2C%0A%20%20%20%20%20%20%20%20show_value%3DTrue%2C%0A%20%20%20%20%20%20%20%20full_width%3DTrue%2C%0A%20%20%20%20)%0A%20%20%20%20noise_std_control%20%3D%20mo.ui.slider(%0A%20%20%20%20%20%20%20%200.0%2C%0A%20%20%20%20%20%20%20%200.10%2C%0A%20%20%20%20%20%20%20%20value%3D0.025%2C%0A%20%20%20%20%20%20%20%20step%3D0.005%2C%0A%20%20%20%20%20%20%20%20label%3D%22Noise%20standard%20deviation%22%2C%0A%20%20%20%20%20%20%20%20show_value%3DTrue%2C%0A%20%20%20%20%20%20%20%20full_width%3DTrue%2C%0A%20%20%20%20)%0A%20%20%20%20forward_seed_control%20%3D%20mo.ui.slider(%0A%20%20%20%20%20%20%20%200%2C%0A%20%20%20%20%20%20%20%2099%2C%0A%20%20%20%20%20%20%20%20value%3D3%2C%0A%20%20%20%20%20%20%20%20step%3D1%2C%0A%20%20%20%20%20%20%20%20label%3D%22Random%20seed%22%2C%0A%20%20%20%20%20%20%20%20show_value%3DTrue%2C%0A%20%20%20%20%20%20%20%20full_width%3DTrue%2C%0A%20%20%20%20)%0A%20%20%20%20mo.vstack(%5B%0A%20%20%20%20%20%20%20%20signal_sign_control%2C%0A%20%20%20%20%20%20%20%20mo.hstack(%5Bspike_count_control%2C%20blur_width_control%5D)%2C%0A%20%20%20%20%20%20%20%20mo.hstack(%5Bnoise_std_control%2C%20forward_seed_control%5D)%2C%0A%20%20%20%20%5D)%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20blur_width_control%2C%0A%20%20%20%20%20%20%20%20forward_seed_control%2C%0A%20%20%20%20%20%20%20%20noise_std_control%2C%0A%20%20%20%20%20%20%20%20signal_sign_control%2C%0A%20%20%20%20%20%20%20%20spike_count_control%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20forward_plot(%0A%20%20%20%20apply_blur%2C%0A%20%20%20%20blur_width_control%2C%0A%20%20%20%20circular_blur_spectrum%2C%0A%20%20%20%20forward_seed_control%2C%0A%20%20%20%20noise_std_control%2C%0A%20%20%20%20np%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20signal_sign_control%2C%0A%20%20%20%20spike_count_control%2C%0A)%3A%0A%20%20%20%20signal_length%20%3D%20192%0A%20%20%20%20signal_index%20%3D%20np.arange(signal_length)%0A%20%20%20%20_forward_rng%20%3D%20np.random.default_rng(int(forward_seed_control.value))%0A%20%20%20%20spike_positions%20%3D%20np.sort(%0A%20%20%20%20%20%20%20%20_forward_rng.choice(%0A%20%20%20%20%20%20%20%20%20%20%20%20np.arange(8%2C%20signal_length%20-%208)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20size%3Dint(spike_count_control.value)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20replace%3DFalse%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20)%0A%20%20%20%20_spike_amplitudes%20%3D%20_forward_rng.uniform(0.55%2C%201.15%2C%20size%3Dlen(spike_positions))%0A%20%20%20%20if%20signal_sign_control.value%20%3D%3D%20%22Signed%20spikes%22%3A%0A%20%20%20%20%20%20%20%20_spike_amplitudes%20%3D%20_spike_amplitudes%20*%20_forward_rng.choice(%0A%20%20%20%20%20%20%20%20%20%20%20%20np.array(%5B-1.0%2C%201.0%5D)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20size%3Dlen(spike_positions)%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20true_signal%20%3D%20np.zeros(signal_length)%0A%20%20%20%20true_signal%5Bspike_positions%5D%20%3D%20_spike_amplitudes%0A%20%20%20%20blur_kernel%2C%20blur_spectrum%20%3D%20circular_blur_spectrum(%0A%20%20%20%20%20%20%20%20signal_length%2C%0A%20%20%20%20%20%20%20%20float(blur_width_control.value)%2C%0A%20%20%20%20)%0A%20%20%20%20clean_observation%20%3D%20apply_blur(true_signal%2C%20blur_spectrum)%0A%20%20%20%20noise_standard_deviation%20%3D%20float(noise_std_control.value)%0A%20%20%20%20observed_signal%20%3D%20clean_observation%20%2B%20_forward_rng.normal(%0A%20%20%20%20%20%20%20%200.0%2C%0A%20%20%20%20%20%20%20%20noise_standard_deviation%2C%0A%20%20%20%20%20%20%20%20signal_length%2C%0A%20%20%20%20)%0A%20%20%20%20blur_gain%20%3D%20np.abs(blur_spectrum)%0A%20%20%20%20weak_frequency_fraction%20%3D%20float(np.mean(blur_gain%20%3C%200.05))%0A%20%20%20%20_blur_frequency%20%3D%20np.fft.fftshift(np.fft.fftfreq(signal_length))%0A%20%20%20%20_blur_gain_shifted%20%3D%20np.fft.fftshift(blur_gain)%0A%20%20%20%20fig_forward%2C%20_forward_axes%20%3D%20plt.subplots(1%2C%203%2C%20figsize%3D(12.0%2C%203.55))%0A%20%20%20%20_forward_axes%5B0%5D.axhline(0%2C%20color%3D%22%23777777%22%2C%20lw%3D0.8)%0A%20%20%20%20_forward_axes%5B0%5D.vlines(signal_index%2C%200%2C%20true_signal%2C%20color%3D%22%23285f8f%22%2C%20lw%3D1.8)%0A%20%20%20%20_forward_axes%5B0%5D.scatter(spike_positions%2C%20true_signal%5Bspike_positions%5D%2C%20color%3D%22%23285f8f%22%2C%20s%3D25%2C%20zorder%3D3)%0A%20%20%20%20_forward_axes%5B0%5D.set(%0A%20%20%20%20%20%20%20%20title%3Df%22Hidden%20signal%3A%20%7Blen(spike_positions)%7D%20active%20samples%22%2C%0A%20%20%20%20%20%20%20%20xlabel%3D%22sample%20index%22%2C%0A%20%20%20%20%20%20%20%20ylabel%3D%22amplitude%22%2C%0A%20%20%20%20)%0A%20%20%20%20_forward_axes%5B1%5D.plot(signal_index%2C%20clean_observation%2C%20color%3D%22%23f28e2b%22%2C%20lw%3D2.0%2C%20label%3D%22blurred%22)%0A%20%20%20%20_forward_axes%5B1%5D.plot(signal_index%2C%20observed_signal%2C%20color%3D%22%234c78a8%22%2C%20lw%3D1.1%2C%20alpha%3D0.78%2C%20label%3D%22blurred%20%2B%20noise%22)%0A%20%20%20%20_forward_axes%5B1%5D.set(title%3D%22What%20the%20sensor%20observes%22%2C%20xlabel%3D%22sample%20index%22%2C%20ylabel%3D%22amplitude%22)%0A%20%20%20%20_forward_axes%5B1%5D.legend(frameon%3DTrue%2C%20fontsize%3D8)%0A%20%20%20%20_forward_axes%5B2%5D.semilogy(_blur_frequency%2C%20np.maximum(_blur_gain_shifted%2C%201e-12)%2C%20color%3D%22%237b4fa3%22%2C%20lw%3D2.1)%0A%20%20%20%20_forward_axes%5B2%5D.axhline(0.05%2C%20color%3D%22%23d1495b%22%2C%20ls%3D%22--%22%2C%20lw%3D1.4%2C%20label%3D%225%25%20gain%22)%0A%20%20%20%20_forward_axes%5B2%5D.set(%0A%20%20%20%20%20%20%20%20title%3Df%22Blur%20suppresses%20%7B100%20*%20weak_frequency_fraction%3A.0f%7D%25%20of%20frequencies%20below%205%25%22%2C%0A%20%20%20%20%20%20%20%20xlabel%3D%22normalized%20frequency%22%2C%0A%20%20%20%20%20%20%20%20ylabel%3Dr%22gain%20%24%7CH%5Bk%5D%7C%24%22%2C%0A%20%20%20%20%20%20%20%20ylim%3D(1e-8%2C%201.5)%2C%0A%20%20%20%20)%0A%20%20%20%20_forward_axes%5B2%5D.legend(frameon%3DTrue%2C%20fontsize%3D8)%0A%20%20%20%20for%20_forward_axis%20in%20_forward_axes%3A%0A%20%20%20%20%20%20%20%20_forward_axis.grid(alpha%3D0.24)%0A%20%20%20%20fig_forward.suptitle(%22Forward%20model%3A%20sparse%20cause%20%E2%86%92%20blurred%2C%20noisy%20effect%22%2C%20fontweight%3D%22bold%22%2C%20y%3D1.03)%0A%20%20%20%20fig_forward.tight_layout()%0A%20%20%20%20fig_forward%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20blur_spectrum%2C%0A%20%20%20%20%20%20%20%20observed_signal%2C%0A%20%20%20%20%20%20%20%20signal_index%2C%0A%20%20%20%20%20%20%20%20signal_length%2C%0A%20%20%20%20%20%20%20%20spike_positions%2C%0A%20%20%20%20%20%20%20%20true_signal%2C%0A%20%20%20%20%20%20%20%20weak_frequency_fraction%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20forward_summary(mo%2C%20weak_frequency_fraction)%3A%0A%20%20%20%20mo.callout(%0A%20%20%20%20%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22The%20blur%20leaves%20**%7B100%20*%20weak_frequency_fraction%3A.1f%7D%25**%20of%20Fourier%20directions%20with%20less%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%22than%205%25%20gain.%20Those%20directions%20carry%20little%20information%20about%20the%20hidden%20signal%2C%20so%20an%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%22inverse%20must%20decide%20whether%20to%20amplify%20them%20or%20trust%20a%20prior.%22%0A%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20kind%3D%22warn%22%20if%20weak_frequency_fraction%20%3E%200.25%20else%20%22info%22%2C%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20prior_explanation(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%202%20%C2%B7%20Prior%20shape%20becomes%20an%20optimization%20penalty%0A%0A%20%20%20%20With%20Gaussian%20noise%2C%20MAP%20estimation%20minimizes%20a%20data-misfit%20term%20plus%20the%20negative%20log-prior.%0A%20%20%20%20Ignoring%20constants%2C%0A%0A%20%20%20%20%24%24%0A%20%20%20%20-%5Clog%20p_%7B%5Cmathrm%20G%7D(x)%20%5Cpropto%20x%5E2%2C%0A%20%20%20%20%5Cqquad%0A%20%20%20%20-%5Clog%20p_%7B%5Cmathrm%20L%7D(x)%20%5Cpropto%20%7Cx%7C.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20The%20quadratic%20Gaussian%20penalty%20is%20smooth%20and%20shrinks%20every%20coefficient.%20The%20Laplacian%20penalty%0A%20%20%20%20has%20a%20cusp%20at%20zero%3B%20its%20proximal%20update%20is%20soft%20thresholding%2C%20which%20sends%20a%20whole%20interval%20of%0A%20%20%20%20small%20coefficients%20to%20**exactly%20zero**.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20prior_controls(mo)%3A%0A%20%20%20%20shrinkage_strength_control%20%3D%20mo.ui.slider(%0A%20%20%20%20%20%20%20%200.0%2C%0A%20%20%20%20%20%20%20%201.5%2C%0A%20%20%20%20%20%20%20%20value%3D0.8%2C%0A%20%20%20%20%20%20%20%20step%3D0.05%2C%0A%20%20%20%20%20%20%20%20label%3D%22Illustrative%20shrinkage%20strength%22%2C%0A%20%20%20%20%20%20%20%20show_value%3DTrue%2C%0A%20%20%20%20%20%20%20%20full_width%3DTrue%2C%0A%20%20%20%20)%0A%20%20%20%20shrinkage_strength_control%0A%20%20%20%20return%20(shrinkage_strength_control%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20prior_plot(np%2C%20plt%2C%20shrinkage_strength_control%2C%20soft_threshold)%3A%0A%20%20%20%20_prior_x%20%3D%20np.linspace(-3.5%2C%203.5%2C%20801)%0A%20%20%20%20_prior_gaussian%20%3D%20np.exp(-0.5%20*%20_prior_x**2)%20%2F%20np.sqrt(2%20*%20np.pi)%0A%20%20%20%20_prior_laplace%20%3D%20np.exp(-np.sqrt(2)%20*%20np.abs(_prior_x))%20%2F%20np.sqrt(2)%0A%20%20%20%20_prior_gaussian_penalty%20%3D%200.5%20*%20_prior_x**2%0A%20%20%20%20_prior_laplace_penalty%20%3D%20np.sqrt(2)%20*%20np.abs(_prior_x)%0A%20%20%20%20_shrinkage_strength%20%3D%20float(shrinkage_strength_control.value)%0A%20%20%20%20_shrinkage_input%20%3D%20np.linspace(-3.0%2C%203.0%2C%20601)%0A%20%20%20%20_gaussian_shrinkage%20%3D%20_shrinkage_input%20%2F%20(1.0%20%2B%20_shrinkage_strength)%0A%20%20%20%20_laplacian_shrinkage%20%3D%20soft_threshold(_shrinkage_input%2C%20_shrinkage_strength)%0A%20%20%20%20fig_priors%2C%20_prior_axes%20%3D%20plt.subplots(1%2C%203%2C%20figsize%3D(12.0%2C%203.55))%0A%20%20%20%20_prior_axes%5B0%5D.plot(_prior_x%2C%20_prior_gaussian%2C%20color%3D%22%23f28e2b%22%2C%20lw%3D2.2%2C%20label%3D%22Gaussian%22)%0A%20%20%20%20_prior_axes%5B0%5D.plot(_prior_x%2C%20_prior_laplace%2C%20color%3D%22%23285f8f%22%2C%20lw%3D2.2%2C%20label%3D%22Laplacian%22)%0A%20%20%20%20_prior_axes%5B0%5D.fill_between(_prior_x%2C%20_prior_laplace%2C%20color%3D%22%237ba7cc%22%2C%20alpha%3D0.18)%0A%20%20%20%20_prior_axes%5B0%5D.set(title%3D%22Equal-variance%20priors%22%2C%20xlabel%3Dr%22coefficient%20%24x%24%22%2C%20ylabel%3D%22density%22)%0A%20%20%20%20_prior_axes%5B0%5D.legend(frameon%3DTrue%2C%20fontsize%3D9)%0A%20%20%20%20_prior_axes%5B1%5D.plot(_prior_x%2C%20_prior_gaussian_penalty%2C%20color%3D%22%23f28e2b%22%2C%20lw%3D2.2%2C%20label%3Dr%22quadratic%20%24x%5E2%2F2%24%22)%0A%20%20%20%20_prior_axes%5B1%5D.plot(_prior_x%2C%20_prior_laplace_penalty%2C%20color%3D%22%23285f8f%22%2C%20lw%3D2.2%2C%20label%3Dr%22absolute%20%24%5Csqrt%7B2%7D%7Cx%7C%24%22)%0A%20%20%20%20_prior_axes%5B1%5D.set(title%3D%22Negative%20log-priors%22%2C%20xlabel%3Dr%22coefficient%20%24x%24%22%2C%20ylabel%3D%22penalty%22)%0A%20%20%20%20_prior_axes%5B1%5D.legend(frameon%3DTrue%2C%20fontsize%3D9)%0A%20%20%20%20_prior_axes%5B2%5D.plot(_shrinkage_input%2C%20_shrinkage_input%2C%20color%3D%22%23888888%22%2C%20ls%3D%22%3A%22%2C%20lw%3D1.3%2C%20label%3D%22no%20shrinkage%22)%0A%20%20%20%20_prior_axes%5B2%5D.plot(_shrinkage_input%2C%20_gaussian_shrinkage%2C%20color%3D%22%23f28e2b%22%2C%20lw%3D2.2%2C%20label%3D%22Gaussian%20%2F%20ridge%22)%0A%20%20%20%20_prior_axes%5B2%5D.plot(_shrinkage_input%2C%20_laplacian_shrinkage%2C%20color%3D%22%23285f8f%22%2C%20lw%3D2.2%2C%20label%3D%22Laplacian%20%2F%20soft%20threshold%22)%0A%20%20%20%20_prior_axes%5B2%5D.axvspan(-_shrinkage_strength%2C%20_shrinkage_strength%2C%20color%3D%22%2359a14f%22%2C%20alpha%3D0.13%2C%20label%3D%22mapped%20to%20zero%22)%0A%20%20%20%20_prior_axes%5B2%5D.set(%0A%20%20%20%20%20%20%20%20title%3Df%22Coefficient%20update%20at%20strength%20%7B_shrinkage_strength%3A.2f%7D%22%2C%0A%20%20%20%20%20%20%20%20xlabel%3D%22input%20coefficient%22%2C%0A%20%20%20%20%20%20%20%20ylabel%3D%22updated%20coefficient%22%2C%0A%20%20%20%20)%0A%20%20%20%20_prior_axes%5B2%5D.legend(frameon%3DTrue%2C%20fontsize%3D8)%0A%20%20%20%20for%20_prior_axis%20in%20_prior_axes%3A%0A%20%20%20%20%20%20%20%20_prior_axis.grid(alpha%3D0.24)%0A%20%20%20%20fig_priors.suptitle(%22Same%20data%20term%2C%20different%20prior%20geometry%22%2C%20fontweight%3D%22bold%22%2C%20y%3D1.03)%0A%20%20%20%20fig_priors.tight_layout()%0A%20%20%20%20fig_priors%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20prior_summary(mo%2C%20shrinkage_strength_control)%3A%0A%20%20%20%20_threshold%20%3D%20float(shrinkage_strength_control.value)%0A%20%20%20%20mo.callout(%0A%20%20%20%20%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20%20%20%20%20%22The%20Gaussian%20update%20is%20a%20straight%20line%3A%20every%20nonzero%20input%20stays%20nonzero.%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22The%20Laplacian%20update%20maps%20the%20entire%20interval%20**%5B%7B-_threshold%3A.2f%7D%2C%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22%7B_threshold%3A.2f%7D%5D**%20to%20zero.%20That%20threshold%20is%20the%20mechanism%20that%20can%20recover%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%22a%20sparse%20support%20rather%20than%20merely%20a%20small-amplitude%20dense%20signal.%22%0A%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20kind%3D%22info%22%2C%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20reconstruction_explanation(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%203%20%C2%B7%20Solve%20the%20same%20inverse%20problem%20with%20two%20priors%0A%0A%20%20%20%20The%20Gaussian-prior%20estimate%20solves%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cwidehat%7B%5Cmathbf%20x%7D_%7B%5Cmathrm%20G%7D%0A%20%20%20%20%3D%20%5Carg%5Cmin_%7B%5Cmathbf%20x%7D%0A%20%20%20%20%5Cfrac%7B1%7D%7B2%7D%5C%7C%5Cmathbf%20H%5Cmathbf%20x-%5Cmathbf%20y%5C%7C_2%5E2%0A%20%20%20%20%2B%20%5Cfrac%7B%5Clambda_2%7D%7B2%7D%5C%7C%5Cmathbf%20x%5C%7C_2%5E2%2C%0A%20%20%20%20%24%24%0A%0A%20%20%20%20with%20a%20closed-form%20Fourier-domain%20solution.%20The%20Laplacian-prior%20MAP%20estimate%20replaces%20the%0A%20%20%20%20quadratic%20penalty%20by%20%24%5Clambda_1%5C%7C%5Cmathbf%20x%5C%7C_1%24.%20There%20is%20no%20linear%20closed%20form%2C%20so%20the%20notebook%0A%20%20%20%20uses%20FISTA%3A%20gradient%20steps%20on%20the%20data%20term%20alternated%20with%20soft%20thresholding.%0A%0A%20%20%20%20The%20two%20regularization%20parameters%20are%20intentionally%20independent.%20Compare%20solutions%20at%20equal%0A%20%20%20%20values%20first%2C%20then%20move%20each%20to%20see%20its%20own%20fit-versus-prior%20trade-off.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20reconstruction_controls(mo)%3A%0A%20%20%20%20gaussian_log_lambda_control%20%3D%20mo.ui.slider(%0A%20%20%20%20%20%20%20%20-4.0%2C%0A%20%20%20%20%20%20%20%20-0.5%2C%0A%20%20%20%20%20%20%20%20value%3D-1.7%2C%0A%20%20%20%20%20%20%20%20step%3D0.1%2C%0A%20%20%20%20%20%20%20%20label%3D%22Gaussian%20prior%3A%20log10%20%CE%BB%E2%82%82%22%2C%0A%20%20%20%20%20%20%20%20show_value%3DTrue%2C%0A%20%20%20%20%20%20%20%20full_width%3DTrue%2C%0A%20%20%20%20)%0A%20%20%20%20laplacian_log_lambda_control%20%3D%20mo.ui.slider(%0A%20%20%20%20%20%20%20%20-4.0%2C%0A%20%20%20%20%20%20%20%20-0.5%2C%0A%20%20%20%20%20%20%20%20value%3D-1.7%2C%0A%20%20%20%20%20%20%20%20step%3D0.1%2C%0A%20%20%20%20%20%20%20%20label%3D%22Laplacian%20prior%3A%20log10%20%CE%BB%E2%82%81%22%2C%0A%20%20%20%20%20%20%20%20show_value%3DTrue%2C%0A%20%20%20%20%20%20%20%20full_width%3DTrue%2C%0A%20%20%20%20)%0A%20%20%20%20fista_iterations_control%20%3D%20mo.ui.slider(%0A%20%20%20%20%20%20%20%2050%2C%0A%20%20%20%20%20%20%20%20600%2C%0A%20%20%20%20%20%20%20%20value%3D350%2C%0A%20%20%20%20%20%20%20%20step%3D25%2C%0A%20%20%20%20%20%20%20%20label%3D%22FISTA%20iterations%22%2C%0A%20%20%20%20%20%20%20%20show_value%3DTrue%2C%0A%20%20%20%20%20%20%20%20full_width%3DTrue%2C%0A%20%20%20%20)%0A%20%20%20%20mo.vstack(%5B%0A%20%20%20%20%20%20%20%20mo.hstack(%5Bgaussian_log_lambda_control%2C%20laplacian_log_lambda_control%5D)%2C%0A%20%20%20%20%20%20%20%20fista_iterations_control%2C%0A%20%20%20%20%5D)%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20fista_iterations_control%2C%0A%20%20%20%20%20%20%20%20gaussian_log_lambda_control%2C%0A%20%20%20%20%20%20%20%20laplacian_log_lambda_control%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20reconstruction_plot(%0A%20%20%20%20apply_blur%2C%0A%20%20%20%20blur_spectrum%2C%0A%20%20%20%20fista_iterations_control%2C%0A%20%20%20%20fista_l1_deconvolution%2C%0A%20%20%20%20gaussian_log_lambda_control%2C%0A%20%20%20%20laplacian_log_lambda_control%2C%0A%20%20%20%20naive_deconvolution%2C%0A%20%20%20%20np%2C%0A%20%20%20%20observed_signal%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20ridge_deconvolution%2C%0A%20%20%20%20signal_index%2C%0A%20%20%20%20signal_rmse%2C%0A%20%20%20%20spike_positions%2C%0A%20%20%20%20true_signal%2C%0A)%3A%0A%20%20%20%20gaussian_regularization%20%3D%2010.0%20**%20float(gaussian_log_lambda_control.value)%0A%20%20%20%20laplacian_regularization%20%3D%2010.0%20**%20float(laplacian_log_lambda_control.value)%0A%20%20%20%20naive_estimate%20%3D%20naive_deconvolution(observed_signal%2C%20blur_spectrum)%0A%20%20%20%20gaussian_estimate%20%3D%20ridge_deconvolution(%0A%20%20%20%20%20%20%20%20observed_signal%2C%0A%20%20%20%20%20%20%20%20blur_spectrum%2C%0A%20%20%20%20%20%20%20%20gaussian_regularization%2C%0A%20%20%20%20)%0A%20%20%20%20laplacian_estimate%2C%20fista_objective%20%3D%20fista_l1_deconvolution(%0A%20%20%20%20%20%20%20%20observed_signal%2C%0A%20%20%20%20%20%20%20%20blur_spectrum%2C%0A%20%20%20%20%20%20%20%20laplacian_regularization%2C%0A%20%20%20%20%20%20%20%20int(fista_iterations_control.value)%2C%0A%20%20%20%20)%0A%20%20%20%20naive_rmse%20%3D%20signal_rmse(true_signal%2C%20naive_estimate)%0A%20%20%20%20gaussian_rmse%20%3D%20signal_rmse(true_signal%2C%20gaussian_estimate)%0A%20%20%20%20laplacian_rmse%20%3D%20signal_rmse(true_signal%2C%20laplacian_estimate)%0A%20%20%20%20gaussian_residual_rms%20%3D%20signal_rmse(%0A%20%20%20%20%20%20%20%20observed_signal%2C%0A%20%20%20%20%20%20%20%20apply_blur(gaussian_estimate%2C%20blur_spectrum)%2C%0A%20%20%20%20)%0A%20%20%20%20laplacian_residual_rms%20%3D%20signal_rmse(%0A%20%20%20%20%20%20%20%20observed_signal%2C%0A%20%20%20%20%20%20%20%20apply_blur(laplacian_estimate%2C%20blur_spectrum)%2C%0A%20%20%20%20)%0A%20%20%20%20support_threshold%20%3D%200.15%0A%20%20%20%20gaussian_active_count%20%3D%20int(np.sum(np.abs(gaussian_estimate)%20%3E%20support_threshold))%0A%20%20%20%20laplacian_active_count%20%3D%20int(np.sum(np.abs(laplacian_estimate)%20%3E%20support_threshold))%0A%20%20%20%20_true_support%20%3D%20set(spike_positions.tolist())%0A%20%20%20%20_gaussian_support%20%3D%20set(np.flatnonzero(np.abs(gaussian_estimate)%20%3E%20support_threshold).tolist())%0A%20%20%20%20_laplacian_support%20%3D%20set(np.flatnonzero(np.abs(laplacian_estimate)%20%3E%20support_threshold).tolist())%0A%20%20%20%20gaussian_support_precision%20%3D%20len(_true_support%20%26%20_gaussian_support)%20%2F%20max(len(_gaussian_support)%2C%201)%0A%20%20%20%20laplacian_support_precision%20%3D%20len(_true_support%20%26%20_laplacian_support)%20%2F%20max(len(_laplacian_support)%2C%201)%0A%20%20%20%20gaussian_support_recall%20%3D%20len(_true_support%20%26%20_gaussian_support)%20%2F%20len(_true_support)%0A%20%20%20%20laplacian_support_recall%20%3D%20len(_true_support%20%26%20_laplacian_support)%20%2F%20len(_true_support)%0A%20%20%20%20_reconstruction_limit%20%3D%201.25%20*%20max(float(np.max(np.abs(true_signal)))%2C%200.5)%0A%20%20%20%20fig_reconstruction%2C%20_reconstruction_axes%20%3D%20plt.subplots(2%2C%202%2C%20figsize%3D(11.5%2C%206.5)%2C%20sharex%3DTrue)%0A%20%20%20%20_reconstruction_axes%5B0%2C%200%5D.vlines(signal_index%2C%200%2C%20true_signal%2C%20color%3D%22%23777777%22%2C%20lw%3D1.2%2C%20alpha%3D0.75%2C%20label%3D%22truth%22)%0A%20%20%20%20_reconstruction_axes%5B0%2C%200%5D.plot(signal_index%2C%20observed_signal%2C%20color%3D%22%234c78a8%22%2C%20lw%3D1.4%2C%20label%3D%22data%22)%0A%20%20%20%20_reconstruction_axes%5B0%2C%200%5D.set(title%3D%22Blurred%2C%20noisy%20data%22%2C%20ylabel%3D%22amplitude%22)%0A%20%20%20%20_reconstruction_axes%5B0%2C%200%5D.legend(frameon%3DTrue%2C%20fontsize%3D8)%0A%20%20%20%20_reconstruction_axes%5B0%2C%201%5D.vlines(signal_index%2C%200%2C%20true_signal%2C%20color%3D%22%23777777%22%2C%20lw%3D1.2%2C%20alpha%3D0.75)%0A%20%20%20%20_reconstruction_axes%5B0%2C%201%5D.plot(signal_index%2C%20naive_estimate%2C%20color%3D%22%23d1495b%22%2C%20lw%3D1.0)%0A%20%20%20%20_reconstruction_axes%5B0%2C%201%5D.set(%0A%20%20%20%20%20%20%20%20title%3Df%22Direct%20inverse%3A%20RMSE%20%7Bnaive_rmse%3A.1f%7D%20(off%20scale)%22%2C%0A%20%20%20%20%20%20%20%20ylabel%3D%22amplitude%22%2C%0A%20%20%20%20)%0A%20%20%20%20_reconstruction_axes%5B1%2C%200%5D.vlines(signal_index%2C%200%2C%20true_signal%2C%20color%3D%22%23777777%22%2C%20lw%3D1.2%2C%20alpha%3D0.75%2C%20label%3D%22truth%22)%0A%20%20%20%20_reconstruction_axes%5B1%2C%200%5D.plot(signal_index%2C%20gaussian_estimate%2C%20color%3D%22%23f28e2b%22%2C%20lw%3D1.8%2C%20label%3D%22Gaussian-prior%20MAP%22)%0A%20%20%20%20_reconstruction_axes%5B1%2C%200%5D.set(%0A%20%20%20%20%20%20%20%20title%3Df%22Quadratic%20%2F%20Gaussian%3A%20RMSE%20%7Bgaussian_rmse%3A.3f%7D%22%2C%0A%20%20%20%20%20%20%20%20xlabel%3D%22sample%20index%22%2C%0A%20%20%20%20%20%20%20%20ylabel%3D%22amplitude%22%2C%0A%20%20%20%20)%0A%20%20%20%20_reconstruction_axes%5B1%2C%200%5D.legend(frameon%3DTrue%2C%20fontsize%3D8)%0A%20%20%20%20_reconstruction_axes%5B1%2C%201%5D.vlines(signal_index%2C%200%2C%20true_signal%2C%20color%3D%22%23777777%22%2C%20lw%3D1.2%2C%20alpha%3D0.75%2C%20label%3D%22truth%22)%0A%20%20%20%20_reconstruction_axes%5B1%2C%201%5D.plot(signal_index%2C%20laplacian_estimate%2C%20color%3D%22%23285f8f%22%2C%20lw%3D1.8%2C%20label%3D%22Laplacian-prior%20MAP%22)%0A%20%20%20%20_reconstruction_axes%5B1%2C%201%5D.set(%0A%20%20%20%20%20%20%20%20title%3Df%22Absolute%20%2F%20Laplacian%3A%20RMSE%20%7Blaplacian_rmse%3A.3f%7D%22%2C%0A%20%20%20%20%20%20%20%20xlabel%3D%22sample%20index%22%2C%0A%20%20%20%20%20%20%20%20ylabel%3D%22amplitude%22%2C%0A%20%20%20%20)%0A%20%20%20%20_reconstruction_axes%5B1%2C%201%5D.legend(frameon%3DTrue%2C%20fontsize%3D8)%0A%20%20%20%20for%20_reconstruction_axis%20in%20_reconstruction_axes.flat%3A%0A%20%20%20%20%20%20%20%20_reconstruction_axis.set_ylim(-_reconstruction_limit%2C%20_reconstruction_limit)%0A%20%20%20%20%20%20%20%20_reconstruction_axis.grid(alpha%3D0.24)%0A%20%20%20%20fig_reconstruction.suptitle(%0A%20%20%20%20%20%20%20%20%22Regularization%20suppresses%20unstable%20directions%3B%20the%20prior%20chooses%20the%20shape%22%2C%0A%20%20%20%20%20%20%20%20fontweight%3D%22bold%22%2C%0A%20%20%20%20%20%20%20%20y%3D1.01%2C%0A%20%20%20%20)%0A%20%20%20%20fig_reconstruction.tight_layout()%0A%20%20%20%20fig_reconstruction%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20fista_objective%2C%0A%20%20%20%20%20%20%20%20gaussian_active_count%2C%0A%20%20%20%20%20%20%20%20gaussian_estimate%2C%0A%20%20%20%20%20%20%20%20gaussian_residual_rms%2C%0A%20%20%20%20%20%20%20%20gaussian_rmse%2C%0A%20%20%20%20%20%20%20%20gaussian_support_precision%2C%0A%20%20%20%20%20%20%20%20gaussian_support_recall%2C%0A%20%20%20%20%20%20%20%20laplacian_active_count%2C%0A%20%20%20%20%20%20%20%20laplacian_estimate%2C%0A%20%20%20%20%20%20%20%20laplacian_residual_rms%2C%0A%20%20%20%20%20%20%20%20laplacian_rmse%2C%0A%20%20%20%20%20%20%20%20laplacian_support_precision%2C%0A%20%20%20%20%20%20%20%20laplacian_support_recall%2C%0A%20%20%20%20%20%20%20%20support_threshold%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20reconstruction_metrics(%0A%20%20%20%20gaussian_active_count%2C%0A%20%20%20%20gaussian_residual_rms%2C%0A%20%20%20%20gaussian_rmse%2C%0A%20%20%20%20gaussian_support_precision%2C%0A%20%20%20%20gaussian_support_recall%2C%0A%20%20%20%20laplacian_active_count%2C%0A%20%20%20%20laplacian_residual_rms%2C%0A%20%20%20%20laplacian_rmse%2C%0A%20%20%20%20laplacian_support_precision%2C%0A%20%20%20%20laplacian_support_recall%2C%0A%20%20%20%20mo%2C%0A%20%20%20%20support_threshold%2C%0A)%3A%0A%20%20%20%20mo.md(rf%22%22%22%0A%20%20%20%20%7C%20Estimate%20%7C%20Signal%20RMSE%20%7C%20Data%20residual%20RMS%20%7C%20Active%20coefficients%20%7C%20Support%20precision%20%7C%20Support%20recall%20%7C%0A%20%20%20%20%7C---%7C---%3A%7C---%3A%7C---%3A%7C---%3A%7C---%3A%7C%0A%20%20%20%20%7C%20Gaussian%20%2F%20%24%5Cell_2%24%20%7C%20%7Bgaussian_rmse%3A.4f%7D%20%7C%20%7Bgaussian_residual_rms%3A.4f%7D%20%7C%20%7Bgaussian_active_count%7D%20%7C%20%7Bgaussian_support_precision%3A.2f%7D%20%7C%20%7Bgaussian_support_recall%3A.2f%7D%20%7C%0A%20%20%20%20%7C%20Laplacian%20%2F%20%24%5Cell_1%24%20%7C%20%7Blaplacian_rmse%3A.4f%7D%20%7C%20%7Blaplacian_residual_rms%3A.4f%7D%20%7C%20%7Blaplacian_active_count%7D%20%7C%20%7Blaplacian_support_precision%3A.2f%7D%20%7C%20%7Blaplacian_support_recall%3A.2f%7D%20%7C%0A%0A%20%20%20%20%E2%80%9CActive%E2%80%9D%20means%20%24%7C%5Cwidehat%20x%5Bk%5D%7C%3E%7Bsupport_threshold%3A.2f%7D%24.%20A%20lower%20residual%20only%20means%20a%20closer%20fit%0A%20%20%20%20to%20the%20noisy%20data%3B%20the%20hidden-signal%20RMSE%20and%20support%20columns%20reveal%20whether%20that%20fit%20recovered%0A%20%20%20%20the%20sparse%20cause.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20solver_diagnostics(%0A%20%20%20%20fista_objective%2C%0A%20%20%20%20gaussian_estimate%2C%0A%20%20%20%20laplacian_estimate%2C%0A%20%20%20%20np%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20signal_length%2C%0A%20%20%20%20support_threshold%2C%0A)%3A%0A%20%20%20%20_fista_iterations%20%3D%20np.arange(1%2C%20len(fista_objective)%20%2B%201)%0A%20%20%20%20_gaussian_sorted%20%3D%20np.sort(np.abs(gaussian_estimate))%5B%3A%3A-1%5D%0A%20%20%20%20_laplacian_sorted%20%3D%20np.sort(np.abs(laplacian_estimate))%5B%3A%3A-1%5D%0A%20%20%20%20_coefficient_rank%20%3D%20np.arange(1%2C%20signal_length%20%2B%201)%0A%20%20%20%20fig_diagnostics%2C%20_diagnostic_axes%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(10.8%2C%203.8))%0A%20%20%20%20_diagnostic_axes%5B0%5D.plot(_fista_iterations%2C%20fista_objective%2C%20color%3D%22%23285f8f%22%2C%20lw%3D1.9)%0A%20%20%20%20_diagnostic_axes%5B0%5D.set(%0A%20%20%20%20%20%20%20%20title%3Df%22FISTA%20objective%20after%20%7Blen(fista_objective)%7D%20iterations%22%2C%0A%20%20%20%20%20%20%20%20xlabel%3D%22iteration%22%2C%0A%20%20%20%20%20%20%20%20ylabel%3D%22data%20term%20%2B%20Laplacian%20penalty%22%2C%0A%20%20%20%20)%0A%20%20%20%20_diagnostic_axes%5B1%5D.semilogy(_coefficient_rank%2C%20_gaussian_sorted%20%2B%201e-10%2C%20color%3D%22%23f28e2b%22%2C%20lw%3D2.0%2C%20label%3D%22Gaussian%22)%0A%20%20%20%20_diagnostic_axes%5B1%5D.semilogy(_coefficient_rank%2C%20_laplacian_sorted%20%2B%201e-10%2C%20color%3D%22%23285f8f%22%2C%20lw%3D2.0%2C%20label%3D%22Laplacian%22)%0A%20%20%20%20_diagnostic_axes%5B1%5D.axhline(support_threshold%2C%20color%3D%22%23d1495b%22%2C%20ls%3D%22--%22%2C%20lw%3D1.3%2C%20label%3D%22active%20threshold%22)%0A%20%20%20%20_diagnostic_axes%5B1%5D.set(%0A%20%20%20%20%20%20%20%20title%3D%22Ranked%20reconstruction%20coefficients%22%2C%0A%20%20%20%20%20%20%20%20xlabel%3D%22coefficient%20rank%22%2C%0A%20%20%20%20%20%20%20%20ylabel%3D%22absolute%20amplitude%22%2C%0A%20%20%20%20%20%20%20%20ylim%3D(1e-9%2C%20max(_gaussian_sorted%5B0%5D%2C%20_laplacian_sorted%5B0%5D)%20*%201.4)%2C%0A%20%20%20%20)%0A%20%20%20%20_diagnostic_axes%5B1%5D.legend(frameon%3DTrue%2C%20fontsize%3D9)%0A%20%20%20%20for%20_diagnostic_axis%20in%20_diagnostic_axes%3A%0A%20%20%20%20%20%20%20%20_diagnostic_axis.grid(alpha%3D0.24)%0A%20%20%20%20fig_diagnostics.suptitle(%22Numerical%20convergence%20and%20recovered%20sparsity%22%2C%20fontweight%3D%22bold%22%2C%20y%3D1.03)%0A%20%20%20%20fig_diagnostics.tight_layout()%0A%20%20%20%20fig_diagnostics%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20reconstruction_summary(%0A%20%20%20%20gaussian_active_count%2C%0A%20%20%20%20gaussian_rmse%2C%0A%20%20%20%20laplacian_active_count%2C%0A%20%20%20%20laplacian_rmse%2C%0A%20%20%20%20mo%2C%0A)%3A%0A%20%20%20%20_better_prior%20%3D%20%22Laplacian%22%20if%20laplacian_rmse%20%3C%20gaussian_rmse%20else%20%22Gaussian%22%0A%20%20%20%20mo.callout(%0A%20%20%20%20%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22At%20the%20current%20settings%2C%20the%20**%7B_better_prior%7D%20prior**%20has%20lower%20reconstruction%20RMSE.%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22The%20Gaussian%20solution%20retains%20%7Bgaussian_active_count%7D%20active%20coefficients%3B%20the%20Laplacian%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20f%22solution%20retains%20%7Blaplacian_active_count%7D.%20The%20point%20is%20not%20that%20one%20prior%20always%20wins%3A%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%22the%20prior%20should%20match%20the%20structure%20believed%20to%20have%20generated%20the%20unknown%20signal.%22%0A%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20kind%3D%22success%22%2C%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20recap(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20What%20to%20take%20away%0A%0A%20%20%20%20-%20Blur%20makes%20inversion%20ill-conditioned%20by%20suppressing%20frequency%20directions%3B%20direct%20inversion%0A%20%20%20%20%20%20amplifies%20noise%20in%20precisely%20those%20directions.%0A%20%20%20%20-%20Quadratic%20Tikhonov%20regularization%20is%20the%20MAP%20estimator%20for%20a%20Gaussian%20prior.%20It%20has%20a%20linear%2C%0A%20%20%20%20%20%20closed-form%20solution%20and%20produces%20smooth%20shrinkage.%0A%20%20%20%20-%20An%20%24%5Cell_1%24%20penalty%20is%20the%20MAP%20estimator%20for%20a%20Laplacian%20prior.%20Its%20cusp%20produces%20soft%0A%20%20%20%20%20%20thresholding%20and%20exact%20zeros%2C%20but%20the%20estimate%20must%20be%20found%20numerically.%0A%20%20%20%20-%20A%20sparse%20signal%20rewards%20the%20sparse%20prior.%20Change%20to%20weaker%20blur%20or%20stronger%20regularization%20and%0A%20%20%20%20%20%20watch%20when%20that%20advantage%20disappears.%0A%0A%20%20%20%20A%20prior%20is%20therefore%20not%20merely%20a%20numerical%20stabilizer.%20It%20is%20a%20claim%20about%20what%20kinds%20of%20hidden%0A%20%20%20%20signals%20are%20plausible.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A