How to perform quadratic optimization

5 次查看（过去 30 天）
Parag
Parag 2016-4-27
回答： Jack 2019-10-1
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

请先登录，再回答此问题。

采纳的回答

John D'Errico
John D'Errico 2016-4-27
编辑：John D'Errico 2016-4-27
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. :)
  3 个评论
显示 1更早的评论
John D'Errico
John D'Errico 2016-4-28
编辑：John D'Errico 2016-4-28
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
Parag 2016-4-28
Thank you so much for your help and brief explanation. I learned a lot

请先登录，再进行评论。

更多回答（2 个）

Teja Muppirala
Teja Muppirala 2016-4-28
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);
  1 个评论
Parag
Parag 2016-4-28
Thank you so much for your help. It works really well by using quadprog

请先登录，再进行评论。

Jack
Jack 2019-10-1
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:
Esge Profile Projected on Normal GFX50s110mmf28.png
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:
Binned Edge Profile Sigmoids vs Quad Prog Monotonic_.png
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.

类别

在 Help Center 和 File Exchange 中查找有关 Linear Least Squares 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by