Interpolation Optimization Problem
Using few tricks the equation can be transformed into classic Quadratic Programming problem:

$ E = \arg \max_{E} \sum_{\mathbf{r}} \left( E(\mathbf{r}) - \sum_{\mathbf{s} \in N(\mathbf{r})} {w}_{\mathbf{rs}} E(\mathbf{s}) \right)^2 $

Defining $ \mathbf{e} = vec \left ( E \right ) $, Namely vectorizing the matrix, and defining a matrix $ A $ defined the proper usage of the weights $ {w}_{\mathbf{rs}} $, Namely put them in the correct entry at each row per pixel the above could be rewritten as:

$$ \mathbf{e} = \arg \max_{\mathbf{e}} {\left ( \mathbf{e} - A\mathbf{e} \right ) }^{T} \left ( \mathbf{e} - A\mathbf{e} \right ) = \arg \max_{\mathbf{e}} {\left( \left ( I - A \right ) \mathbf{e} \right )}^{T} \left( \left ( I - A \right ) \mathbf{e} \right ) $$

Defining $ B = I - A $ yields:

$$ \mathbf{e} = \arg \max_{\mathbf{e}} {\left( \left ( I - A \right ) \mathbf{e} \right )}^{T} \left( \left ( I - A \right ) \mathbf{e} \right ) = \arg \max_{\mathbf{e}} {\mathbf{e}}^{T} {B}^{T} B \mathbf{e} $$

Defining $ L = {B}^{T} B $ , namely $ L , is a Symmetric Positive Definite Matrix yields:

$$ \mathbf{e} = \arg \max_{\mathbf{e}} {\mathbf{e}}^{T} L \mathbf{e}  $$

Adding the equality constraints using $ C \mathbf{e} = \mathbf{t} $ where $ \mathbf{t} $ is vectorization of all "Reference Pixels"

All in all this yields a Classic Quadratic Programming problem given by:

$$ \mathbf{e} = \arg \max_{\mathbf{e}} {\mathbf{e}}^{T} L \mathbf{e}, \; s.t. \; C \mathbf{e} = \mathbf{t} $$

Few notes:


The matrix $ L $ is Symmetric and Positive Definite matrix and usually sparse (Assuming "Small Neighborhood").
The matrix $ L $ wast built by the Matrix $ I - A $ which each row of it was normalized to be zero (Each row of $ A $ is normalized to 1 ).


Now two things are needed:


Optimized fast algorithm to solve this problem utilizing its properties.
Assuming the image is large, are the methods to solve this problem by blocks to avoid memory limits?