Constrained Weighted Least Squares (Image Interpolation)

Hello,

I'd like your assistance with solving a constrained weighted least squares problem.

I can't use `lsqlin` since the size of the images creates matrices beyond the abilities of my computer (Or any PC for that matter).

We're given a matrix I of M by N. On this matrix set of L pixels S is set to be anchors. All the other pixels should be recalculated based on that pixels using the following constrains:

Which is equivalent of:

Where the matrix A just "Selects" the anchor pixels from the interpolated vector.

Now, the trick for most pixels, outside the relevant pixel neighborhood, the weights are zeros. Namely, the equation above could be written as:

Now, assume the neighborhood is something like K x K neighborhood. I think this should allow moving into a sparse problem.

Assuming we have all the W_rs needed, how can I built this problem in a manner MATLAB's can solve it?

Thank You.

Giving an example for w_rs:

Where I is the original values of the pixels at pixel r and s respectively.

Of course it is normalized per pixel.

[Matt J]

Is W_rs separable in r and s? In other words, can it be decomposed as

 W_rs= w_r * w_s

Or might it have another kind of separability, e.g., in x and y similar to

 W(x,y,x',y')=w(x-x')*w(y-y')

[Matt J]

What about shift-invariance?

 W_rs=W_(r-s)

[Royi Avital]

Hi Matt,

Yes, W_rs = W_sr.

Yet it is spatially variant.

I added a definition (One of many) at the question itself.

Thank You.

[Matt J]

OK. Just to be clear, I didn't ask if W_rs=W_sr. I asked if it has the separable form

 W_rs= f_s *f_r

Your recently posted example for W_rs does not have this form.

Answer 1

[Matt J]

Here's what I'd propose. Note that your constraints are highly trivial and shouldn't require an iterative solver. Just treat the E(p)=I(p) as knowns (which is what they are) and apply the solver exclusively to the unknowns E(~p).

    %%Fake initial data
    M=200;
    N=M;
    K=3;
    I=rand(M,N);
    sig=1;
    p=false(1,M*N); 
      p(1:10:end)=1;

   tic;  %%Set up equation coefficient matrix

    xx=1:M;
    yy=1:N;

    delta=floor((-K:K)/2); %assume K odd

    [x,xdelta,y,ydelta]=ndgrid(xx,delta,yy,delta);

    xd=x+xdelta;
    yd=y+ydelta;
    idx=(xd>=1 & xd<=M & yd>=1 & yd<=N); %legitimate indices

    x=x(idx); xd=xd(idx);
    y=y(idx); yd=yd(idx);

    r=sub2ind([M,N],x,y); r=r(:);
    s=sub2ind([M,N],xd,yd); s=s(:);

    W=spdiags(exp(-(I(r)-I(s))/2/sig^2),0, length(r),length(r));

     A=W*(sparse(1:length(r),r,1)-sparse(1:length(s),s,1));

   toc;

   tic; %%Solve equations

     E=zeros(M*N,1);
     E(p)=I(p);
     E(~p) = -A(:,~p)\(A*E);  %solve

   toc;

[Royi Avital]

Hi Matt, Do you think MATLAB would be able to handle images with resolution of [4000 x 4000]?

Thank You.

[Matt J]

Not sure. It partly depends on K.

[Royi Avital]

Let's say K is 7.

Moreover, how would you change your solution if the function would be as described here:

http://www.mathworks.com/matlabcentral/answers/125746-solution-for-an-optimization-cost-function-interpolation-optimization

Thank You.

[Matt J]

Well, assuming my code works as I think it does, it generates vectors r,s containing all the pairs r(i) and s(i) for which W(r(i),s(i))~=0. That's all that is required for both problem formulations.

I think N=4000 will be too big, however.

[Royi Avital]

Hi Matt,

Is there a chance that:

delta=floor((-K:K)/2); %assume K odd

Should be:

delta = -floor(K / 2):floor(K / 2); %assume K odd

Thanks.

[Royi Avital]

Hi Matt,

I think you forgot the power of 2 in the exponent.

[Matt J]

Yep. You're corrections are correct in both places.

I also never did the normalization sum_s Wrs=1

[Royi Avital]

Since I don't really understand how did you achieve this solution, How can I incorporate the normalization?

Thank You for your assistance.

Moreover, How come this solution solves the other problem I linked to?

[Matt J]

I think the best thing would be to compute W_rs as a sparse matrix

 W=sparse(r,s,exp(-(I(r)-I(s).^2)/2/sig^2,  N*M,N*M);

Now you can just normalize the rows to sum to 1

 sW=sum(W,2);
  nzRows=sW>0;

 W(nzRows,:)=bsxfun(@rdivide,W(nzRows,:),sW(nzrows));

and now

 [r,s,w]=find(W);

Although, I think that normally all(nzRows)=true.

Note that once you have the data triples r,s,w and/or the matrix W, then either of your objective functions can be expressed almost as written and fed to fmincon,e.g.,

function f=objfun1(E,r,s,w)

    f=sum(w.*(E(r)-E(s)).^2);

end

function f=objfun2(E,W)

     f=norm(E-W*E).^2;

end

Answer 2

[Matt J]

Maybe a more transparent way to see the problem is to notice that it can be expressed in the equivalent matrix/vector form,

 min E.'*(I-W)*E  (*)

ignoring your constraints for a moment. To see the above, take the given form in your post (divided by two) and do the following massaging, where W is the sparse matrix with components W(r,s).

 f/2=sum_rs(w_rs.*(E(r)-E(s)).^2)/2;
    =sum_rs(w_rs.*(E(r)^2)/2 +sum_rs(w_rs*E(s).^2)/2-sum(w_rs*E(r)*E(s))
    =norm(E).^2 - E.'*W*E
    =E.'(I-W)*E

which is exactly (*). Thus, once you've generated W using the sparse command, everything is expressible as a simple quadratic program.

Similarly the other version of the problem can be expressed

 min. f=norm(E-W*E)^2

or equivalently,

 min E*((I-W)^2)*E

So, the only difference between the two problems is that the Hessian 2*(I-W) in the first gets squared in the second 2*(I-W)^2.

To deal with "constraints" in either problem, write the problem generically as

 min f(E) = E.'*Q*E/2  (**)

and rewrite E in the partitioned form E=[Eu;Ep] where Eu are the unknown E(r) and Ep are the known E(p). Similarly, parition Q as

   Q=[Quu, Qup;Qup.' Qpp]

Substituting into () will lead to an unconstrained quadratic program in Eu

 min f(Eu) = Eu.'*Quu*Eu/2 -(Qup*Ep).'*Eu

whcih can be minimized analytically

   Eu=Quu\(Qup*Ep)

How computationally hard this all is, depends on the sparsity of Q and hence the sparsity of W. The number of bytes needed to hold the sparse matrix W will be approximately 24*N*M*K^2.

For N-M=4000 and K=7, this comes out to about 17GB. So, you'll need to cut back somewhere, I think.