import%20marimo%0A%0A__generated_with%20%3D%20%220.24.0%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Blind%20Deconvolution%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20warnings%0A%20%20%20%20from%20pathlib%20import%20Path%0A%0A%20%20%20%20warnings.filterwarnings(%22ignore%22)%0A%0A%20%20%20%20import%20cvxpy%20as%20cp%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%20from%20dbcp%20import%20BiconvexProblem%2C%20convolve%0A%0A%20%20%20%20_example_directory%20%3D%20Path(__file__).resolve().parent%0A%20%20%20%20plt.style.use(_example_directory%20%2F%20%22zhlatex.mplstyle%22)%0A%20%20%20%20figure_directory%20%3D%20_example_directory%20%2F%20%22figures%22%0A%20%20%20%20figure_directory.mkdir(parents%3DTrue%2C%20exist_ok%3DTrue)%0A%0A%20%20%20%20np.random.seed(10015)%0A%20%20%20%20return%20BiconvexProblem%2C%20convolve%2C%20cp%2C%20figure_directory%2C%20mo%2C%20np%2C%20plt%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Introduction%0A%0A%20%20%20%20Blind%20deconvolution%20is%20a%20technique%20used%20to%20recover%20some%20sharp%20signal%20or%20image%20from%20a%20blurred%20observation%0A%20%20%20%20when%20the%20blur%20itself%20is%20unknown.%0A%20%20%20%20It%20jointly%20estimates%20both%20the%20original%20signal%20and%20the%20blur%20kernel%2C%20with%20some%20prior%20knowledge%20about%20their%20structures.%0A%0A%20%20%20%20Suppose%20we%20are%20given%20a%20data%20vector%20%24d%20%5Cin%20%5Cmathbf%7BR%7D%5E%7Bm%20%2B%20n%20-%201%7D%24%2C%20which%20is%20the%20convolution%20of%20an%20unknown%0A%20%20%20%20sparse%20signal%20%24x%20%5Cin%20%5Cmathbf%7BR%7D%5En%24%20and%20an%20unknown%20smooth%20vector%20%24y%20%5Cin%20%5Cmathbf%7BR%7D%5Em%24%20with%20bounded%0A%20%20%20%20%24%5Cell_%5Cinfty%24-norm%20(i.e.%2C%20bounded%20largest%20entry).%0A%20%20%20%20Additionally%2C%20we%20have%20the%20prior%20knowledge%20that%20both%20the%20vectors%20%24x%24%20and%20%24y%24%20are%20nonnegative.%0A%20%20%20%20The%20corresponding%20blind%20deconvolution%20problem%20can%20be%20formulated%20as%20the%20following%20biconvex%20optimization%20problem%3A%0A%0A%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%5Cbegin%7Barray%7D%7Bll%7D%0A%20%20%20%20%20%20%20%20%20%20%20%20%5Ctext%7Bminimize%7D%20%26%20%7B%5C%7Cx%20%5Cotimes%20%20y%20-%20d%5C%7C%7D_2%5E2%20%2B%20%5Calpha_%7B%5Crm%20sp%7D%20%7B%5C%7Cx%5C%7C%7D_1%20%2B%20%5Calpha_%7B%5Crm%20sm%7D%20%7B%5C%7CDy%5C%7C%7D_2%5E2%5C%5C%0A%20%20%20%20%20%20%20%20%20%20%20%20%5Ctext%7Bsubject%20to%7D%20%26%20x%20%5Csucceq%200%2C%5Cquad%20y%20%5Csucceq%200%5C%5C%0A%20%20%20%20%20%20%20%20%20%20%20%20%26%20%7B%5C%7Cy%5C%7C%7D_%5Cinfty%20%5Cleq%20%5Cbeta%0A%20%20%20%20%20%20%20%20%5Cend%7Barray%7D%0A%20%20%20%20%5C%5D%0A%0A%20%20%20%20with%20variables%20%24x%24%20and%20%24y%24%2C%20where%20%24%5Calpha_%7B%5Crm%20sp%7D%2C%20%5Calpha_%7B%5Crm%20sm%7D%20%3E%200%24%20are%20the%20regularization%20parameters%0A%20%20%20%20for%20the%20sparsity%20of%20%24x%24%20and%20smoothness%20of%20%24y%24%2C%20respectively%2C%20and%20%24%5Cbeta%20%3E%200%24%20is%20the%20bound%20on%20the%0A%20%20%20%20%24%5Cell_%5Cinfty%24-norm%20of%20the%20vector%20%24y%24.%0A%20%20%20%20The%20matrix%20%24D%20%5Cin%20%5Cmathbf%7BR%7D%5E%7B(m%20-%201)%20%5Ctimes%20m%7D%24%20is%20the%20first-order%20difference%20operator%2C%20given%20by%2C%0A%0A%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20D%20%3D%20%5Cleft%5B%5Cbegin%7Barray%7D%7Bccccc%7D%0A%20%20%20%20%20%20%20%20%20%20%20%201%20%26%20-1%20%26%26%26%5C%5C%0A%20%20%20%20%20%20%20%20%20%20%20%20%26%201%20%26%20-1%20%26%26%5C%5C%0A%20%20%20%20%20%20%20%20%20%20%20%20%26%26%20%5Cddots%20%26%20%5Cddots%20%26%5C%5C%0A%20%20%20%20%20%20%20%20%20%20%20%20%26%26%26%201%20%26%20-1%0A%20%20%20%20%20%20%20%20%5Cend%7Barray%7D%5Cright%5D%20%5Cin%20%5Cmathbf%7BR%7D%5E%7B(m%20-%201)%20%5Ctimes%20m%7D%2C%0A%20%20%20%20%5C%5D%0A%0A%20%20%20%20so%20that%20%24Dy%24%20computes%20the%20vector%20of%20successive%20differences%20of%20%24y%24.%0A%20%20%20%20The%20convolution%20%24x%20%5Cotimes%20y%24%20of%20the%20vectors%20%24x%24%20and%20%24y%24%20is%20given%20by%0A%0A%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%7B(x%20%5Cotimes%20y)%7D_k%20%3D%20%5Csum_%7Bi%20%2B%20j%20%3D%20k%20%2B%201%7D%20x_i%20y_j%2C%5Cquad%20k%20%3D%201%2C%20%5Cldots%2C%20m%20%2B%20n%20-%201.%0A%20%20%20%20%5C%5D%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Generate%20problem%20data%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20n%20%3D%20120%0A%20%20%20%20m%20%3D%2040%0A%0A%20%20%20%20x0%20%3D%20np.zeros(n)%0A%20%20%20%20x0%5B6%5D%20%3D%201%0A%20%20%20%20y0%20%3D%20np.exp(-np.square(np.linspace(-2%2C%202%2C%20m))%20*%202)%0A%20%20%20%20d%20%3D%20np.convolve(x0%2C%20y0)%0A%20%20%20%20return%20d%2C%20m%2C%20n%2C%20x0%2C%20y0%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Specify%20and%20solve%20the%20problem%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(BiconvexProblem%2C%20convolve%2C%20cp%2C%20d%2C%20m%2C%20n)%3A%0A%20%20%20%20alpha_sp%20%3D%200.1%0A%20%20%20%20alpha_sm%20%3D%200.2%0A%20%20%20%20beta%20%3D%201%0A%0A%20%20%20%20x%20%3D%20cp.Variable(n%2C%20nonneg%3DTrue)%0A%20%20%20%20y%20%3D%20cp.Variable(m%2C%20nonneg%3DTrue)%0A%20%20%20%20obj%20%3D%20cp.Minimize(%0A%20%20%20%20%20%20%20%20cp.sum_squares(convolve(x%2C%20y)%20-%20d)%20%2B%20alpha_sp%20*%20cp.norm1(x)%20%2B%20alpha_sm%20*%20cp.sum_squares(cp.diff(y))%0A%20%20%20%20)%0A%20%20%20%20constr%20%3D%20%5Bcp.norm(y%2C%20%22inf%22)%20%3C%3D%20beta%5D%0A%20%20%20%20prob%20%3D%20BiconvexProblem(obj%2C%20%5Bx%5D%2C%20%5By%5D%2C%20constr)%0A%20%20%20%20prob.solve(cp.CLARABEL%2C%20abs_tol%3D1e-5%2C%20max_iter%3D200)%0A%20%20%20%20return%20x%2C%20y%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Plot%20the%20results%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(d%2C%20figure_directory%2C%20np%2C%20plt%2C%20x%2C%20x0%2C%20y%2C%20y0)%3A%0A%20%20%20%20fig%2C%20axs%20%3D%20plt.subplots(1%2C%201%2C%20figsize%3D(6%2C%204.5))%0A%20%20%20%20axs.plot(x0%2C%20linestyle%3D%22--%22%2C%20color%3D%22C3%22%2C%20linewidth%3D2)%0A%20%20%20%20axs.plot(y0%2C%20linestyle%3D%22--%22%2C%20color%3D%22C1%22%2C%20linewidth%3D2)%0A%20%20%20%20axs.plot(d%2C%20linestyle%3D%22--%22%2C%20color%3D%22k%22%2C%20linewidth%3D2)%0A%20%20%20%20axs.plot(x.value%2C%20color%3D%22C0%22%2C%20marker%3D%22.%22%2C%20markersize%3D10)%0A%20%20%20%20axs.plot(y.value%2C%20color%3D%22C2%22%2C%20marker%3D%22s%22)%0A%20%20%20%20axs.plot(np.convolve(x.value%2C%20y.value)%2C%20marker%3D%22D%22%2C%20color%3D%22C4%22%2C%20zorder%3D-1)%0A%0A%20%20%20%20axs.legend(%0A%20%20%20%20%20%20%20%20%5B%22ground%20truth%20%24x%24%22%2C%20%22ground%20truth%20%24y%24%22%2C%20%22ground%20truth%20%24d%24%22%2C%20%22recovered%20%24x%24%22%2C%20%22recovered%20%24y%24%22%2C%20%22recovered%20%24d%24%22%5D%2C%0A%20%20%20%20%20%20%20%20frameon%3DFalse%2C%0A%20%20%20%20%20%20%20%20fontsize%3D12%2C%0A%20%20%20%20)%0A%20%20%20%20axs.set_xlim(0%2C%2060)%0A%20%20%20%20axs.set_xlabel(%22%24i%24%22)%0A%0A%20%20%20%20fig.tight_layout()%0A%20%20%20%20fig.savefig(figure_directory%20%2F%20%22blind_deconv.pdf%22%2C%20bbox_inches%3D%22tight%22)%0A%20%20%20%20plt.show()%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
e320b4031176308e21c2c134b1414050