Complex Least Squares

$ % macros % MathJax \newcommand{\bigtimes}{\mathop{\vcenter{\hbox{$\Huge\times\normalsize$}}}} % prefix with \! and postfix with \!\! % sized grouping symbols \renewcommand {\aa} [1] {\left\langle {#1} \right\rangle} % <> angle brackets \newcommand {\bb} [1] {\left[ {#1} \right]} % [] brackets \newcommand {\cc} [1] {\left\{ {#1} \right\}} % {} curly braces \newcommand {\mm} [1] {\lVert {#1} \rVert} % || || double norm bars \newcommand {\nn} [1] {\lvert {#1} \rvert} % || norm bars \newcommand {\pp} [1] {\left( {#1} \right)} % () parentheses % unit \newcommand {\unit} [1] {\bb{\mathrm{#1}}} % measurement unit % math \newcommand {\fn} [1] {\mathrm{#1}} % function name % sets \newcommand {\setZ} {\mathbb{Z}} \newcommand {\setQ} {\mathbb{Q}} \newcommand {\setR} {\mathbb{R}} \newcommand {\setC} {\mathbb{C}} % arithmetic \newcommand {\q} [2] {\frac{#1}{#2}} % quotient, because fuck \frac % trig \newcommand {\acos} {\mathrm{acos}} % \mathrm{???}^{-1} is misleading \newcommand {\asin} {\mathrm{asin}} \newcommand {\atan} {\mathrm{atan}} \newcommand {\atantwo} {\mathrm{atan2}} % at angle = atan2(y, x) \newcommand {\asec} {\mathrm{asec}} \newcommand {\acsc} {\mathrm{acsc}} \newcommand {\acot} {\mathrm{acot}} % complex numbers \newcommand {\z} [1] {\tilde{#1}} \newcommand {\conj} [1] {{#1}^\ast} \renewcommand {\Re} {\mathfrak{Re}} % real part \renewcommand {\Im} {\mathrm{I}\mathfrak{m}} % imaginary part % quaternions \newcommand {\quat} [1] {\tilde{\mathbf{#1}}} % quaternion symbol \newcommand {\uquat} [1] {\check{\mathbf{#1}}} % versor symbol \newcommand {\gquat} [1] {\tilde{\boldsymbol{#1}}} % greek quaternion symbol \newcommand {\guquat}[1] {\check{\boldsymbol{#1}}} % greek versor symbol % vectors \renewcommand {\vec} [1] {\mathbf{#1}} % vector symbol \newcommand {\uvec} [1] {\hat{\mathbf{#1}}} % unit vector symbol \newcommand {\gvec} [1] {\boldsymbol{#1}} % greek vector symbol \newcommand {\guvec} [1] {\hat{\boldsymbol{#1}}} % greek unit vector symbol % special math vectors \renewcommand {\r} {\vec{r}} % r vector [m] \newcommand {\R} {\vec{R}} % R = r - r' difference vector [m] \newcommand {\ur} {\uvec{r}} % r unit vector [#] \newcommand {\uR} {\uvec{R}} % R unit vector [#] \newcommand {\ux} {\uvec{x}} % x unit vector [#] \newcommand {\uy} {\uvec{y}} % y unit vector [#] \newcommand {\uz} {\uvec{z}} % z unit vector [#] \newcommand {\urho} {\guvec{\rho}} % rho unit vector [#] \newcommand {\utheta} {\guvec{\theta}} % theta unit vector [#] \newcommand {\uphi} {\guvec{\phi}} % phi unit vector [#] \newcommand {\un} {\uvec{n}} % unit normal vector [#] % vector operations \newcommand {\inner} [2] {\left\langle {#1} , {#2} \right\rangle} % <,> \newcommand {\outer} [2] {{#1} \otimes {#2}} \newcommand {\norm} [1] {\mm{#1}} \renewcommand {\dot} {\cdot} % dot product \newcommand {\cross} {\times} % cross product % matrices \newcommand {\mat} [1] {\mathbf{#1}} % matrix symbol \newcommand {\gmat} [1] {\boldsymbol{#1}} % greek matrix symbol % ordinary derivatives \newcommand {\od} [2] {\q{d #1}{d #2}} % ordinary derivative \newcommand {\odn} [3] {\q{d^{#3}{#1}}{d{#2}^{#3}}} % nth order od \newcommand {\odt} [1] {\q{d{#1}}{dt}} % time od % partial derivatives \newcommand {\de} {\partial} % partial symbol \newcommand {\pd} [2] {\q{\de{#1}}{\de{#2}}} % partial derivative \newcommand {\pdn} [3] {\q{\de^{#3}{#1}}{\de{#2}^{#3}}} % nth order pd \newcommand {\pdd} [3] {\q{\de^2{#1}}{\de{#2}\de{#3}}} % 2nd order mixed pd \newcommand {\pdt} [1] {\q{\de{#1}}{\de{t}}} % time pd % vector derivatives \newcommand {\del} {\nabla} % del \newcommand {\grad} {\del} % gradient \renewcommand {\div} {\del\dot} % divergence \newcommand {\curl} {\del\cross} % curl % differential vectors \newcommand {\dL} {d\vec{L}} % differential vector length [m] \newcommand {\dS} {d\vec{S}} % differential vector surface area [m^2] % special functions \newcommand {\Hn} [2] {H^{(#1)}_{#2}} % nth order Hankel function \newcommand {\hn} [2] {H^{(#1)}_{#2}} % nth order spherical Hankel function % transforms \newcommand {\FT} {\mathcal{F}} % fourier transform \newcommand {\IFT} {\FT^{-1}} % inverse fourier transform % signal processing \newcommand {\conv} [2] {{#1}\ast{#2}} % convolution \newcommand {\corr} [2] {{#1}\star{#2}} % correlation % abstract algebra \newcommand {\lie} [1] {\mathfrak{#1}} % lie algebra % other \renewcommand {\d} {\delta} % optimization %\DeclareMathOperator* {\argmin} {arg\,min} %\DeclareMathOperator* {\argmax} {arg\,max} \newcommand {\argmin} {\fn{arg\,min}} \newcommand {\argmax} {\fn{arg\,max}} % waves \renewcommand {\l} {\lambda} % wavelength [m] \renewcommand {\k} {\vec{k}} % wavevector [rad/m] \newcommand {\uk} {\uvec{k}} % unit wavevector [#] \newcommand {\w} {\omega} % angular frequency [rad/s] \renewcommand {\TH} {e^{j \w t}} % engineering time-harmonic function [#] % classical mechanics \newcommand {\F} {\vec{F}} % force [N] \newcommand {\p} {\vec{p}} % momentum [kg m/s] % \r % position [m], aliased \renewcommand {\v} {\vec{v}} % velocity vector [m/s] \renewcommand {\a} {\vec{a}} % acceleration [m/s^2] \newcommand {\vGamma} {\gvec{\Gamma}} % torque [N m] \renewcommand {\L} {\vec{L}} % angular momentum [kg m^2 / s] \newcommand {\mI} {\mat{I}} % moment of inertia tensor [kg m^2/rad] \newcommand {\vw} {\gvec{\omega}} % angular velocity vector [rad/s] \newcommand {\valpha} {\gvec{\alpha}} % angular acceleration vector [rad/s^2] % electromagnetics % fields \newcommand {\E} {\vec{E}} % electric field intensity [V/m] \renewcommand {\H} {\vec{H}} % magnetic field intensity [A/m] \newcommand {\D} {\vec{D}} % electric flux density [C/m^2] \newcommand {\B} {\vec{B}} % magnetic flux density [Wb/m^2] % potentials \newcommand {\A} {\vec{A}} % vector potential [Wb/m], [C/m] % \F % electric vector potential [C/m], aliased % sources \newcommand {\I} {\vec{I}} % line current density [A] , [V] \newcommand {\J} {\vec{J}} % surface current density [A/m] , [V/m] \newcommand {\K} {\vec{K}} % volume current density [A/m^2], [V/m^2] % \M % magnetic current [V/m^2], aliased, obsolete % materials \newcommand {\ep} {\epsilon} % permittivity [F/m] % \mu % permeability [H/m], aliased \renewcommand {\P} {\vec{P}} % polarization [C/m^2], [Wb/m^2] % \p % electric dipole moment [C m], aliased \newcommand {\M} {\vec{M}} % magnetization [A/m], [V/m] \newcommand {\m} {\vec{m}} % magnetic dipole moment [A m^2] % power \renewcommand {\S} {\vec{S}} % poynting vector [W/m^2] \newcommand {\Sa} {\aa{\vec{S}}_t} % time averaged poynting vector [W/m^2] % quantum mechanics \newcommand {\bra} [1] {\left\langle {#1} \right|} % <| \newcommand {\ket} [1] {\left| {#1} \right\rangle} % |> \newcommand {\braket} [2] {\left\langle {#1} \middle| {#2} \right\rangle} $

Table of Contents

Problem Statement

Given matrix $\mat{A} \in \setC^{m \times n}$ and vector $\vec{b} \in \setC^m$, find vector $\vec{x} \in \setC^n$ such that $$ \label{eqn:lsq_Axb} \mat{A} \vec{x} = \vec{b} . $$ Usually this is an ill-posed problem having no solutions or infinite solutions.

While it is not always possible to satisfy Eq. $\eqref{eqn:lsq_Axb}$, it is always possible to take the difference of both sides. Define the residual as $$ \label{eqn:residual} \boxed{\vec{r}(\vec{x}) \equiv \vec{b} - \mat{A} \vec{x}} . $$ A reasonable, although somewhat arbitrary, way to compare the fitness of any $\vec{x}$ as an approximate solution to Eq. $\eqref{eqn:lsq_Axb}$ is by the residual norm (squared) $$ \label{eqn:cost} f(\vec{x}) = \norm{\vec{b} - \mat{A} \vec{x}}^2 $$ where the norm is taken to be that induced by the standard complex inner product $$ \aa{\vec{x}, \vec{y}} = \vec{x}^H \vec{y} \quad\Rightarrow\quad \norm{\vec{x}}^2 = \aa{\vec{x}, \vec{x}} = \vec{x}^H \vec{x} . $$ A smaller residual norm is interpreted as indicating a better approximation to Eq. $\eqref{eqn:lsq_Axb}$, and zero residual norm an exact solution. Function $f : \setC^n \rightarrow \setR$ is a cost function for which a global minimum is sought. This is the complex least squares problem. A solution vector $\vec{x}_s \in \setC^n$ is defined to minimize the cost function $$ \vec{x}_s = \argmin_{\vec{x} \in \setC^n} \mm{\vec{b} - \mat{A} \vec{x}}^2 . $$

Gradient-based methods are often used to solve optimization problems []. For a real-valued function of real variables $f : \setR^n \rightarrow \setR$ the gradient points in the direction of greatest change, and extrema are necessarily at critical points with zero gradient $$ \grad f (\vec{x}_s) = \vec{0} . $$ It is desirable to adapt these ideas to complex functions. Recall for $\setR^3$ the gradient is motivated by writing the total derivative of $f$ in terms of a dot product (a kind of inner product) $$ \label{eqn:grad} df = \pd{f}{x} dx + \pd{f}{y} dy + \pd{f}{z} dz = \bb{\ux \pd{f}{x} + \uy \pd{f}{y} + \uz \pd{f}{z}} \dot \bb{\ux dx + \uy dy + \uz dz} = \grad{f} \dot d \r . $$ It would seem a complex gradient could be defined in terms of complex derivatives and the standard complex inner product, and this is true. However the complex least squares cost function is not complex-differentiable despite being continuous. Additional tools of calculus are needed to generalize the gradient to this case.

Solution

Recall the complex least squares cost function to be minimized is the residual norm $$ f(\vec{x}) = \norm{\vec{b} - \mat{A} \vec{x}}^2 . $$ Take the CR gradient of the cost function $$ \grad f(\vec{x}) = \grad \mm{\vec{b} - \mat{A} \vec{x}}^2 $$ convert to summation notation to evaluate

$$ \label{eqn:lsq_grad} \begin{split} \conj{\bb{\pd{f}{x_k}}} &= \conj{\bb{\pd{}{x_k} \sum_{j=1}^m \nn{b_j - \sum_{i=1}^n A_{ji} x_i}^2}} \\ &= \conj{\bb{\sum_{j=1}^m \pd{}{x_k} \nn{b_j - \sum_{i=1}^n A_{ji} x_i}^2}} \\ &= \conj{\bb{\sum_{j=1}^m \conj{\bb{b_j - \sum_{i=1}^n A_{ji} x_i}} \pd{}{x_k} \bb{b_j -\sum_{i=1}^n A_{ji} x_i}}} \\ &= \conj{\bb{\sum_{j=1}^m \conj{\bb{b_j - \sum_{i=1}^n A_{ji} x_i}} \bb{\pd{b_j}{x_k} - \sum_{i=1}^n A_{ji} \pd{x_i}{x_k}}}} \\ &= \conj{\bb{\sum_{j=1}^m \conj{\bb{b_j - \sum_{i=1}^n A_{ji} x_i}} \bb{0 - \sum_{i=1}^n A_{ji} \delta_{ik}}}} \\ &= \conj{\bb{\sum_{j=1}^m \conj{\bb{b_j - \sum_{i=1}^n A_{ji} x_i}} \bb{-A_{jk}}}} \\ &= -\sum_{j=1}^m \conj{A_{jk}} \bb{b_j - \sum_{i=1}^n A_{ji} x_i} \end{split} $$

and convert back to matrix notation to be concise $$ \label{eqn:grad_cost} \grad f(\vec{x}) = -\mat{A}^H \bb{\vec{b} - \mat{A} \vec{x}} . $$ Optimize the cost function $f$ by setting its CR gradient to zero $$ \grad f(\vec{x}_s) = - \mat{A}^H [\vec{b} - \mat{A} \vec{x}_s] = \vec{0} $$ to derive the complex least squares equation $$ \label{eqn:lsq} \boxed{\mat{A}^H \mat{A} \vec{x}_s = \mat{A}^H \vec{b}} . $$ This is equivalent to the modified linear system $$ \mat{A}' \vec{x}_s = \vec{b}' $$ where $\mat{A}' = \mat{A}^H \mat{A}$ and $\vec{b}' = \mat{A}^H \vec{b}$. These are also known as normal equations []. Geometrically a complex least squares solution vector $\vec{x}_s$ makes the solution residual vector $\vec{r}_s$ orthogonal to the columnspace (output space) of $\mat{A}$ $$ \mat{A}^H [\vec{b} - \mat{A} \vec{x}_s] = \mat{A}^H \vec{r}_s = \vec{0} $$ meaning the solution residual only has components in the space inaccessible to $\mat{A}$.

Of passing interest is the real least squares cost function $f : \setR^n \rightarrow \setR$ $$ f(\vec{x}) = \q{1}{2} \norm{\vec{b} - \mat{A} \vec{x}}^2 $$ where a factor of $1/2$ is introduced to compensate for differences between the real and CR gradients. The real gradient is computed in a similar way as Eq. $\eqref{eqn:lsq_grad}$ except instead of the CR chain rule the real chain rule is used to differentiate the square and cancel the $1/2$, resulting in $$ \grad f(\vec{x}) = - \mat{A}^T [\vec{b} - \mat{A} \vec{x}] . $$ Optimize function $f$ by setting the gradient to zero $$ \grad f(\vec{x}_s) = - \mat{A}^T [\vec{b} - \mat{A} \vec{x}_s] = \vec{0} $$ to form the real least squares equation $$ \label{eqn:real_lsq} \mat{A}^T \mat{A} \vec{x}_s = \mat{A}^T \vec{b} . $$ The only difference between Eq. $\eqref{eqn:lsq}$ and Eq. $\eqref{eqn:real_lsq}$ is that the conjugate transposes are replaced by transposes. The similarity between these results belies the important theoretical differences in their derivation.

Regularization

A problem with the complex least squares cost function Eq. $\eqref{eqn:cost}$ is that infinite solution vectors exist if $\mat{A}$ has a non-trivial nullspace. To see this, evaluate the cost function $f$ at a solution vector $\vec{x}_s$ displaced by any non-zero nullspace vector $\vec{e} \in \fn{nullspace}(\mat{A}) - \cc{\vec{0}} \subset \setC^n$ such that $$ f(\vec{x}_s + \vec{e}) = \norm{\vec{b} - \mat{A} [\vec{x}_s + \vec{e}]}^2 %= \norm{\vec{b} - \mat{A} \vec{x}_s - \mat{A} \vec{e}}^2 = \norm{\vec{b} - \mat{A} \vec{x}_s - \vec{0}}^2 = \norm{\vec{b} - \mat{A} \vec{x}_s}^2 = f(\vec{x}_s) . $$ Thus any non-zero nullspace vector can be added to a solution vector to yield a different solution vector with the same minimal cost function value.

Some additional condition must be imposed to guarantee a unique solution. This process is called regularization. In general the cost function is modified by adding a regularization function $R : \setC^n \rightarrow \setR$ that depends on $\vec{x}$ $$ f(\vec{x}) = \norm{\vec{b} - \mat{A} \vec{x}}^2 + R(\vec{x}) . $$ The regularization function $R$ is somewhat arbitrary and many choices exist. Since components of non-trivial nullspaces can be arbitrarily large, it is reasonable to penalize the norm of solution vectors to prevent nullspaces from superfluously contributing to solution vectors. As with regularization functions, several norms exist.

An obvious choice is the norm induced by the standard complex inner product $$ \norm{\vec{x}}_2 = \aa{\vec{x}, \vec{x}}^{\q{1}{2}} = [\vec{x}^H \vec{x}]^{\q{1}{2}} . $$ This is the $L^2$ norm, and a subscript 2 is used here to distinguish it from other norms. Setting $R$ to the $L^2$ norm squared results in $\mathbf{L^2}$ (Tikhonov) regularization $$ R(\vec{x}) = \norm{\vec{x}}_2^2 . $$ Prior information about the solution can be incorporated into the $L^2$ regularizer by penalizing the difference between an expected solution vector $\vec{x}_p \in \setC^n$ and vector $\vec{x}$ $$ R(\vec{x}) = \norm{\vec{x}_p - \vec{x}}_2^2 . $$ Weighting can also be incorporated into the $L^2$ regularizer by multiplying the vectors in the norm by invertible matrix $\mat{\Gamma} \in \setC^{n \times n}$ $$ R(\vec{x}) = \norm{\mat{\Gamma} \bb{\vec{x}_p - \vec{x}}}_2^2 . $$ A common weighting choice is $\mat{\Gamma} = \gamma^{\q{1}{2}} \mat{1}$ where $\gamma \in \setR^{+}$ such that $$ R(\vec{x}) = \gamma \norm{\vec{x}_p - \vec{x}}_2^2 $$ which provides a way to adjust the relative importance of the residual norm versus the regularization function in the cost function. Unfortunately choosing the optimal $\gamma$ is problem-specific and usually done empirically [].

Other choices of regularization function exist. The generalization of the $L^2$ norm is the $L^p$ norm, which changes the shape of the unit sphere $$ \label{eqn:lpnorm} \mm{\vec{x}}_p \equiv \bigg[ \sum_{i=1}^n \nn{x_i}^p \bigg]^{\q{1}{p}} \quad\text{where}\quad p \in [1,\infty] \subset \setR . $$ Inserting $p = 0$ into Eq. $\eqref{eqn:lpnorm}$ results in the $L^0$ "norm" which counts the number of non-zero elements in a vector. $L^0$ is not a true norm (e.g. $\norm{\alpha \vec{x}}_0 \neq \nn{\alpha} \norm{\vec{x}}_0$), but it is an attractive regularizer for enforcing sparsity. Unfortunately $L^0$ minimization is NP-hard []. While it is difficult to use $L^0$ regularization directly, prior information can be used to sparsify a linear system []. More compuationally expensive regularization methods include gradient-based [] and total-variation-based [] schemes.

Left Preconditioning

To precondition a linear system means to transform it into a system with a better condition number before solving it. While this doesn't change the solution, it may make numeric computation more accurate and iterative algorithms converge faster. Left preconditioning can be incorporated into complex least squares by left multiplying the original linear system by invertible matrix $\mat{P} \in \setC^{m \cross m}$ $$ \mat{P} \mat{A} \vec{x} = \mat{P} \vec{b} $$ which results in a cost function similar to the weighted $L^2$ regularizer with prior $$ f(\vec{x}) = \mm{\vec{P} \bb{\vec{b} - \mat{A} \vec{x}}}^2 . $$ In the presence of measurement noise, the covariance matrix $$ \mat{K}_{\vec{x} \vec{y}} \equiv [\fn{cov}(x_i, y_j)] $$ can be used to weight measurement channels based on statistical considerations. Specifically, the inverse of the auto-covariance matrix, the precision matrix, of data vector $\vec{b}$ can be chosen as a preconditioner [] $$ \mat{P} = \mat{K}_{\vec{b}\vec{b}}^{-1} . $$

Preconditioning is best incorporated directly into iterative algorithms to take advantage of clever simplifications, as will be shown, but is included here for generality.

Solution Revisited

Regularization and preconditioning are now applied to formulate a well-posed least squares problem with a unique solution. The original complex least squares cost function Eq. $\eqref{eqn:cost}$ is the residual norm squared $$ f(\mathbf{x}) = \|\mathbf{b} - \mathbf{A} \mathbf{x}\|^2 . $$ This cost function is modified by incorporating left preconditioning and weighted $L^2$ regularization with a prior $$ f(\vec{x}) = \mm{\mat{P} \bb{\vec{b} - \mat{A} \vec{x}}}^2 + \mm{\mat{\Gamma} \bb{\vec{x}_p - \vec{x}}}^2 . $$ Optimize $f$ by taking the CR gradient using Eq. $\eqref{eqn:grad_cost}$ as a template and set to zero $$ \grad f(\vec{x}_s) = -\mat{A}^H \mat{P}^H \mat{P} \bb{\vec{b} - \mat{A} \vec{x}_s} - \mat{\Gamma}^H \mat{\Gamma} \bb{\vec{x}_p - \vec{x}_s} = \vec{0} $$ to derive the regularized complex least squares equation $$ \boxed{ \bb{\mat{A}^H \mat{P}^H \mat{P} \mat{A} + \mat{\Gamma}^H \mat{\Gamma}} \vec{x}_s = \bb{\mat{A}^H \mat{P}^H \mat{P} \vec{b} + \mat{\Gamma}^H \mat{\Gamma} \vec{x}_p} . } $$ While this looks complicated, this is equivalent to the modified linear system $$ \vec{A}' \vec{x}_s = \vec{b}' $$ where $$ \vec{A}' = \mat{A}^H \mat{P}^H \mat{P} \mat{A} + \mat{\Gamma}^H \mat{\Gamma} $$ and $$ \vec{b}' = \mat{A}^H \mat{P}^H \mat{P} \vec{b} + \mat{\Gamma}^H \mat{\Gamma} \vec{x}_p . $$ Matrix $\mat{A}'$ is conjugate-symmetric and positive-definite

Conjugate-symmetric: $$ \begin{equation*} \mat{A}^{\prime H} = [\mat{A}^H \mat{P}^H \mat{P} \mat{A} + \mat{\Gamma}^H \mat{\Gamma}]^H = \mat{A}^H \mat{P}^H \mat{P} \mat{A} + \mat{\Gamma}^H \mat{\Gamma} = \mat{A}' \end{equation*} $$ Positive-definite: $$ \begin{align*} \aa{\vec{x}, \mat{A}' \vec{x}} &= \aa{\vec{x}, [\mat{A}^H \mat{P}^H \mat{P} \mat{A} + \mat{\Gamma}^H \mat{\Gamma}] \vec{x}} \\ &= \aa{\vec{x}, \mat{A}^H \mat{P}^H \mat{P} \mat{A} \vec{x}} + \aa{\vec{x}, \mat{\Gamma}^H \mat{\Gamma} \vec{x}} \\ &= \aa{\mat{P} \mat{A} \vec{x}, \mat{P} \mat{A} \vec{x}} + \aa{\mat{\Gamma} \vec{x}, \mat{\Gamma} \vec{x}} \\ &= \mm{\mat{P} \mat{A} \vec{x}}^2 + \mm{\mat{\Gamma} \vec{x}}^2 \\ &\geq \mm{\mat{\Gamma} \vec{x}}^2 \\ &> 0 \quad\text{when}\quad \vec{x} \neq \vec{0} \end{align*} $$

and thus invertible $$ \vec{x}_s = \mat{A}^{\prime -1} \vec{b} $$ because a square matrix with a trivial nullspace is full rank by the rank-nullity theorem.

For the special case of $\mat{P} = \mat{1}$, $\mat{\Gamma} = \gamma^{\q{1}{2}} \mat{1}$, and $\vec{x}_p = \vec{0}$, the modified complex least squares equation reduces to the form $$ \boxed{ \bb{\mat{A}^H \mat{A} + \gamma \mat{1}} \vec{x}_s = \mat{A}^H \vec{b} } . $$ This is a well-posed inverse problem with a unique solution!

Summary

TODO

References