二次函式的形式如下
設問題的維度n=158,取初始點為全為0的n維向量。
由於問題的形式特殊,所以步長α採用精確線搜尋的顯示表示。
對於函式的引數g,b,採用隨機生成。
function [fun,grad,hess,b]=f(x)
n=158;
a=unidrnd(10,n,1);
g=a*a'+unidrnd(2)*eye(n);
b=0.5*g*ones(n,1);
grad=g*x-b; %由於函式的係數是隨機生成的,所以在下面程式中的迭代都是採用該方法
hess=g;
fun=0.5*x'*g*x-b'*x;
end
function [x_star]=steepest(x0,eps)
k=0;
xk=x0; %選擇初始點
maxit=100000; %最大迭代次數
[~,gradk,g,b]=f(xk); %設定相關引數
while 1
res=norm(gradk); %gradk的二範數
k=k+1; %第k次迭代
if res<=eps||k>maxit %演算法停止條件
if k>maxit
fprintf('');
endbreak;
enddk=-gradk; %確定下降方向
alphak=(-dk'*gradk)/(dk'*g*dk); %計算步長
xk=xk+alphak*dk; %迭代
gradk=g*xk-b;
fprintf('at the %d-th iteration,the residual is %f\n',k,res);
endx_star=xk;
function [x_star]=newton(x0,eps)
k=0;
xk=x0; %選擇初始點
maxit=1000000; %最大迭代次數
[~,gradk,g,b]=f(xk); %設定相關引數
while 1
res=norm(gradk); %gradk的二範數
k=k+1; %第k次迭代
if res<=eps||k>maxit %演算法停止條件
break;
enddk=-inv(g)*gradk; %牛頓法確定下降方向
alphak=(-dk'*gradk)/(dk'*g*dk); %計算步長
xk=xk+alphak*dk; %迭代
gradk=g*xk-b;
fprintf('at the %d-th iteration,the residual is %f\n',k,res);
endx_star=xk;
function [x_star]=bfgs(x0,eps)
k=0;
xk=x0; %選擇初始點
maxit=100000; %最大迭代次數
hk=eye(length(x0)); %hk矩陣初始化
i=eye(length(x0));
[~,gradk,g,b]=f(xk); %設定相關引數
while 1
res=norm(gradk); %gradk的二範數
k=k+1; %第k次迭代
if res<=eps||k>maxit %演算法停止條件
break;
enddk=-hk*gradk; %確定下降方向
alphak=-dk'*gradk/(dk'*g*dk);%計算步長
xk=xk+alphak*dk; %迭代核心
%更新hk矩陣,利用sherman-morrison-woodbury公式
grad_temp=gradk; %計算yk
gradk=g*xk-b;
yk=gradk-grad_temp;
deltak=g\yk; %計算deltak
temp1=i-(deltak*yk')/(deltak'*yk); %更新hk矩陣
temp2=(deltak*deltak')/(deltak'*yk);
hk=temp1*hk*temp1+temp2;
fprintf('at the %d-th iteration,the residual is %f\n',k,res);
endx_star=xk;
function [x_star]=fr(x0,eps)
k=0;
xk=x0; %選擇初始點
maxit=10000; %最大迭代次數
[~,gradk,g,b]=f(xk); %設定相關引數
while 1
res=norm(gradk); %gradk的二範數
k=k+1; %第k次迭代
if res<=eps||k>maxit %演算法停止條件
break;
enddk=-gradk; %確定下降方向
alphak=-gradk'*dk/(dk'*g*dk); %計算步長
xk=xk+alphak*dk; %迭代核心
%計算beta
grad_temp=gradk;
gradk=g*xk-b;
betak=gradk'*gradk/(grad_temp'*grad_temp);
dk=-gradk+betak*dk;
fprintf('at the %d-th iteration,the residual is %f\n',k,res);
endx_star=xk;
執行結果如下:
>> steepest(zeros(158,1),1e-4)
at the 1-th iteration,the residual is 34754.134481
at the 2-th iteration,the residual is 3.002778
at the 3-th iteration,the residual is 1.633110
at the 4-th iteration,the residual is 0.000141
>> newton(zeros(158,1),1e-4)
at the 1-th iteration,the residual is 37774.400535
>> bfgs(zeros(158,1),1e-4)
at the 1-th iteration,the residual is 34511.098765
at the 2-th iteration,the residual is 5.758404
>> fr(zeros(158,1),1e-4)
at the 1-th iteration,the residual is 27961.323766
at the 2-th iteration,the residual is 6.037534
at the 3-th iteration,the residual is 3.306320
at the 4-th iteration,the residual is 0.000714
at the 5-th iteration,the residual is 0.000391
最後極小點為乙個全是0.5的n維向量。
其實根據函式的構造來看,儘管係數是隨機的,但是對於極小點有gx=b,解就很明顯了。
牛頓法一次迭代就能得到結果,得益於使用精確線搜尋,二次函式具有顯示表達。
最速下降法相比其他方法,顯然不是「最速」的。。。
擬牛頓法與牛頓法的區別在於迭代矩陣hk的不同。
求一元二次函式的根
總時間限制 1000ms記憶體限制 65536kb 描述利用公式x1 b sqrt b b 4 a c 2 a x2 b sqrt b b 4 a c 2 a 求一元二次方程ax2 bx c 0的根,其中a不等於0。輸入輸入一行,包含三個浮點數a,b,c 它們之間以乙個空格分開 分別表示方程ax2 ...
matlab 三元二次函式求最值
fmincon函式 今天晚上幫姐姐求乙個方程的最值 果斷用matlab啊 剛開始想得挺簡單的,就是for迴圈 後來一想計算量太大 可是算出來的是最小值 然後又找求最大值的函式 可是找了半天沒找到 最後我機智的把函式前面加了乙個負號 哈哈貼一下函式吧,挺簡單的 比如說y 99.27932 20.217...
凸優化問題,凸二次規劃問題QP,凸函式
約束優化問題 凸函式凸優化問題 凸二次規劃問題 約束優化問題 min w f w min w f w s.t.gi w 0 i 1,k 1 s.t.gi w 0 i 1,k 1 hj w 0 j 1,l 2 hj w 0 j 1,l 2 注 這是乙個最小化問題.不等式約束嚴格執行的含義是 小於等於號...