最優(yōu)化模型與算法-基于Python實(shí)現(xiàn) 教案 ch06 凸優(yōu)化算法

凸優(yōu)化算法1．編程實(shí)現(xiàn)梯度投影算法以求解具有盒子約束的最小二乘問題： (6.7.1)其中，為已知參數(shù)?？勺孕性O(shè)定數(shù)據(jù)及參數(shù)的值，步長(zhǎng)分別使用固定步長(zhǎng)規(guī)則和沿可行方向的Armijo規(guī)則。繪制迭代點(diǎn)處的目標(biāo)函數(shù)值、梯度值的模長(zhǎng)隨迭代次數(shù)的變化曲線。略。2．編程實(shí)現(xiàn)坐標(biāo)下降法以求解嶺回歸模型： (6.7.2)其中，為觀測(cè)樣本，，正則化項(xiàng)，為給定的正則化參數(shù)。坐標(biāo)的選取分別使用隨機(jī)化坐標(biāo)下降法(將指標(biāo)隨機(jī)排序后進(jìn)行遍歷)和順序遍歷(順序遍歷指標(biāo))，比較算法性能是否有差異。繪制迭代點(diǎn)處的目標(biāo)函數(shù)值隨迭代次數(shù)的變化曲線，輸出最優(yōu)解、最優(yōu)值和迭代次數(shù)，并與嶺回歸模型的解析解進(jìn)行比較。略。3．編程實(shí)現(xiàn)坐標(biāo)下降法以求解Lasso模型： (6.7.3)其中，正則化項(xiàng)。其他的設(shè)定與上題中的坐標(biāo)下降法相同。注意Lasso模型不再具有解析解。略。ISTA算法的改進(jìn)。本章介紹了ISTA算法求解Lasso模型，實(shí)際上，該算法可以進(jìn)一步改進(jìn)，A.Beck和M.Teboulle等人在2009年提出了一種加速方法FISTA[7]，顯著加快了ISTA算法的收斂速度。編程實(shí)現(xiàn)FISTA算法，并與ISTA的性能進(jìn)行比較。importnumpyasnp#定義L1范數(shù)的軟閾值函數(shù)defsoft_thresholding(x,threshold):returnnp.sign(x)*np.maximum(np.abs(x)-threshold,0)#FISTA算法deffista_regression(X,y,lambda_,iterations=100,learning_rate=0.01):n_samples,n_features=X.shapew=np.zeros((n_features,1))z=np.zeros((n_features,1))t=1foriinrange(iterations):#計(jì)算梯度gradient=2*np.dot(X.T,(np.dot(X,w)-y))/n_samples#更新權(quán)重w_new=soft_thresholding(z-learning_rate*gradient,lambda_*learning_rate)#使用加速技巧更新z和tt_new=(1+np.sqrt(1+4*t**2))/2z=w_new+(t-1)/t_new*(w_new-w)w=w_newt=t_newreturnw#ISTA算法defista_regression(X,y,lambda_,iterations=100,learning_rate=0.01):n_samples,n_features=X.shapew=np.zeros((n_features,1))foriinrange(iterations):#計(jì)算梯度gradient=2*np.dot(X.T,(np.dot(X,w)-y))/n_samples#更新權(quán)重w=soft_thresholding(w-learning_rate*gradient,lambda_*learning_rate)returnw#示例數(shù)據(jù)X=np.array([[1,2],[2,3],[3,4]])y=np.array([[2],[3],[4]])#設(shè)置超參數(shù)lambda_=0.1iterations=100learning_rate=0.01#使用FISTA算法進(jìn)行回歸fista_weights=fista_regression(X,y,lambda_,iterations,learning_rate)#使用ISTA算法進(jìn)行回歸ista_weights=ista_regression(X,y,lambda_,iterations,learning_rate)print("FISTA算法得到的權(quán)重：")print(fista_weights)print("\nISTA算法得到的權(quán)重：")print(ista_weights)5．編程實(shí)現(xiàn)交替方向乘子(ADMM)算法以求解廣義的Lasso模型： (6.7.4)其中，。為如下形式的差分矩陣： (6.7.5)略。6．編程實(shí)現(xiàn)SVRG和SAGA算法，分別求解嶺回歸模型： (6.7.6)其中，為觀測(cè)樣本，，正則化項(xiàng)，為給定的正則化參數(shù)。importnumpyasnp#定義sigmoid函數(shù)defsigmoid(x):return1/(1+np.exp(-x))#嶺回歸目標(biāo)函數(shù)的梯度計(jì)算defridge_regression_gradient(X,y,w,A):n_samples,n_features=X.shapegradient=2*np.dot(X.T,(np.dot(X,w)-y))/n_samples+2*A*wreturngradient#SVRG算法defsvrg_ridge_regression(X,y,A,learning_rate=0.01,iterations=100,m=100):n_samples,n_features=X.shapew=np.zeros((n_features,1))#計(jì)算所有樣本的梯度full_gradient=ridge_regression_gradient(X,y,w,A)foriinrange(iterations):#從所有樣本中隨機(jī)選擇m個(gè)樣本indices=np.random.choice(n_samples,m,replace=False)X_batch=X[indices]y_batch=y[indices]#計(jì)算當(dāng)前批次樣本的梯度batch_gradient=ridge_regression_gradient(X_batch,y_batch,w,A)#使用全局梯度更新w_oldw_old=w.copy()w=w-learning_rate*full_gradient#使用batch梯度更新wforjinrange(m):index=indices[j]x_j=X[index,:].reshape(-1,1)y_j=y[index]w=w-learning_rate*(ridge_regression_gradient(x_j,y_j,w_old,A)-batch_gradient+full_gradient)returnw#SAGA算法defsaga_ridge_regression(X,y,A,learning_rate=0.01,iterations=100,m=100):n_samples,n_features=X.shapew=np.zeros((n_features,1))#計(jì)算每個(gè)樣本的梯度，并初始化平均梯度為全局梯度gradients=[]foriinrange(n_samples):x_i=X[i,:].reshape(-1,1)y_i=y[i]gradient_i=ridge_regression_gradient(x_i,y_i,w,A)gradients.append(gradient_i)avg_gradient=np.mean(np.array(gradients),axis=0)foriinrange(iterations):#從所有樣本中隨機(jī)選擇m個(gè)樣本indices=np.random.choice(n_samples,m,replace=False)X_batch=X[indices]y_batch=y[indices]#計(jì)算當(dāng)前批次樣本的變化量delta_gradients=[]forjinrange(m):index=indices[j]x_j=X[index,:].reshape(-1,1)y_j=y[index]gradient_j=ridge_regression_gradient(x_j,y_j,w,A)delta_gradient_j=gradient_j-gradients[index]gradients[index]=gradient_jdelta_gradients.append(delta_gradient_j)#更新平均梯度delta_avg_gradient=np.mean(np.array(delta_gradients),axis=0)avg_gradient=avg_gradient+delta_avg_gradient#更新權(quán)重w=w-learning_rate*avg_gradientreturnw#示例數(shù)據(jù)X=np.random.randn(100,10)w_true=np.random.randn(10,1)y=sigmoid(np.dot(X,w_true)+0.1*np.random.randn(100,1))#設(shè)置超參數(shù)A=0.1learning_rate=0.01iterations=1000m=50#使用SVRG算法