最小化二次成本函数：检查算法的有效性以及如何优化它
Minimazing quadratic cost function: check validity of algorithm and how to optimize it
我想在给定输出数据“y”的情况下估计 PSF 的参数和线性系数。关系由下式给出：
我有以下功能
我必须估计仅给出“y”输出数据的所有 6 个参数（a、b、r0、c0、alpha、beta）。
为此，我编写了以下成本函数以最小化：
function [J, Model] = Crit_J(p,D)
% Compute the model corresponding to parameters p
[R,C] = size(D);
[Cols,Rows] = meshgrid(1:C,1:R);
% Model
Model = p(1)*(1+((Rows-p(3)).^2+(Cols-p(4)).^2)/p(5)^2).^(-p(6)) + ...
p(2)*ones(R,C);
model = Model(:);
% Introduce H matrix
d = D(:);
H = [ model, ones(length(model),1)];
% Compute the cost function : taking absolute value
J = abs((d-H*[p(1),p(2)]')'*(d-H*[p(1),p(2)]'));
end
向量 p :
% Parameters theta and nu: initial values
a = 10;
b = 3;
r0 = 0;
c0 = 0;
alpha = 1;
beta = 1;
% Array of parameters
p = [a,b,r0,c0,alpha,beta]';
我必须明确矩阵 H 等于（R=C=20）：
在主程序中，我是这样使用 Crit_J 函数的：
% Initial value pour J cost function
Jmin = Inf;
% Initial value for p vector
init = 10;
% Minimization
%for i=1:1000000
kmax = 1e6
epsilon = 100;
k = 0;
while (k < kmax && Jmin > epsilon)
p = init*randn(6,1);
J = Crit_J(p,D);
if (J < Jmin)
Jmin = J;
pCurrent = p;
end
k=k+1;
end
不幸的是，即使迭代次数很多，估计也没有得到很好的估计。
我不知道我必须为参数取什么初始值（在循环之前）。
如果我使用迭代：while (Jmin > epsilon)，程序永远不会 returns 因为 epsilon 值可能太小了。我更喜欢使用经典的 for i=1:1000000.
这就是为什么我想检查最小化算法是否正确的原因。
如何确认函数 Crit_J 已在主程序中正确定义和使用？
是否存在实现这种最小化的 Matlab 函数？
更新 1
我遵循 @inkert 给出的答案并使用 Crit_J2 函数，如下所示：
function cost = Crit_J2(p,D,alpha,beta,r0,c0)
% Compute the model corresponding to parameters p
[R,C] = size(D);
[Cols,Rows] = meshgrid(1:C,1:R);
% Model
Model = p(1)*(1+((Rows-r0).^2+(Cols-c0).^2)/alpha^2).^(-beta) + ...
p(2)*ones(R,C);
%Model = p(1)*D + p(2)*ones(R,C);
model = Model(:);
d = D(:);
%disp(model);
% Introduce H matrix
H = [ model, ones(length(model),1)];
% Compute the cost function : taking absolute value
J2 = abs((d-H*[p(1),p(2)]')'*(d-H*[p(1),p(2)]'));
cost = J2;
end
不幸的是，即使我在代码开始时使用接近于固定的初始参数，估计也很糟糕。
我想知道我的二次成本函数是否有效；事实上，我使用了以下表达式：
J = abs((d-H*[p(1),p(2)]')'*(d-H*[p(1),p(2)]'));
d(=D(:))是输出数据向量，H是矩阵允许写成矩阵形式：
d = H.p + epsilon
我是否必须在 J 成本函数的定义中添加噪声？
您可以在 : example output.dat
上获取二维输出文件
我得到：
p_estimate =
1.1610
2.3271
13.2684
13.6551
3.3586
0.4541
and expected : 14
5
8.56
10.995
3
2.45
此处为说明图片：
如您所见，我的代码中某处存在问题。
更新 2
Matlab 的函数“ga”也会产生错误的估计。我是这样用的:
% Linear inequalities : A*x < b
A1 = [0, 0, -1, 0, 0, 0];
b1 = 0;
A2 = [0, 0, 0, -1, 0, 0];
b2 = 0;
A3 = [0, 0, 0, 0, -1, 0];
b3 = 0;
A4 = [0, 0, 0, 0, 0, -1];
b4 = -1;
%% Matrix A and b
A = [A1; A2; A3; A4];
b = [b1; b2; b3; b4];
%
% 2 lines below : attempt to find minimum
[pars, Jmin] = ga(@(x)Crit_J(x,D),6,A,b,[],[],[],[],[]);
% Function of cost
function cost = Crit_J(p,D)
% Compute the model corresponding to parameters p
[R,C] = size(D);
[Cols,Rows] = meshgrid(1:C,1:R);
% Model
Model = p(1)*(1+((Rows-p(3)).^2+(Cols-p(4)).^2)/p(5)^2).^(-p(6)) + ...
p(2)*ones(R,C);
model = Model(:);
% Introduce H matrix
d = D(:);
H = [ model, ones(length(model),1)];
% Compute the cost function : taking absolute value
J = abs((d-H*[p(1),p(2)]')'*(d-H*[p(1),p(2)]'));
cost = J;
end
不等式如下：
p = [a,b,r0,c0,alpha,beta];
a=p(1) > 0
b=p(2) > 0
r0=p(3) > 0
c0=p(4) > 0
alpha=p(5) > 0
beta=p(6) > 1
任何人都可能看到这些不等式或使用“ga”的错误？
更新 3
这是我关于最小化成本函数的最后结果：好消息和坏消息。
1) 好消息： 现在，3 种方法（fminsearch、fminunc 和 ga）给出了相同的估计2 个随机参数（a 和 b 到 y = a*PSF+b）。
2) 坏消息： a 和 b 的估计值太小了。
现在，我将 3 个方法（即上面引用的 3 个函数）的小代码放入：
% Amplitude value = a
amplitude = 14;
% Sky background = b
backgvalue = 5;
backg = backgvalue*ones(L,C);
% White Noise sigma
variance = 1;
WhiteNoise = sqrt(variance)*randn(L,C);
% Generate grid
[Cols,Rows] = meshgrid(1:C,1:L);
% Pack parameters
nu = [r0; c0; alpha; beta];
% Create normalizded PSF with amplitude = 1 and without background and %gaussian white noise
PSF = Moffat(nu,Rows,Cols);
% Normalized PSF and general case
Dnorm = PSF;
D = amplitude*PSF + backg + WhiteNoise;
a = 3;
b = 1;
p0 = [a,b];
%%%%%%%%%%%% fminsearch %%%%%%%%%%%%%%%%%%
pars = fminsearch(@(x)Crit_J2(x,D,alpha,beta,r0,c0),p0);
disp('pars : fminsearch');
pars
%%%%%%%%%%%% fminunc %%%%%%%%%%%%%%%%%%
[pars, Jmin] = fminunc(@(x)Crit_J2(x,D,alpha,beta,r0,c0),p0);
disp('pars : fminunc');
pars
%%%%%%%%%%%% ga %%%%%%%%%%%%%%%%%%
%%% impose a > 0 and b > 0
A1lin = [-1, 0];
b1lin = 0;
A2lin = [0, -1];
b2lin = 0;
%%% Matrix A and b
Alin = [A1lin; A2lin];
blin = [b1lin; b2lin];
[pars, Jmin] = ga(@(x)Crit_J2(x,D,alpha,beta,r0,c0),2,Alin,blin);
disp('pars : ga');
pars
估计是：a = 3.74和b = 1.04（而不是起始固定值a = 14和b = 5）
我是否应该断定问题来自我的代码中使用的成本函数的定义？即：
% Function of cost
function cost = Crit_J2(p,D,alpha,beta,r0,c0)
% Compute the model corresponding to parameters p
[R,C] = size(D);
[Cols,Rows] = meshgrid(1:C,1:R);
% Model
Model = p(1)*(1+((Rows-r0).^2+(Cols-c0).^2)/alpha^2).^(-beta) + ...
p(2)*ones(R,C);
%Model = p(1)*D + p(2)*ones(R,C);
model = Model(:);
d = D(:);
%disp(model);
% Introduce H matrix
H = [ model, ones(length(model),1)];
% Compute the cost function : taking absolute value
J2 = abs((d-H*[p(1),p(2)]')'*(d-H*[p(1),p(2)]'));
cost = J2;
end
我正在使用以下成本函数：
J = abs((d-H*[p(1),p(2)]')'*(d-H*[p(1),p(2)]'));
withd估计开始时生成的初始数据向量，H矩阵使得d = H.thetawiththeta=[a,b]=[p(1),p(2)]'
[a,b]（恒星振幅+背景）这个错误估计（under estimation）的解释是什么？
现在您只是随机更改参数。如果您更改为 fminsearch 中实现的优化程序，此算法将为您找到最小值。
% initial parameters
a = 10;
b = 3;
r0 = 0;
c0 = 0;
alpha = 1;
beta = 1;
% Array of parameters
p0 = [a,b,r0,c0,alpha,beta]';
% Minimization of cost function for p
[pars, Jmin] = fminsearch(@(x)Crit_J(x,D), p0);
您必须将 Crit_J 函数更改为仅输出 'cost' 以最小化此功能。然后您可以稍后通过 运行 原始函数获得您的模型。
更新：
也许您正在寻找局部最小值，您可以尝试 fminunc instead of fminsearch, or switch to a generic algorithm，但这对于这个问题来说可能有点矫枉过正。
我想在给定输出数据“y”的情况下估计 PSF 的参数和线性系数。关系由下式给出：
我有以下功能
我必须估计仅给出“y”输出数据的所有 6 个参数（a、b、r0、c0、alpha、beta）。
为此，我编写了以下成本函数以最小化：
function [J, Model] = Crit_J(p,D)
% Compute the model corresponding to parameters p
[R,C] = size(D);
[Cols,Rows] = meshgrid(1:C,1:R);
% Model
Model = p(1)*(1+((Rows-p(3)).^2+(Cols-p(4)).^2)/p(5)^2).^(-p(6)) + ...
p(2)*ones(R,C);
model = Model(:);
% Introduce H matrix
d = D(:);
H = [ model, ones(length(model),1)];
% Compute the cost function : taking absolute value
J = abs((d-H*[p(1),p(2)]')'*(d-H*[p(1),p(2)]'));
end
向量 p :
% Parameters theta and nu: initial values
a = 10;
b = 3;
r0 = 0;
c0 = 0;
alpha = 1;
beta = 1;
% Array of parameters
p = [a,b,r0,c0,alpha,beta]';
我必须明确矩阵 H 等于（R=C=20）：
在主程序中，我是这样使用 Crit_J 函数的：
% Initial value pour J cost function
Jmin = Inf;
% Initial value for p vector
init = 10;
% Minimization
%for i=1:1000000
kmax = 1e6
epsilon = 100;
k = 0;
while (k < kmax && Jmin > epsilon)
p = init*randn(6,1);
J = Crit_J(p,D);
if (J < Jmin)
Jmin = J;
pCurrent = p;
end
k=k+1;
end
不幸的是，即使迭代次数很多，估计也没有得到很好的估计。
我不知道我必须为参数取什么初始值（在循环之前）。
如果我使用迭代：while (Jmin > epsilon)，程序永远不会 returns 因为 epsilon 值可能太小了。我更喜欢使用经典的 for i=1:1000000.
这就是为什么我想检查最小化算法是否正确的原因。
如何确认函数 Crit_J 已在主程序中正确定义和使用？
是否存在实现这种最小化的 Matlab 函数？
更新 1
我遵循 @inkert 给出的答案并使用 Crit_J2 函数，如下所示：
function cost = Crit_J2(p,D,alpha,beta,r0,c0)
% Compute the model corresponding to parameters p
[R,C] = size(D);
[Cols,Rows] = meshgrid(1:C,1:R);
% Model
Model = p(1)*(1+((Rows-r0).^2+(Cols-c0).^2)/alpha^2).^(-beta) + ...
p(2)*ones(R,C);
%Model = p(1)*D + p(2)*ones(R,C);
model = Model(:);
d = D(:);
%disp(model);
% Introduce H matrix
H = [ model, ones(length(model),1)];
% Compute the cost function : taking absolute value
J2 = abs((d-H*[p(1),p(2)]')'*(d-H*[p(1),p(2)]'));
cost = J2;
end
不幸的是，即使我在代码开始时使用接近于固定的初始参数，估计也很糟糕。
我想知道我的二次成本函数是否有效；事实上，我使用了以下表达式：
J = abs((d-H*[p(1),p(2)]')'*(d-H*[p(1),p(2)]'));
d(=D(:))是输出数据向量，H是矩阵允许写成矩阵形式：
d = H.p + epsilon
我是否必须在 J 成本函数的定义中添加噪声？
您可以在 : example output.dat
上获取二维输出文件我得到：
p_estimate =
1.1610
2.3271
13.2684
13.6551
3.3586
0.4541
and expected : 14
5
8.56
10.995
3
2.45
此处为说明图片：
如您所见，我的代码中某处存在问题。
更新 2
Matlab 的函数“ga”也会产生错误的估计。我是这样用的:
% Linear inequalities : A*x < b
A1 = [0, 0, -1, 0, 0, 0];
b1 = 0;
A2 = [0, 0, 0, -1, 0, 0];
b2 = 0;
A3 = [0, 0, 0, 0, -1, 0];
b3 = 0;
A4 = [0, 0, 0, 0, 0, -1];
b4 = -1;
%% Matrix A and b
A = [A1; A2; A3; A4];
b = [b1; b2; b3; b4];
%
% 2 lines below : attempt to find minimum
[pars, Jmin] = ga(@(x)Crit_J(x,D),6,A,b,[],[],[],[],[]);
% Function of cost
function cost = Crit_J(p,D)
% Compute the model corresponding to parameters p
[R,C] = size(D);
[Cols,Rows] = meshgrid(1:C,1:R);
% Model
Model = p(1)*(1+((Rows-p(3)).^2+(Cols-p(4)).^2)/p(5)^2).^(-p(6)) + ...
p(2)*ones(R,C);
model = Model(:);
% Introduce H matrix
d = D(:);
H = [ model, ones(length(model),1)];
% Compute the cost function : taking absolute value
J = abs((d-H*[p(1),p(2)]')'*(d-H*[p(1),p(2)]'));
cost = J;
end
不等式如下：
p = [a,b,r0,c0,alpha,beta];
a=p(1) > 0
b=p(2) > 0
r0=p(3) > 0
c0=p(4) > 0
alpha=p(5) > 0
beta=p(6) > 1
任何人都可能看到这些不等式或使用“ga”的错误？
更新 3
这是我关于最小化成本函数的最后结果：好消息和坏消息。
1) 好消息： 现在，3 种方法（fminsearch、fminunc 和 ga）给出了相同的估计2 个随机参数（a 和 b 到 y = a*PSF+b）。
2) 坏消息： a 和 b 的估计值太小了。
现在，我将 3 个方法（即上面引用的 3 个函数）的小代码放入：
% Amplitude value = a
amplitude = 14;
% Sky background = b
backgvalue = 5;
backg = backgvalue*ones(L,C);
% White Noise sigma
variance = 1;
WhiteNoise = sqrt(variance)*randn(L,C);
% Generate grid
[Cols,Rows] = meshgrid(1:C,1:L);
% Pack parameters
nu = [r0; c0; alpha; beta];
% Create normalizded PSF with amplitude = 1 and without background and %gaussian white noise
PSF = Moffat(nu,Rows,Cols);
% Normalized PSF and general case
Dnorm = PSF;
D = amplitude*PSF + backg + WhiteNoise;
a = 3;
b = 1;
p0 = [a,b];
%%%%%%%%%%%% fminsearch %%%%%%%%%%%%%%%%%%
pars = fminsearch(@(x)Crit_J2(x,D,alpha,beta,r0,c0),p0);
disp('pars : fminsearch');
pars
%%%%%%%%%%%% fminunc %%%%%%%%%%%%%%%%%%
[pars, Jmin] = fminunc(@(x)Crit_J2(x,D,alpha,beta,r0,c0),p0);
disp('pars : fminunc');
pars
%%%%%%%%%%%% ga %%%%%%%%%%%%%%%%%%
%%% impose a > 0 and b > 0
A1lin = [-1, 0];
b1lin = 0;
A2lin = [0, -1];
b2lin = 0;
%%% Matrix A and b
Alin = [A1lin; A2lin];
blin = [b1lin; b2lin];
[pars, Jmin] = ga(@(x)Crit_J2(x,D,alpha,beta,r0,c0),2,Alin,blin);
disp('pars : ga');
pars
估计是：a = 3.74和b = 1.04（而不是起始固定值a = 14和b = 5）
我是否应该断定问题来自我的代码中使用的成本函数的定义？即：
% Function of cost
function cost = Crit_J2(p,D,alpha,beta,r0,c0)
% Compute the model corresponding to parameters p
[R,C] = size(D);
[Cols,Rows] = meshgrid(1:C,1:R);
% Model
Model = p(1)*(1+((Rows-r0).^2+(Cols-c0).^2)/alpha^2).^(-beta) + ...
p(2)*ones(R,C);
%Model = p(1)*D + p(2)*ones(R,C);
model = Model(:);
d = D(:);
%disp(model);
% Introduce H matrix
H = [ model, ones(length(model),1)];
% Compute the cost function : taking absolute value
J2 = abs((d-H*[p(1),p(2)]')'*(d-H*[p(1),p(2)]'));
cost = J2;
end
我正在使用以下成本函数：
J = abs((d-H*[p(1),p(2)]')'*(d-H*[p(1),p(2)]'));
withd估计开始时生成的初始数据向量，H矩阵使得d = H.thetawiththeta=[a,b]=[p(1),p(2)]'
[a,b]（恒星振幅+背景）这个错误估计（under estimation）的解释是什么？
现在您只是随机更改参数。如果您更改为 fminsearch 中实现的优化程序，此算法将为您找到最小值。
% initial parameters
a = 10;
b = 3;
r0 = 0;
c0 = 0;
alpha = 1;
beta = 1;
% Array of parameters
p0 = [a,b,r0,c0,alpha,beta]';
% Minimization of cost function for p
[pars, Jmin] = fminsearch(@(x)Crit_J(x,D), p0);
您必须将 Crit_J 函数更改为仅输出 'cost' 以最小化此功能。然后您可以稍后通过 运行 原始函数获得您的模型。
更新：
也许您正在寻找局部最小值，您可以尝试 fminunc instead of fminsearch, or switch to a generic algorithm，但这对于这个问题来说可能有点矫枉过正。