在此記錄使用matlab作梯度下降法(GD)求函數極值的一個例子：
問題設定：
1. 我們有一個$n$個數據點，每個數據點是一個$d$維的向量，向量組成一個data矩陣$\mathbf{X}\in \mathbb{R}^{n\times d}$，這是我們的輸入特征矩陣。
2. 我們有一個響應的響應向量$\mathbf{y}\in \mathbb{R}^n$。
3. 我們將使用線性模型來fit上述數據。因此我們將優化問題形式化成如下形式：$$\arg\min_{\mathbf{w}}f(\mathbf{w})=\frac{1}{n}\|\mathbf{y}-\mathbf{\overline{X}}\mathbf{w}\|_2^2$$
其中$\mathbf{\overline{X}}=(\mathbf{1,X})\in \mathbb{R}^{n\times (d+1)}$ and $\mathbf{w}=(w_0,w_1,...,w_d)^\top\in \mathbb{R}^{d+1}$
顯然這是一個回歸問題，我們的目標從通俗意義上講就是尋找合適的權重向量$\mathbf{w}$，使得線性模型能夠擬合的更好。
預處理：
1. 按列對數據矩陣進行最大最小歸一化，該操作能夠加快梯度下降的速度，同時保證了輸入的數值都在0和1之間。$\mathbf{x}_i$為$\mathbf{X}$的第i列。 $$z_{ij}\leftarrow \frac{x_{ij}-\min(\mathbf{x}_i)}{\max(\mathbf{x}_i)-\min(\mathbf{x}_i)}$$
這樣我們的優化問題得到了轉化:$$\arg\min_{\mathbf{u}}g(\mathbf{w})=\frac{1}{n}\|\mathbf{y}-\mathbf{\overline{Z}}\mathbf{u}\|_2^2$$
2. 考慮對目標函數的Lipschitz constants進行估計。因為我們使用線性回歸模型，Lipschitz constants可以方便求得，這樣便於我們在梯度下降法是選擇合適的步長。假如非線性模型，可能要用其他方法進行估計（可選）。
問題解決：
使用梯度下降法進行問題解決，算法如下：
我們可以看到，這里涉及到求目標函數$f$對$\mathbf{x}_k$的梯度。顯然在這里，因為是線性模型，梯度的求解十分的簡單：$$\nabla f(\mathbf{x}_k)=-\frac{2}{n}\mathbf{\overline{X}}^\top(\mathbf{y}-\mathbf{\overline{X}}\mathbf{u}_k)$$
進行思考，還有沒有其他辦法可以把這個梯度給弄出來？假如使用Tensorflow，Pytorch這樣可以自動保存計算圖的東東，那么梯度是可以由機器自動求出來的。當然在這里我是用matlab實現，暫時沒有發現這樣的利器，所以我認為假如在這里想求出梯度，那么我們必須要把梯度的閉式解搞出來，不然沒法繼續進行。
下面是一段matlab的代碼：
function [g_result,u_result] = GD(N_Z,y,alpha,u0)
%GD 梯度下降法
% Detailed explanation goes here
[n,~] = size(N_Z);
u = u0;
k = 0;
t = y-N_Z*u;
disp("g(u):");
while(合理的終止條件)
k = k + 1;
u = u - alpha * (-2/n)*N_Z'*t;
t = y-N_Z*u;
if(mod(k,10)==0)
disp(t'*t/n);
end
end
g_result = (y-N_Z * u)' * (y-N_Z * u)/n;
u_result = u;
end
當然假如初始化的時候$u_0$選擇不當，而且因為沒有正則項，以上的算法將會有很大的問題：梯度消失，導致優化到最后的時候非常慢。我花了好多個小時才將loss講到0.19左右，而閉式解法能夠使得loss為0.06幾，運行時間也不會難以忍受。
問題推廣：
在這里，我們的問題是線性模型，回歸問題。能否有更廣的應用？思考后認為，只要需要優化的目標是標量，且該目標函數對輸入向量的梯度容易求得即可。只是因為該算法簡單朴素，可能在實際應用的時候會碰見惱人的梯度消失問題。