How to perform quadratic optimization

Hi, I am trying to perform quadratic optimisation for correction the Edge spread function as shown in the figure.Looking at the figure we can visualize optimize solution that the ESF should not have bumps before start and end of slope but want to optimize using optimization methods.

I never used the optimisation tool before. I want to optimize it by minimizing the equation as shown below and constrain.

I tried following thing but it shooting me error fmincon stopped because it exceeded the function evaluation limit

if true
  % code
A = zeros(length(ESF),length(ESF));
% ind_c = 0;
for i = 1:length(ESF)-1
 A(i,i) = 1;
 A(i,i+1)= -1;
end
X = sym('x',[365,1]);

func = @(y)(sum((y-ESF).^2)); b = zeros(size(A,1),1); y = fmincon(func,Es',A,b); end

i will appreciate your help. Thank you

Answer 1

[John D'Errico]

I'll assume that you have the optimization toolbox. If you don't, you could still do this, but I'm not going to teach you how to write the equivalent of lsqlin, and you need lsqlin to do this efficiently. Ok, yes, you can solve the problem using lsqnonneg.

I assume that you have several hundred points from a curve, and you wish to solve the problem where you minimize the difference between the data and an approximate curve, subject to the constraint that the new curve is truly monotone. I don't see any data attached, so I cannot give you an example where I solve your problem.

The absolutely simplest solution is to use my SLM toolbox, to fit a (piecewise linear) spline through your data, subject to the constraints that the curve is monotone.

(I need to run right now for about an hour. I'll return and see if I can cook up an example of how that will work. If you are able to attach some data, that would help.)

More complex is you could solve it using lsqlin. But that will force you to formulate the problem for that tool. Slightly harder than using SLM, but not a lot harder. BRB...

Back, with an example of what I THINK you are asking to do. I'll start with a smooth curve over a few points in x. I'll use erf, to create a basic sigmoidal shape. Then interpolate it using spline. This will incur ringing in the curve.

x = linspace(-15,15,20);
xi = -15:.1:15;
yi = spline(x,erf(x),xi);
plot(xi,yi)

So you can see the ringing/edge effects in the curve.

slm = slmengine(xi,yi,'knots',xi,'degree',1,'plot','on','increasing','on');

The plot with the points and the fitted curve is hard to see, but I'll zoom into the shoulder area next.

As you see, the curve follows the original data exactly wherever it was monotone already. It only modifies the shape when it needs to do so, using the requirement of monotonicity.

You can extract the actual fitted y values (for this curve type) directly from the struct it produces as:

yihat = slm.coef';

So this is a very simple way to solve your problem. Simple, because slmengine did all the work for you. Find the SLM toolbox on the file exchange at this link .

I could also have approximated the curve using a nice smooth spline, again using monotonicity as a requirement. But that was not your stated goal.

And yes, I can show you how to solve this problem using other tools, like lsqlin or lsqnonneg. But do I really need to do so? I am being lazy and hoping SLM will suffice. :)

[Parag]

Thank you so much for you help. I am attching the mat file here and trying to impliment using Isqlin. I tried to solve using fmincon. Here is some code i tried. I really appreciate you help.

if true
  A = zeros(length(ESF),length(ESF));
% ind_c = 0;
for i = 1:length(ESF)-1
 A(i,i) = 1;
 A(i,i+1)= -1;
end
%%X = sym('x',[365,1]);

func = @(y)(sum((y-ESF).^2)); b = zeros(size(A,1),1); y = fmincon(func,ESF,A,b'); end

[John D'Errico]

DON'T USE FMINUNC! Why would you use a nonlinear tool for an essentially linear problem? Using fmincon here is the equivalent to the use of a Mack truck to take a single pea to Boston.

Anyway, I have no idea why you want to write it yourself since I showed you how to solve it in one line already using SLM. I suppose if your goal is to understand how to write it, oh well...

I'll assume that you have a list of points as (x,y) pairs, so two vectors of the same size. I'll use the data you supplied instead of the data I cooked up before.

y = ESF;
n = numel(y);
x = (1:n)';

plot(x,y)
grid on

As a problem for lsqlin, the "objective" is a simple one. lsqlin solves a linear least squares problem. Our unknowns here will be the vector yhat.

The objective for lsqlin will be essentially a set of linear equations of the form C*yhat=y. Here C is a simple matrix to build

C = eye(n,n);

The constraints on yhat are also trivial to write. They will simply be a b-diagonal matrix, of size n-1 by n. Create that matrix by subtracting the two diagonals. I could also have written C in a call to toeplitz, but this was trivial as is.

Aineq = [eye(n-1),zeros(n-1,1)] - [zeros(n-1,1),eye(n-1)];
bineq = zeros(n-1,1);

Solve, using lsqlin. Set a few options first, and give it what we know is a good place to start, thus yhat=y initially.

options = optimset('lsqlin');
options.Display = 'none';
options.Algorithm = 'active-set';
[yhat,resnorm,residual,exitflag] = lsqlin(C,y,Aineq,bineq,[],[],[],[],y,options);

plot(x,y,'b-',x,yhat,'r-')
grid on

The solution replicates y exactly in the transition region, and is constant, or at worst, monotonically increasing in the tails. Since your data was just white noise in the tails, yhat will be essentially constant there.

Note that this will be essentially the same solution that my SLM toolbox would have generated.

[Parag]

Thank you so much for your help and brief explanation. I learned a lot

Answer 2

[Teja Muppirala]

Here is how you might solve it using QUADPROG. Note that you don't need to do anything symbolic using "sym".

L = length(ESF);

H = eye(L);
f = -ESF;

A = zeros(L-1,L);
for i = 1:L-1
    A(i,i) = 1;
    A(i,i+1)= -1;
end
b = zeros(L-1,1);

yopt = quadprog(H,f,A,b);

figure
plot(ESF);
hold on;
plot(yopt);

[Parag]

Thank you so much for your help. It works really well by using quadprog

Answer 3

[Jack]

The physical Edge Spread Function is monotonic in normal conditions, so the undershoots/overshoots are due to processing (demosaicing or, more likely, sharpening). If one uses raw data, there never are under/overshoots.

This is actual raw data projected onto the normal, from the green channels of a 310x248px edge capture, slanted 13.45 degrees off the vertical:

We would like to get an ESF out of that without doing much filtering, if any. Quadprod does a decent job of it with Teja's code (1 and -1 switched in the for loop for this orientation, correct?) compared to a fit by the sum of three sigmoids:

There is a little more nuance in quadprod's ESF and that nuance is important for higher frequency accuracy in the OTF.

Thank you Teja.

Jack

PS Projected raw data attached for anyone who would like to try it or improve upon this (pr is location on the projected axis, esfr is the raw intensity at pr, mid is an educated guess at the center of the edge.