How to use lsqnonlin on a function which is mostly linear?

I have a (large) least squares minimization problem which is almost linear.

I want to find a set of weights for a small (7x7) kernel which best filters noise from a greyscale image. I have a reference no-noise image which has spatial correlation from both lens blur and structure in the scene. I also know how much measurement noise (approximately gaussian) will be added to each pixel.

My approach is to produce a matrix A with hundreds of thousands of rows (n), one for each pixel, and 49 columns, one for each pixel in the neighborhood of a filtered pixel. A contains the reference no-noise image data. I also have the desired output Y, which has the hundreds of thousands of rows (one for each pixel) and a single column (the center pixel), also no-noise image data. I'm looking for a 49x1 weight vector x.

When taking a picture, every pixel will have error added with variance V to each pixel.

So I'm looking for the x that least squares minimizes sum((A*x - Y).^2) + n * x'*x * V

The trivial filter weights are zero for everything but the center pixel. These give A*x-Y=0, and sees n*V summed variance across the entire picture. I use this as my starting point. I have the optimization working now with lsqnonlin, in which the optimizer figures out how much to weight the surrounding pixels to trade reduction in variance for an increase in estimation inaccuracy.

Yes, there are linear constraints on the weights (must sum to 1, center of mass must be centered). I actually have fewer than 49 optimization variables that I transform into the weights to ensure these constraints, but (a) that part works, and (b) it's an unnecessary complication to the problem I'd like help with.

The trouble is that my objective function is giving lsqnonlin a gigantic result with n+1 terms (one extra term is sqrt(n*x'*x*V)), so that the Jacobian is enormous (40-odd optimization variables by hundreds of thousands of rows). The optimizer (trust-region-reflective) is stopping after it exceeds the function evaluation limit (97 iterations, 4400 evaluations), while it still has a largish stepsize, so I don't think I've found a good minimum yet. It's also slow and I'd like to run thousands of these and expand the number of rows a hundredfold.

I think my approach is very inefficient and there is probably a better way. What is it?

Answer 1

[Matt J]

I think you should just use quadprog. No Jacobians involved and it will be easier to add in your linear constraints later on.

H=A.'*A + n*V*eye(49); %fixed
f=-(Y.'*A).';

x = quadprog(H,f);

[Iain McClatchie]

Thanks for the quick response! I'm verifying that this works now. Can you just explain the difference between

A.' * A

and

A' * A

At the command prompt, they appear to be identical.

[Matt J]

A' is conjugate tranpose

https://www.mathworks.com/help/matlab/ref/ctranspose.html

while A.' is the ordinary matrix transpose. They will be identical for real A.

[Iain McClatchie]

It works! Thanks, Matt.

I have two minor questions about how you derived your H and f.

The expression I wanted to optimize was

sum((A*x - Y).^2) + n * x'*x * V

If I transform that into what quadprog() wants, which is what you did, I get

=(A*x - Y).' * (A*x - Y) + n * x.' * x * V
=x.'*A.'*A*x - 2*Y.'*A*x +Y.'*Y + x.'*n*V*eye(49)*x
=x.'*(A.'*A + n*V*eye(49))*x - 2*Y.'*A*x + Y.'*Y

Now Y.'*Y is a constant, and since we're just looking for the minimum we can throw that out, and divide by 2 to make it look like what quadprog() wants.

min 0.5*x.'*H*x + f.'*x

so

H = A.'*A + n*V*eye(49);
f = -A.'*Y;

This isn't exactly what you suggested.

You had suggested I needed to multiply A.'*A by 0.5. I don't see why that is true. Can you show me what I did wrong?

Second, my f is the transpose of your, since the documentation for quadprog() says it wants a column vector. I tried it with a row vector and got the same answer, so your suggestion works, but if I'm reading correctly my formulation is consistent with the documentation.

[Matt J]

Yes, I think your derivation is the correct one. I obviously had one or two mistakes.

[Iain McClatchie]

Great. Now I believe I understand what I'm doing. Thanks for the help. And you were right, dropping the linear constraints in was a snap, way easier than the transform I was doing with lsqnonlin().