在接下來的內容中,我們將提供一些關於數值最佳化演算法的筆記。由於存在大量方法 ,我們將限制自己討論所謂的梯度方法 。我們之所以將此類方法作為數值最佳化演算法的自然起點,主要有以下兩個原因。一方面,這些方法在該領域確實是主力軍,因此它們在實踐中的頻繁使用證明了它們在這裡被討論的合理性。另一方面,這種方法非常直觀,因為它在一定程度上是自然地從眾所周知的最優點性質 中推匯出來的。我們將重點關注此類方法的三個例子:牛頓法 、最速下降法 以及可變度量法 ,其中包括擬牛頓法 。
在我們開始之前,我們還是要強調一點,似乎並沒有一個“唯一”的演算法,特定演算法的效能總是取決於待解決的具體問題。因此,經驗和“試錯”在應用工作中非常重要。為了闡明這一點,我們將提供一些應用示例,其中可以透過圖形方式比較不同演算法的效能。此外,最後還提供了一個關於最大似然估計 的具體示例。對於統計學家和計量經濟學家 來說,最大似然估計可能是在實踐中必須依賴數值最佳化演算法的最重要示例。
任何數值最佳化演算法都需要解決找到函式的“可觀察”屬性的問題,以便計算機程式知道何時達到解。由於我們正在處理最佳化問題,因此兩個眾所周知的結果似乎是尋找這些屬性的合理起點。
如果 f 可微且 x ⋆ {\displaystyle x^{\star }} 是 (區域性) 最小值,那麼
( 1 a ) D f ( x ⋆ ) = 0 {\displaystyle (1a)\quad Df(x^{\star })=0}
即雅可比矩陣 D f ( x ) {\displaystyle Df(x)} 等於零
以及
如果 f 二階可微且 x ⋆ {\displaystyle x^{\star }} 是 (區域性) 最小值,那麼
( 1 b ) x T D 2 f ( x ⋆ ) x ≥ 0 {\displaystyle (1b)\quad x^{T}D^{2}f(x^{\star })x\geq 0}
即海森矩陣 D 2 f ( x ) {\displaystyle D^{2}f(x)} 是正半定 的。
在下文中,我們將始終用 x ⋆ {\displaystyle x^{\star }} 表示最小值。雖然這兩個條件似乎代表了幫助找到最優 x ⋆ {\displaystyle x^{\star }} 的陳述,但有一個小小的陷阱,它們給出了 x ⋆ {\displaystyle x^{\star }} 是函式 f {\displaystyle f} 的最優值的含義。但為了我們的目的,我們需要相反的含義,即最終我們想得到一個形式為:“如果某個條件 g ( f ( x ⋆ ) ) {\displaystyle g(f(x^{\star }))} 成立,則 x ⋆ {\displaystyle x^{\star }} 是最小值”的陳述。但上面兩個條件顯然不足以達到這一點(例如,考慮 f ( x ) = x 3 {\displaystyle f(x)=x^{3}} 的情況,其中 D f ( 0 ) = D 2 f ( 0 ) = 0 {\displaystyle Df(0)=D^{2}f(0)=0} ,但 x ⋆ ≠ 0 {\displaystyle x^{\star }\neq 0} )。因此,我們必須考慮 x ⋆ {\displaystyle x^{\star }} 的整個鄰域,如以下檢測最優值的充分條件所述。
如果 D f ( x ⋆ ) = 0 {\displaystyle Df(x^{\star })=0} 且 x T D 2 f ( z ) x ≥ 0 , ∀ x ∈ R n {\displaystyle x^{T}D^{2}f(z)x\geq 0,\forall x\in \mathbb {R} ^{n}} 且 z ∈ B ( x ⋆ , δ ) {\displaystyle z\in {\mathcal {B}}(x^{\star },\delta )} ,則: x ⋆ {\displaystyle x^{\star }} 是一個區域性最小值。
證明:對於 x ∈ B ( x ⋆ , δ ) {\displaystyle x\in {\mathcal {B}}(x^{\star },\delta )} 令 z = x ⋆ + t ( x − x ⋆ ) ∈ B {\displaystyle z=x^{\star }+t(x-x^{\star })\in {\mathcal {B}}} 。根據 泰勒近似 得到: f ( x ) − f ( x ⋆ ) = 0 + 1 2 ( x − x ⋆ ) T D 2 f ( z ) ( x − x ⋆ ) ≥ 0 {\displaystyle f(x)-f(x^{\star })=0+{\frac {1}{2}}(x-x^{\star })^{T}D^{2}f(z)(x-x^{\star })\geq 0} ,其中 B ( x ⋆ , δ ) {\displaystyle {\mathcal {B}}(x^{\star },\delta )} 表示以 x ⋆ {\displaystyle x^{\star }} 為中心的開球,即 B ( x ⋆ , δ ) = { x : | | x − x ⋆ | | ≤ δ } {\displaystyle {\mathcal {B}}(x^{\star },\delta )=\{x:||x-x^{\star }||\leq \delta \}} 對於 δ > 0 {\displaystyle \delta >0} 。
與上述兩個條件相反,此條件足以用於檢測最優值 - 考慮以下兩個簡單示例
f ( x ) = x 3 {\displaystyle f(x)=x^{3}} 其中 D f ( x ⋆ = 0 ) = 0 {\displaystyle Df(x^{\star }=0)=0} 但 x T D 2 f ( z ) x = 6 z x 2 ≱ 0 ( e . g . z = − δ 2 ) {\displaystyle x^{T}D^{2}f(z)x=6zx^{2}\not \geq 0\quad (e.g.\,\,z=-{\frac {\delta }{2}})}
以及
f ( x ) = x 4 {\displaystyle f(x)=x^{4}} 其中 D f ( x ⋆ = 0 ) = 0 {\displaystyle Df(x^{\star }=0)=0} 以及 x T D 2 f ( z ) x = 12 z 2 x 2 ≥ 0 ∀ z {\displaystyle x^{T}D^{2}f(z)x=12z^{2}x^{2}\geq 0\quad \forall z} 。
牢記這個小注意事項,我們現在可以轉向數值最佳化程式。
以下所有演算法都依賴於以下假設
(A1) 集合 N ( f , f ( x ( 0 ) ) = { x ∈ R n | f ( x ) ≤ f ( x ( 0 ) ) } {\displaystyle N(f,f(x^{(0)})=\{x\in \mathbb {R} ^{n}|f(x)\leq f(x^{(0)})\}} 是 緊緻 的。
其中 x ( 0 ) {\displaystyle x^{(0)}} 是演算法給定的某個初始值。這個假設的重要性體現在 魏爾斯特拉斯定理 中,該定理指出每個緊緻集都包含其 上確界 和 下確界 。因此,(A1) 保證了 N ( f , f ( x ( 0 ) ) {\displaystyle N(f,f(x^{(0)})} 中存在某個解。並且,在這個全域性最小值 x ⋆ {\displaystyle x^{\star }} 處,當然滿足 D ( f ( x ⋆ ) ) = 0 {\displaystyle D(f(x^{\star }))=0} 。因此 - 考慮到上面的討論 - 最佳化問題基本上歸結為求解方程組 D ( f ( x ⋆ ) ) = 0 {\displaystyle D(f(x^{\star }))=0} 。
當然,這種方法的問題是, D ( f ( x ⋆ ) ) = 0 {\displaystyle D(f(x^{\star }))=0} 也適用於 極大值和鞍點 。因此,好的演算法應該確保將極大值和鞍點排除為潛在解。透過要求 f ( x ( k + 1 ) ) < f ( x ( k ) ) {\displaystyle f(x^{(k+1)})<f(x^{(k)})} ,即我們限制在 序列 { x ( k ) } k {\displaystyle \{x^{(k)}\}_{k}} 上,使得函式值在每一步都減小,可以很容易地排除極大值。問題是,這是否總是可能的。幸運的是,答案是肯定的。為什麼是這樣的基本見解是:當構建對映 x ( k + 1 ) = φ ( x ( k ) ) {\displaystyle x^{(k+1)}=\varphi (x^{(k)})} (即我們如何從 x ( k ) {\displaystyle x^{(k)}} 到 x ( k + 1 ) {\displaystyle x^{(k+1)}} 的規則)時,我們有兩個自由度,即
方向 d ( k ) {\displaystyle d^{(k)}} 和
步長 σ ( k ) {\displaystyle \sigma ^{(k)}} 。
因此,我們可以選擇要移動的方向以到達 x ( k + 1 ) {\displaystyle x^{(k+1)}} 以及 此移動的距離。 因此,如果我們選擇 d ( k ) {\displaystyle d^{(k)}} 和 σ ( k ) {\displaystyle \sigma ^{(k)}} 以“正確的方式”,我們可以有效地確保函式值減小。 下面提供了這種推理的正式表示。
引理:如果 d ( k ) ∈ R n {\displaystyle d^{(k)}\in \mathbb {R} ^{n}} 且 D f ( x ) T d ( k ) < 0 {\displaystyle Df(x)^{T}d^{(k)}<0} ,那麼: ∃ σ ¯ > 0 {\displaystyle \exists {\bar {\sigma }}>0} 使得
f ( x + σ ( k ) d ( k ) ) < f ( x ) ∀ σ ∈ ( 0 , σ ¯ ) {\displaystyle f(x+\sigma ^{(k)}d^{(k)})<f(x)\quad \forall \sigma \in (0,{\bar {\sigma }})}
證明:由於 D f ( x ) T d ( k ) < 0 {\displaystyle Df(x)^{T}d^{(k)}<0} 且 D f ( x ) T d ( k ) = lim σ → 0 f ( x + σ ( k ) d ( k ) ) − f ( x ) σ ( k ) {\displaystyle Df(x)^{T}d^{(k)}=\lim _{\sigma \rightarrow 0}{\frac {f(x+\sigma ^{(k)}d^{(k)})-f(x)}{\sigma ^{(k)}}}} ,則對於足夠小的 σ ( k ) {\displaystyle \sigma ^{(k)}} , f ( x + σ ( k ) d ( k ) ) < f ( x ) {\displaystyle f(x+\sigma ^{(k)}d^{(k)})<f(x)} 成立。
滿足此條件的方向向量 d ( k ) {\displaystyle d^{(k)}} 被稱為下降方向 。 在實踐中,這個引理允許我們使用以下步驟來數值求解最佳化問題。
1. 定義序列 { x ( k ) } k {\displaystyle \{x^{(k)}\}_{k}} 透過遞推公式 x ( k + 1 ) = x ( k ) + σ ( k ) d ( k ) {\displaystyle x^{(k+1)}=x^{(k)}+\sigma ^{(k)}d^{(k)}} 。
2. 從點 x ( k ) {\displaystyle x^{(k)}} 的區域性資訊中選擇方向 d ( k ) {\displaystyle d^{(k)}} 。
3. 選擇一個步長 σ ( k ) {\displaystyle \sigma ^{(k)}} ,以確保演算法的收斂。
4. 如果 | f ( x ( k + 1 ) ) − f ( x ( k ) ) | < ϵ {\displaystyle |f(x^{(k+1)})-f(x^{(k)})|<\epsilon } ,其中 ϵ > 0 {\displaystyle \epsilon >0} 是某個選擇的最小值容差,則停止迭代。
此過程已經暗示了 d ( k ) {\displaystyle d^{(k)}} 和 σ ( k ) {\displaystyle \sigma ^{(k)}} 的選擇並非獨立的,而是相互依賴的。特別是需要注意的是,即使該方法是下降方法(即 d ( k ) {\displaystyle d^{(k)}} 和 σ ( k ) {\displaystyle \sigma ^{(k)}} 是根據引理 1 選擇的),也不能保證收斂到最小值。乍一看,這似乎有點令人費解。如果我們找到了一個序列 { x ( k ) } k {\displaystyle \{x^{(k)}\}_{k}} ,使得函式值在每一步都下降,人們可能會認為在某個階段,即在k 趨於無窮大的極限中,我們應該得到解。為什麼事實並非如此,可以從以下借鑑自 W. Alt (2002, p. 76) 的例子中看出。
考慮以下例子,它雖然明顯是下降的,但並不收斂。設判據函式由以下給出
f ( x ) = x 2 {\displaystyle f(x)=x^{2}} ,設初始值為 x ( 0 ) = 1 {\displaystyle x^{(0)}=1} ,考慮一個(常數)方向向量 d ( k ) = − 1 {\displaystyle d^{(k)}=-1} ,並選擇步長為 σ ( k ) = ( 1 2 ) k + 2 {\displaystyle \sigma ^{(k)}=({\frac {1}{2}})^{k+2}} 。因此,序列 { x ( k ) } k {\displaystyle \{x^{(k)}\}_{k}} 的遞迴定義如下
( 2 ) x ( k + 1 ) = x ( k ) + ( 1 2 ) k + 2 ( − 1 ) = x ( k − 1 ) − ( 1 2 ) k + 1 − ( 1 2 ) k + 2 = x ( 0 ) − ∑ j = 0 k ( 1 2 ) j + 2 . {\displaystyle (2)\quad x^{(k+1)}=x^{(k)}+({\frac {1}{2}})^{k+2}(-1)=x^{(k-1)}-\left({\frac {1}{2}}\right)^{k+1}-\left({\frac {1}{2}}\right)^{k+2}=x^{(0)}-\sum _{j=0}^{k}\left({\frac {1}{2}}\right)^{j+2}.}
注意, x ( k ) > 0 ∀ k {\displaystyle x^{(k)}>0\,\,\forall \,\,k} ,因此 f ( x ( k + 1 ) ) < f ( x ( k ) ) ∀ k {\displaystyle f(x^{(k+1)})<f(x^{(k)})\,\,\forall \,\,k} ,因此它顯然是一個下降法。然而,我們發現
( 3 ) lim k → ∞ x ( k ) = lim k → ∞ x ( 0 ) − ∑ j = 0 k − 1 ( 1 2 ) j + 2 = lim k → ∞ 1 − 1 4 ( 1 − ( 1 2 ) k 1 2 ) = lim k → ∞ 1 2 + ( 1 2 ) k + 1 = 1 2 ≠ 0 = x ⋆ . {\displaystyle (3)\quad \lim _{k\rightarrow \infty }x^{(k)}=\lim _{k\rightarrow \infty }x^{(0)}-\sum _{j=0}^{k-1}\left({\frac {1}{2}}\right)^{j+2}=\lim _{k\rightarrow \infty }1-{\frac {1}{4}}\left({\frac {1-\left({\frac {1}{2}}\right)^{k}}{\frac {1}{2}}}\right)=\lim _{k\rightarrow \infty }{\frac {1}{2}}+\left({\frac {1}{2}}\right)^{k+1}={\frac {1}{2}}\neq 0=x^{\star }.}
這種不收斂的原因必須歸因於步長 σ ( k ) {\displaystyle \sigma ^{(k)}} 下降得太快。對於較大的k ,步長 x ( k + 1 ) − x ( k ) {\displaystyle x^{(k+1)}-x^{(k)}} 變得如此之小,以至於收斂被排除在外。因此,我們必須將步長與下降方向 d ( k ) {\displaystyle d^{(k)}} 聯絡起來。
這種聯絡的明顯思路是要求實際下降與一階近似成正比,即選擇 σ ( k ) {\displaystyle \sigma ^{(k)}} 使得存在常數 c 1 > 0 {\displaystyle c_{1}>0} 使得
( 4 ) f ( x ( k ) + σ ( k ) d ( k ) ) − f ( x ( k ) ) ≤ c 1 σ ( k ) D ( f ( x ( k ) ) ) d ( k ) < 0. {\displaystyle (4)\quad f(x^{(k)}+\sigma ^{(k)}d^{(k)})-f(x^{(k)})\leq c_{1}\sigma ^{(k)}D(f(x^{(k)}))d^{(k)}<0.}
注意,我們仍然只關注下降方向,因此 D f ( x ( k ) ) T d ( k ) < 0 {\displaystyle Df(x^{(k)})^{T}d^{(k)}<0} ,如上文引理 1 中所述。因此, N ( f , f ( x ( k ) ) ) {\displaystyle N(f,f(x^{(k)}))} 的緊緻性意味著 LHS 的收斂 ,並且根據 (4)
( 5 ) lim k → ∞ σ ( k ) D ( f ( x ( k ) ) ) d ( k ) = 0. {\displaystyle (5)\quad \lim _{k\rightarrow \infty }\sigma ^{(k)}D(f(x^{(k)}))d^{(k)}=0.}
最後,我們希望選擇一個序列 { x ( k ) } k {\displaystyle \{x^{(k)}\}_{k}} 使得 lim k → ∞ D ( f ( x ( k ) ) ) = 0 {\displaystyle \lim _{k\rightarrow \infty }D(f(x^{(k)}))=0} ,因為這正是我們要解決的一階必要條件。在什麼條件下,(5) 實際上意味著 lim k → ∞ D ( f ( x ( k ) ) ) = 0 {\displaystyle \lim _{k\rightarrow \infty }D(f(x^{(k)}))=0} ? 首先,步長 σ ( k ) {\displaystyle \sigma ^{(k)}} 不能過快地趨於零。這正是我們在上面例子中遇到的情況。因此,似乎有道理從下界限制步長,要求
( 6 ) σ ( k ) ≥ − c 2 D f ( x ( k ) ) T d ( k ) | | d ( k ) | | 2 > 0 {\displaystyle (6)\quad \sigma ^{(k)}\geq -c_{2}{\frac {Df(x^{(k)})^{T}d^{(k)}}{||d^{(k)}||^{2}}}>0}
對於某個常數 c 2 > 0 {\displaystyle c_{2}>0} 。將 (6) 代入 (5) 最終得到
( 7 ) f ( x ( k ) + σ ( k ) d ( k ) ) − f ( x ( k ) ) ≤ − c ( D f ( x ( k ) ) T d ( k ) | | d ( k ) | | ) 2 , c = c 1 c 2 {\displaystyle (7)\quad f(x^{(k)}+\sigma ^{(k)}d^{(k)})-f(x^{(k)})\leq -c\left({\frac {Df(x^{(k)})^{T}d^{(k)}}{||d^{(k)}||}}\right)^{2},\quad c=c_{1}c_{2}}
其中,緊緻性 N ( f , f ( x ( k ) ) ) {\displaystyle N(f,f(x^{(k)}))} 確保了左側的收斂 ,因此
( 8 ) lim k → ∞ − c ( D f ( x ( k ) ) T d ( k ) | | d ( k ) | | ) 2 = lim k → ∞ D f ( x ( k ) ) T d ( k ) | | d ( k ) | | = 0 {\displaystyle (8)\quad \lim _{k\rightarrow \infty }-c({\frac {Df(x^{(k)})^{T}d^{(k)}}{||d^{(k)}||}})^{2}=\lim _{k\rightarrow \infty }{\frac {Df(x^{(k)})^{T}d^{(k)}}{||d^{(k)}||}}=0}
滿足 (4) 和 (6) 的步長被稱為有效步長 ,並用 σ E ( k ) {\displaystyle \sigma _{E}^{(k)}} 表示。條件 (6) 的重要性在以下示例 1 的延續中得到說明。
需要注意的是,正是 (6) 的失敗導致示例 1 未能收斂。將示例中的步長代入 (6) 中得到
( 6.1 ) σ ( k ) = ( 1 2 ) ( k + 2 ) ≥ − c 2 2 x ( k ) ( − 1 ) 1 = c 2 ⋅ 2 ( 1 2 + ( 1 2 ) k + 1 ) ⇔ 1 4 ( 1 + 2 ( k ) ) ≥ c 2 > 0 {\displaystyle (6.1)\quad \sigma ^{(k)}=\left({\frac {1}{2}}\right)^{(k+2)}\geq -c_{2}{\frac {2x^{(k)}(-1)}{1}}=c_{2}\cdot 2({\frac {1}{2}}+\left({\frac {1}{2}})^{k+1}\right)\Leftrightarrow {\frac {1}{4(1+2^{(k)})}}\geq c_{2}>0}
因此,對於所有k ,沒有 常數 c 2 > 0 {\displaystyle c_{2}>0} 滿足 (6) 中要求的這個不等式。因此,步長沒有從下方界定,並且下降速度過快。為了真正認識到 (6) 的重要性,讓我們稍微改變一下這個例子,假設 σ ( k ) = ( 1 2 ) k + 1 {\displaystyle \sigma ^{(k)}=\left({\frac {1}{2}}\right)^{k+1}} 。然後我們發現
( 6.2 ) lim k → ∞ x ( k + 1 ) = lim k → ∞ x ( 0 ) − 1 2 ∑ i ( 1 2 ) i = lim k → ∞ ( 1 2 ) k + 1 = 0 = x ⋆ , {\displaystyle (6.2)\quad \lim _{k\rightarrow \infty }x^{(k+1)}=\lim _{k\rightarrow \infty }x^{(0)}-{\frac {1}{2}}\sum _{i}\left({\frac {1}{2}}\right)^{i}=\lim _{k\rightarrow \infty }\left({\frac {1}{2}}\right)^{k+1}=0=x^{\star },}
也就是說,收斂 確實發生了。此外,認識到這個例子實際上確實滿足條件 (6),因為
( 6.3 ) σ ( k ) = ( 1 2 ) ( k + 1 ) ≥ − c 2 2 x ( k ) ( − 1 ) 1 = c 2 ⋅ 2 ( 1 2 ) k ⇔ 1 4 ≥ c 2 > 0. {\displaystyle (6.3)\quad \sigma ^{(k)}=\left({\frac {1}{2}}\right)^{(k+1)}\geq -c_{2}{\frac {2x^{(k)}(-1)}{1}}=c_{2}\cdot 2\left({\frac {1}{2}}\right)^{k}\Leftrightarrow {\frac {1}{4}}\geq c_{2}>0.}
我們已經論證過,選擇 σ ( k ) {\displaystyle \sigma ^{(k)}} 和 d ( k ) {\displaystyle d^{(k)}} 是交織在一起的。因此,選擇“正確”的 d ( k ) {\displaystyle d^{(k)}} 總是取決於相應的步長 σ ( k ) {\displaystyle \sigma ^{(k)}} 。那麼在這個語境下,“正確”是什麼意思呢?上面我們在公式 (8) 中表明,選擇一個有效的步長 意味著
( 8 ′ ) lim k → ∞ − c ( D f ( x ( k ) ) T d ( k ) | | d ( k ) | | ) 2 = lim k → ∞ D f ( x ( k ) ) T d ( k ) | | d ( k ) | | = 0. {\displaystyle (8')\quad \lim _{k\rightarrow \infty }-c({\frac {Df(x^{(k)})^{T}d^{(k)}}{||d^{(k)}||}})^{2}=\lim _{k\rightarrow \infty }{\frac {Df(x^{(k)})^{T}d^{(k)}}{||d^{(k)}||}}=0.}
因此,“正確”的方向向量將保證 (8') 意味著
( 9 ) lim k → ∞ D f ( x ( k ) ) = 0 {\displaystyle (9)\quad \lim _{k\rightarrow \infty }Df(x^{(k)})=0}
因為 (9) 是選擇序列 { x ( k ) } k {\displaystyle \{x^{(k)}\}_{k}} 收斂的條件。所以讓我們探究一下哪些方向可以得到 (9)。假設步長 σ k {\displaystyle \sigma _{k}} 是有效 的,並且定義
( 10 ) β ( k ) = D f ( x ( k ) ) T d ( k ) | | D f ( x ( k ) ) | | | | d ( k ) | | ⇔ β ( k ) | | D f ( x ( k ) ) | | = D f ( x ( k ) ) T d ( k ) | | d ( k ) | | {\displaystyle (10)\quad \beta ^{(k)}={\frac {Df(x^{(k)})^{T}d^{(k)}}{||Df(x^{(k)})||||d^{(k)}||}}\quad \Leftrightarrow \quad \beta ^{(k)}||Df(x^{(k)})||={\frac {Df(x^{(k)})^{T}d^{(k)}}{||d^{(k)}||}}}
根據 (8') 和 (10) 我們知道
( 11 ) lim k → ∞ β ( k ) | | D f ( x ( k ) ) | | = 0. {\displaystyle (11)\quad \lim _{k\rightarrow \infty }\beta ^{(k)}||Df(x^{(k)})||=0.}
所以如果我們從下方對 β ( k ) {\displaystyle \beta ^{(k)}} 進行邊界設定(即 β ( k ) ≤ − δ < 0 {\displaystyle \beta ^{(k)}\leq -\delta <0} ),(11) 意味著
( 12 ) lim k → ∞ β ( k ) | | D f ( x ( k ) ) | | = lim k → ∞ | | D f ( x ( k ) ) | | = lim k → ∞ D f ( x ( k ) ) = 0 , {\displaystyle (12)\quad \lim _{k\rightarrow \infty }\beta ^{(k)}||Df(x^{(k)})||=\lim _{k\rightarrow \infty }||Df(x^{(k)})||=\lim _{k\rightarrow \infty }Df(x^{(k)})=0,}
其中,(12) 只給出了序列 { x ( k ) } k {\displaystyle \{x^{(k)}\}_{k}} 收斂到解 x ⋆ {\displaystyle x^{\star }} 的條件。由於 (10) 透過 β ( k ) {\displaystyle \beta ^{(k)}} 隱式地定義了方向向量 d ( k ) {\displaystyle d^{(k)}} ,因此對 β ( k ) {\displaystyle \beta ^{(k)}} 的要求直接轉化為對 d ( k ) {\displaystyle d^{(k)}} 的要求。
當考慮對 β ( k ) {\displaystyle \beta ^{(k)}} 的條件時,很明顯 “梯度方法” 這個術語的來源。對於由
β k = D ( f ( x ) ) d ( k ) | | D f ( x ( k ) ) | | | | d ( k ) | | = cos ( D f ( x ( k ) ) , d ( k ) ) {\displaystyle \beta _{k}={\frac {D(f(x))d^{(k)}}{||Df(x^{(k)})||||d^{(k)}||}}=\cos(Df(x^{(k)}),d^{(k)})}
給出的 β ( k ) {\displaystyle \beta ^{(k)}} ,我們有以下結果:
假設 σ ( k ) {\displaystyle \sigma ^{(k)}} 已經被有效地選擇,並且 d ( k ) {\displaystyle d^{(k)}} 滿足
( 13 ) cos ( D f ( x ( k ) ) , d ( k ) ) = β k ≤ − δ < 0 {\displaystyle (13)\quad \cos(Df(x^{(k)}),d^{(k)})=\beta _{k}\leq -\delta <0}
我們有
( 14 ) lim k → ∞ D f ( x ( k ) ) → 0 {\displaystyle (14)\quad \lim _{k\rightarrow \infty }Df(x^{(k)})\rightarrow 0}
因此:如果 x ( k ) {\displaystyle x^{(k)}} 處的負梯度與方向 d ( k ) {\displaystyle d^{(k)}} 之間的夾角始終小於直角,那麼就會發生收斂。依賴於滿足 (13) 的 d ( k ) {\displaystyle d^{(k)}} 的方法被稱為梯度方法。
換句話說:只要不沿梯度的正交 方向移動,並且有效地選擇步長,梯度方法 就能保證收斂到解 x ⋆ {\displaystyle x^{\star }} 。
現在讓我們探索這類別中三種具體的演算法,它們在各自的 d ( k ) {\displaystyle d^{(k)}} 的選擇方面有所不同。
牛頓法 是該領域中最流行的方法。它是一種求解所有型別方程根 的著名方法,因此可以輕鬆地應用於最佳化問題。牛頓法的主要思想是將方程組線性化,得到
( 15 ) g ( x ) = g ( x ^ ) + D g ( x ^ ) T ( x − x ^ ) = 0. {\displaystyle (15)\quad g(x)=g({\hat {x}})+Dg({\hat {x}})^{T}(x-{\hat {x}})=0.}
(15) 可以很容易地解出 x ,因為解僅由(假設 D g ( x ^ ) T {\displaystyle Dg({\hat {x}})^{T}} 為非奇異 )給出
( 16 ) x = x ^ − [ D g ( x ^ ) T ] − 1 g ( x ^ ) . {\displaystyle (16)\quad x={\hat {x}}-[Dg({\hat {x}})^{T}]^{-1}g({\hat {x}}).}
為了我們的目的,我們只選擇 g ( x ) {\displaystyle g(x)} 為梯度 D f ( x ) {\displaystyle Df(x)} ,得到
( 17 ) d N ( k ) = x ( k + 1 ) − x ( k ) = − [ D 2 f ( x ( k ) ) ] − 1 D f ( x ( k ) ) {\displaystyle (17)\quad d_{N}^{(k)}=x^{(k+1)}-x^{(k)}=-[D^{2}f(x^{(k)})]^{-1}Df(x^{(k)})}
其中 d N ( k ) {\displaystyle d_{N}^{(k)}} 被稱為 *牛頓方向*。
分析 (17) 可以揭示牛頓法的主要性質
如果 D 2 f ( x ( k ) ) {\displaystyle D^{2}f(x^{(k)})} 是一個 正定矩陣 ,那麼 d N k {\displaystyle d_{N}^{k}} 是 *引理 1* 中提到的一個下降方向。
牛頓法使用一階 *和* 二階導數的區域性資訊來計算 d N k {\displaystyle d_{N}^{k}} 。
( 18 ) x ( k + 1 ) = x ( k ) + d N ( k ) {\displaystyle (18)\quad x^{(k+1)}=x^{(k)}+d_{N}^{(k)}}
牛頓法使用固定步長 σ ( k ) = 1 {\displaystyle \sigma ^{(k)}=1} 。因此牛頓法 *不一定是* *引理 1* 中提到的下降法。原因是固定步長 σ ( k ) = 1 {\displaystyle \sigma ^{(k)}=1} 可能大於 *引理 1* 中給出的臨界步長 σ k ¯ {\displaystyle {\bar {\sigma _{k}}}} 。下面我們將給出 *牛頓法* 不下降的 *羅森布羅克函式* 作為示例。
該方法可能非常耗時,因為每一步 *k* 都要計算 [ D 2 f ( x ( k ) ) ] − 1 {\displaystyle [D^{2}f(x^{(k)})]^{-1}} 可能很麻煩。在實際應用中,我們可以考慮使用近似方法。例如,我們可以每隔 *s* 步更新一次 Hessian 矩陣,或者我們可以依靠 *區域性近似*。這被稱為 *擬牛頓法*,將在討論 *變數度量法* 的部分進行介紹。
為了確保該方法是下降法,我們可以使用一個有效的步長 σ E ( k ) {\displaystyle \sigma _{E}^{(k)}} 並設定
( 19 ) x ( k + 1 ) = x ( k ) − σ E ( k ) d N ( k ) = x ( k ) − σ E ( k ) [ D 2 f ( x k ) ] − 1 D f ( x ( k ) ) {\displaystyle (19)\quad x^{(k+1)}=x^{(k)}-\sigma _{E}^{(k)}d_{N}^{(k)}=x^{(k)}-\sigma _{E}^{(k)}[D^{2}f(x^{k})]^{-1}Df(x^{(k)})}
另一個常用的方法是最速下降法 。這種方法的思路是選擇方向 d ( k ) {\displaystyle d^{(k)}} ,使得函式值f 的下降最大。雖然乍一看這個過程似乎很有道理,但它實際上利用的資訊比牛頓法 少,因為它忽略了Hessian關於函式曲率的資訊。尤其是在下面的應用中,我們將看到一些關於這個問題的例子。
最速下降法 的方向向量由下式給出:
( 20 ) d S D ( k ) = a r g m a x d : | | d | | = r { − D f ( x ( k ) ) T d } = a r g m i n d : | | d | | = r { D f ( x ( k ) ) T d } = − r D f ( x ) | | D f ( x ) | | {\displaystyle (20)\quad d_{SD}^{(k)}=argmax_{d:||d||=r}\{-Df(x^{(k)})^{T}d\}=argmin_{d:||d||=r}\{Df(x^{(k)})^{T}d\}=-r{\frac {Df(x)}{||Df(x)||}}}
證明:根據柯西-施瓦茨不等式 ,可得
( 21 ) D f ( x ) T d | | D f ( x ) | | | | d | | ≥ − 1 ⇔ D f ( x ) T d ≥ − r | | D f ( x ) | | . {\displaystyle (21)\quad {\frac {Df(x)^{T}d}{||Df(x)||||d||}}\geq -1\quad \Leftrightarrow \quad Df(x)^{T}d\geq -r||Df(x)||.}
顯然 (21) 等式對 d ( k ) = d S D ( k ) {\displaystyle d^{(k)}=d_{SD}^{(k)}} 在 (20) 中給出。
特別要注意,對於 r = | | D f ( x ) | | {\displaystyle r=||Df(x)||} ,我們有 d S D ( k ) = − D f ( x ( k ) ) {\displaystyle d_{SD}^{(k)}=-Df(x^{(k)})} ,即我們只是在負梯度的方向上“行走”。與 *牛頓法* 不同的是,*最速下降法* 不使用固定步長,而是選擇一個有效的步長 σ E ( k ) {\displaystyle \sigma _{E}^{(k)}} 。因此,*最速下降法* 定義了序列 { x ( k ) } k {\displaystyle \{x^{(k)}\}_{k}} 為
( 22 ) x ( k + 1 ) = x ( k ) + σ E ( k ) d S D ( k ) , {\displaystyle (22)\quad x^{(k+1)}=x^{(k)}+\sigma _{E}^{(k)}d_{SD}^{(k)},}
其中 σ E ( k ) {\displaystyle \sigma _{E}^{(k)}} 是一個有效的步長, d S D ( k ) {\displaystyle d_{SD}^{(k)}} 是在 (20) 中給出的最速下降方向。
對於 d S D ( k ) = − r D f ( x ) | | D f ( x ) | | {\displaystyle d_{SD}^{(k)}=-r{\frac {Df(x)}{||Df(x)||}}} ,最速下降法定義了一個下降方向,從 *引理 1* 的意義上說,因為
D f ( x ) T d S D ( k ) = D f ( x ) T ( − r D f ( x ) | | D f ( x ) | | ) = − r | | D f ( x ) | | D f ( x ) T D f ( x ) < 0. {\displaystyle Df(x)^{T}d_{SD}^{(k)}=Df(x)^{T}\left(-r{\frac {Df(x)}{||Df(x)||}}\right)=-{\frac {r}{||Df(x)||}}Df(x)^{T}Df(x)<0.}
*最速下降法* 只是區域性上可行的,因為它忽略了二階資訊。
特別是當目標函式很平坦時(即解 x ⋆ {\displaystyle x^{\star }} 位於“山谷”中),最速下降法定義的序列會劇烈波動(參見下面的應用,特別是 Rosenbrock 函式的例子)。
由於它不需要 Hessian 矩陣,因此 *最速下降法* 的計算和實現很簡單快捷。
比 *牛頓法* 和 *最速下降法* 更通用的方法是 *可變度量法* 類。這類方法依賴於更新公式
( 23 ) x k + 1 = x k − σ E ( k ) [ A k ] − 1 D f ( x k ) . {\displaystyle (23)\quad x^{k+1}=x^{k}-\sigma _{E}^{(k)}[A^{k}]^{-1}Df(x^{k}).}
如果 A k {\displaystyle A^{k}} 是一個 對稱 且 正定 矩陣,(23) 定義了一個下降方法,其中 [ A k ] − 1 {\displaystyle [A^{k}]^{-1}} 是正定的當且僅當 A k {\displaystyle A^{k}} 也是正定的。要理解這一點,只需考慮 譜分解
( 24 ) A k = Γ Λ Γ T {\displaystyle (24)\quad A^{k}=\Gamma \Lambda \Gamma ^{T}}
其中 Γ {\displaystyle \Gamma } 和 Λ {\displaystyle \Lambda } 分別是包含 特徵向量 和 特徵值 的矩陣。如果 A k {\displaystyle A^{k}} 是正定的,所有特徵值 λ i {\displaystyle \lambda _{i}} 嚴格為正。因此它們的逆 λ i − 1 {\displaystyle \lambda _{i}^{-1}} 也是正的,因此 [ A k ] − 1 = Γ Λ − 1 Γ T {\displaystyle [A^{k}]^{-1}=\Gamma \Lambda ^{-1}\Gamma ^{T}} 明顯是正定的。但是,將 d ( k ) = [ A k ] − 1 D f ( x k ) {\displaystyle d^{(k)}=[A^{k}]^{-1}Df(x^{k})} 代入,得到
( 25 ) D f ( x k ) T d ( k ) = − D f ( x k ) T [ A k ] − 1 D f ( x k ) ≡ − v T [ A k ] − 1 v ≤ 0 , {\displaystyle (25)\quad Df(x^{k})^{T}d^{(k)}=-Df(x^{k})^{T}[A^{k}]^{-1}Df(x^{k})\equiv -v^{T}[A^{k}]^{-1}v\leq 0,}
也就是說,該方法確實是下降的。到目前為止,我們還沒有指定矩陣 A k {\displaystyle A^{k}} ,但很容易看出,對於兩種特定的選擇,可變度量法 恰好與最速下降法 和牛頓法 重合。
對於 A k = I {\displaystyle A^{k}={\mathcal {I}}} (其中 I {\displaystyle {\mathcal {I}}} 是 單位矩陣 ),可以得出
( 22 ′ ) x k + 1 = x k − σ E ( k ) D f ( x k ) {\displaystyle (22')\quad x^{k+1}=x^{k}-\sigma _{E}^{(k)}Df(x^{k})}
這僅僅是*最速下降法*。
對於 A k = D 2 f ( x k ) {\displaystyle A^{k}=D^{2}f(x^{k})} 可以得出
( 19 ′ ) x k + 1 = x k − σ E ( k ) [ D 2 f ( x k ) ] − 1 D f ( x k ) {\displaystyle (19')\quad x^{k+1}=x^{k}-\sigma _{E}^{(k)}[D^{2}f(x^{k})]^{-1}Df(x^{k})}
這僅僅是使用步長 σ E ( k ) {\displaystyle \sigma _{E}^{(k)}} 的*牛頓法*。
對於*變數度量法*,另一個自然候選者是*擬牛頓法*。與標準*牛頓法*相比,它使用有效的步長,使其成為下降法;與*最速下降法*相比,它不會完全忽略函式曲率的區域性資訊。因此,*擬牛頓法*由矩陣 A k {\displaystyle A^{k}} 上的兩個要求定義:
A k {\displaystyle A^{k}} 應該近似於Hessian D 2 f ( x k ) {\displaystyle D^{2}f(x^{k})} 以利用關於曲率的資訊,並且
更新 A k → A k + 1 {\displaystyle A^{k}\rightarrow A^{k+1}} 應該很簡單,以便演算法仍然相對快速(即使在高維度)。
為了確保第一個要求, A k + 1 {\displaystyle A^{k+1}} 應該滿足所謂的*擬牛頓方程*
( 26 ) A k + 1 ( x ( k + 1 ) − x ( k ) ) = D f ( x ( k + 1 ) ) − D f ( x ( k ) ) {\displaystyle (26)\quad A^{k+1}(x^{(k+1)}-x^{(k)})=Df(x^{(k+1)})-Df(x^{(k)})}
因為所有滿足 (26) 的 A k {\displaystyle A^{k}} 都反映了Hessian 的資訊。為了說明這一點,考慮定義為
( 27 ) g ( x ) = f ( x k + 1 ) + D f ( x k + 1 ) T ( x − x k + 1 ) + 1 2 ( x − x k + 1 ) T A k + 1 ( x − x k + 1 ) . {\displaystyle (27)\quad g(x)=f(x^{k+1})+Df(x^{k+1})^{T}(x-x^{k+1})+{\frac {1}{2}}(x-x^{k+1})^{T}A^{k+1}(x-x^{k+1}).}
顯然, g ( x k + 1 ) = f ( x k + 1 ) {\displaystyle g(x^{k+1})=f(x^{k+1})} 以及 D g ( x k + 1 ) = D f ( x k + 1 ) {\displaystyle Dg(x^{k+1})=Df(x^{k+1})} 。因此 g ( x ) {\displaystyle g(x)} 和 f ( x ) {\displaystyle f(x)} 在 x ( k + 1 ) {\displaystyle x^{(k+1)}} 的鄰域內相當相似。為了確保 g ( x ) {\displaystyle g(x)} 在 x ( k ) {\displaystyle x^{(k)}} 處也是一個好的近似值,我們希望選擇 A k + 1 {\displaystyle A^{k+1}} 使得 x ( k ) {\displaystyle x^{(k)}} 處的梯度相同。有了
( 28 ) D g ( x k ) = D f ( x k + 1 ) − A k + 1 ( x k + 1 − x k ) {\displaystyle (28)\quad Dg(x^{k})=Df(x^{k+1})-A^{k+1}(x^{k+1}-x^{k})}
很明顯,如果 A k + 1 {\displaystyle A^{k+1}} 滿足 (26) 中給出的擬牛頓方程 ,則 D g ( x k ) = D f ( x k ) {\displaystyle Dg(x^{k})=Df(x^{k})} 。但是,隨之而來的是
( 29 ) A k + 1 ( x k + 1 − x k ) = D f ( x k + 1 ) − D g ( x k ) = D f ( x k + 1 ) − D f ( x k ) = D 2 f ( λ x ( k ) + ( 1 − λ ) x ( k + 1 ) ) ( x k + 1 − x k ) . {\displaystyle (29)\quad A^{k+1}(x^{k+1}-x^{k})=Df(x^{k+1})-Dg(x^{k})=Df(x^{k+1})-Df(x^{k})=D^{2}f(\lambda x^{(k)}+(1-\lambda )x^{(k+1)})(x^{k+1}-x^{k}).}
因此,只要 x ( k + 1 ) {\displaystyle x^{(k+1)}} 和 x ( k ) {\displaystyle x^{(k)}} 彼此相差不遠, A k + 1 {\displaystyle A^{k+1}} 滿足(26)是 D 2 f ( x ( k ) ) {\displaystyle D^{2}f(x^{(k)})} 的一個良好近似。
現在讓我們來談談第二個要求,即 A k {\displaystyle A^{k}} 的更新應該很容易。一種特定的演算法是所謂的 BFGS 演算法 。該演算法的主要優點是它只使用已經計算出的元素 { x ( k ) } k {\displaystyle \{x^{(k)}\}_{k}} 和 { D f ( x ( k ) ) } k {\displaystyle \{Df(x^{(k)})\}_{k}} 來構建更新 A ( k + 1 ) {\displaystyle A^{(k+1)}} 。因此,不需要計算任何新的實體,只需跟蹤 _x_ 序列和梯度序列即可。作為 BFGS 演算法的起點,可以提供任何正定矩陣(例如,在 x ( 0 ) {\displaystyle x^{(0)}} 處的 Hessian 或單位矩陣)。然後,_BFGS 更新公式_ 為
( 30 ) A k = A k − 1 − ( A k − 1 ) T γ k − 1 T γ k − 1 A k − 1 γ k − 1 T A k − 1 γ k − 1 + Δ k − 1 Δ k − 1 T Δ k − 1 T γ k − 1 {\displaystyle (30)\quad A^{k}=A^{k-1}-{\frac {(A^{k-1})^{T}\gamma _{k-1}^{T}\gamma _{k-1}A^{k-1}}{\gamma _{k-1}^{T}A^{k-1}\gamma _{k-1}}}+{\frac {\Delta _{k-1}\Delta _{k-1}^{T}}{\Delta _{k-1}^{T}\gamma _{k-1}}}}
其中 Δ k − 1 = D f ( x ( k ) ) − D f ( x ( k − 1 ) ) {\displaystyle \Delta _{k-1}=Df(x^{(k)})-Df(x^{(k-1)})} 以及 γ k − 1 = x ( k ) − x ( k − 1 ) {\displaystyle \gamma _{k-1}=x^{(k)}-x^{(k-1)}} 。此外,(30) 確保所有 A k {\displaystyle A^{k}} 均為正定矩陣,正如 *Variable Metric Methods* 所要求的那樣,它們是下降的。
它利用了關於 f ( x ) {\displaystyle f(x)} 曲率的二階資訊,因為矩陣 A k {\displaystyle A^{k}} 與 Hessian 矩陣 D 2 f ( x ) {\displaystyle D^{2}f(x)} 相關。
儘管如此,它確保了簡單快速的更新(例如,透過 BFGS 演算法),因此它比標準牛頓法更快。
它是一種下降方法,因為 A k {\displaystyle A^{k}} 是正定的。
它相對容易實現,因為 BFGS 演算法在大多數數值或統計軟體包中都可用。
為了比較這些方法並說明演算法之間的差異,我們現在將評估 *最速下降法*、標準 *牛頓法* 和使用高效步長的 *擬牛頓法* 的效能。我們在這個領域使用兩個經典函式,即 Himmelblau 函式和 Rosenbrock 函式。
Himmelblau 函式定義為
( 31 ) f ( x , y ) = ( x 2 + y − 11 ) 2 + ( x + y 2 − 7 ) 2 {\displaystyle (31)\quad f(x,y)=(x^{2}+y-11)^{2}+(x+y^{2}-7)^{2}}
這個四階多項式有四個極小值、四個鞍點和一個極大值,因此演算法有足夠多的可能性失敗。在下圖中,我們展示了不同起始值的函式的 *等高線圖* 和 3D 圖。
在圖 1 中,我們展示了函式以及三種方法在起始值為 ( 2 , − 4 ) {\displaystyle (2,-4)} 時的路徑。顯然,三種方法沒有找到相同的極小值。原因當然是因為 *最速下降法* 的方向向量不同——透過忽略關於曲率的資訊,它選擇的方向與兩個 *牛頓法* 完全不同(特別是在圖 1 的右圖中)。
圖 1:兩種牛頓法收斂到同一個極小值,而最速下降法收斂到另一個極小值。
現在考慮起始值為 ( 4.5 , − 0.5 ) {\displaystyle (4.5,-0.5)} ,如圖 2 所示。當然,最重要的是現在 *所有* 方法都找到了不同的解。*最速下降法* 找到了與兩個 *牛頓法* 不同的解,這並不奇怪。但是,兩個 *牛頓法* 收斂到不同的解,這表明了步長 σ {\displaystyle \sigma } 的重要性。*擬牛頓法* 在第一次迭代中選擇了一個高效的步長,因此兩種方法在第一次迭代之後的所有迭代中都具有不同的步長 *和* 方向向量。正如圖片所示:結果可能相當顯著。
圖 2:即使所有方法都找到了不同的解。
Rosenbrock 函式定義為
( 32 ) f ( x , y ) = 100 ( y − x 2 ) 2 + ( 1 − x ) 2 {\displaystyle (32)\quad f(x,y)=100(y-x^{2})^{2}+(1-x)^{2}\ }
儘管該函式只有一個最小值,但對於最佳化問題而言,它仍然是一個有趣的函式。原因是該U形函式的谷底非常平坦(參見圖3和圖4的右側面板)。對於計量經濟學家而言,這個函式可能非常有趣,因為在最大似然估計的情況下,平坦的準則函式非常常見。因此,下面圖3和圖4中顯示的結果對於具有該問題的函式來說似乎是比較通用的。
我在使用這個函式和演算法時的經驗是,圖3(給定一個初始值 ( 2 , − 5 ) {\displaystyle (2,-5)} )似乎非常具有代表性。與上面的 Himmelblau 函式相比,所有演算法都找到了相同的解,並且由於只有一個最小值,所以這是可以預期的。更重要的是,不同方法選擇的路徑反映了各自方法的不同性質。可以看到,*最速下降法*波動得比較劇烈。這是由於該方法不使用關於曲率的資訊,而是來回跳動於與谷底相鄰的“山丘”之間。兩種牛頓方法選擇更直接的路徑,因為它們使用了二階資訊。兩種牛頓方法之間的主要區別當然是步長。圖3顯示,*擬牛頓法*在谷底穿梭時使用非常小的步長。相比之下,*牛頓法*的步長是固定的,因此它直接跳向解的方向。儘管人們可能會認為這是*擬牛頓法*的一個缺點,但請注意,通常這些較小的步長具有更高的穩定性的優點,即演算法不太可能跳到另一個解。這一點可以在圖4中看到。
圖3:所有方法都找到了相同的解,但最速下降法波動很大。
圖4考慮的是一個初始值 ( − 2 , − 2 ) {\displaystyle (-2,-2)} ,展示了*牛頓法*使用固定步長的主要問題——該方法可能會“過沖”,即它沒有下降。在第一步中,*牛頓法*(圖中顯示為紫色線)跳出谷底,然後在下一步中反彈回來。在這種情況下,仍然可以收斂到最小值,因為每側的梯度都指向中心唯一的谷底,但可以很容易地想象出不符合這種情況的函式。跳躍的原因是二階導數非常小,以至於步長 [ D f ( x ( k ) ) ] − 1 D f ( x ( k ) ) ) {\displaystyle [Df(x^{(k)})]^{-1}Df(x^{(k)}))} 由於 Hessian 的逆而變得非常大。根據我的經驗,我建議使用有效的步長來更好地控制各自方法選擇的路徑。
圖2:由於固定步長導致的牛頓法過沖。
對於計量經濟學家和統計學家來說,*最大似然估計器*可能是數值最佳化演算法最重要的應用。因此,我們將簡要介紹估計過程如何融入上面開發的框架。
像往常一樣,設
( 33 ) f ( Y | X ; θ ) {\displaystyle (33)\quad f(Y|X;\theta )}
為給定 *X* 的 *Y* 的*條件密度*,其中引數為 θ {\displaystyle \theta } ,並且
( 34 ) l ( θ ; Y | X ) {\displaystyle (34)\quad l(\theta ;Y|X)}
為引數 θ {\displaystyle \theta } 的*條件似然函式*。
如果我們假設資料是*獨立同分布(iid)*,則樣本對數似然函式如下所示:
( 35 ) L ( θ ; Y 1 , … , Y N ) = ∑ i N L ( θ ; Y i ) = ∑ i N log ( l ( θ ; Y i ) ) . {\displaystyle (35)\quad {\mathcal {L}}(\theta ;Y_{1},\dots ,Y_{N})=\sum _{i}^{N}{\mathcal {L}}(\theta ;Y_{i})=\sum _{i}^{N}\log(l(\theta ;Y_{i})).}
因此,最大似然估計歸結為關於引數 θ {\displaystyle \theta } 最大化 (35)。為了簡單起見,如果我們決定使用 *牛頓法* 來解決這個問題,序列 { θ ( k ) } k {\displaystyle \{\theta ^{(k)}\}_{k}} 由遞迴定義:
( 36 ) D θ L ( θ ( k + 1 ) ) = D θ L ( θ ( k ) ) + D θ θ L ( θ ( k ) ) ( θ ( k + 1 ) − θ ( k ) ) = 0 ⇔ θ ( k + 1 ) = θ ( k ) − [ D θ θ L ( θ ( k ) ) ] − 1 D θ L ( θ ( k ) ) {\displaystyle (36)\quad D_{\theta }{\mathcal {L}}(\theta ^{(k+1)})=D_{\theta }{\mathcal {L}}(\theta ^{(k)})+D_{\theta \theta }{\mathcal {L}}(\theta ^{(k)})(\theta ^{(k+1)}-\theta ^{(k)})=0\Leftrightarrow \theta ^{(k+1)}=\theta ^{(k)}-[D_{\theta \theta }{\mathcal {L}}(\theta ^{(k)})]^{-1}D_{\theta }{\mathcal {L}}(\theta ^{(k)})}
其中 D θ L {\displaystyle D_{\theta }{\mathcal {L}}} 和 D θ θ L {\displaystyle D_{\theta \theta }{\mathcal {L}}} 分別表示關於引數向量 θ {\displaystyle \theta } 的一階和二階導數,而 [ D θ θ L ( θ ( k ) ) ] − 1 D θ L ( θ ( k ) ) {\displaystyle [D_{\theta \theta }{\mathcal {L}}(\theta ^{(k)})]^{-1}D_{\theta }{\mathcal {L}}(\theta ^{(k)})} 定義了 (17) 中給出的 *牛頓方向*。由於最大似然估計始終假設條件密度(即誤差項的分佈)在引數 θ {\displaystyle \theta } 範圍內是已知的,因此上述方法可以很容易地應用。
假設一個簡單的線性模型
( 37 a ) Y i = β 1 + β x X i + U i {\displaystyle (37a)\quad Y_{i}=\beta _{1}+\beta _{x}X_{i}+U_{i}}
其中 θ = ( β 1 , β 2 ) ′ {\displaystyle \theta =(\beta _{1},\beta _{2})'} . 然後條件分佈 Y 由 U 的分佈決定,即
( 37 b ) p ( Y i − β 1 − β x X i ) ≡ p ∣ X i ( Y i ) = p ( U i ) , {\displaystyle (37b)\quad p(Y_{i}-\beta _{1}-\beta _{x}X_{i})\equiv p_{\mid X_{i}}(Y_{i})=p(U_{i}),}
其中 p {\displaystyle p} 表示密度函式 。一般來說,最大化 (35) 沒有閉式解(至少如果 U 碰巧不是正態分佈 ),因此必須使用數值方法。因此假設 U 遵循學生 t 分佈 ,具有 m 自由度 ,因此 (35) 由下式給出
( 38 ) L ( θ ; Y ∣ X ) = ∑ log ( Γ ( m + 1 2 ) π m Γ ( m 2 ) ( 1 + ( y i − x i T β ) 2 m ) − m + 1 2 ) {\displaystyle (38)\quad {\mathcal {L}}(\theta ;Y_{\mid X})=\sum \log \left({\frac {\Gamma ({\frac {m+1}{2}})}{{\sqrt {\pi m}}\Gamma ({\frac {m}{2}})}}(1+{\frac {(y_{i}-x_{i}^{T}\beta )^{2}}{m}})^{-{\frac {m+1}{2}}}\right)}
這裡我們只使用了 t 分佈的密度函式的定義。 (38) 可以簡化為
( 39 ) L ( θ ; Y ∣ X ) = N [ log ( Γ ( m + 1 2 ) ) − log ( π m Γ ( m 2 ) ) ] − m + 1 2 ∑ log ( 1 + ( y i − x i T β ) 2 m ) {\displaystyle (39)\quad {\mathcal {L}}(\theta ;Y_{\mid X})=N\left[\log \left(\Gamma \left({\frac {m+1}{2}}\right)\right)-\log \left({\sqrt {\pi m}}\Gamma \left({\frac {m}{2}}\right)\right)\right]-{\frac {m+1}{2}}\sum \log \left(1+{\frac {(y_{i}-x_{i}^{T}\beta )^{2}}{m}}\right)}
因此(如果我們假設自由度 m 是已知的)
( 40 ) argmax { L ( θ ; Y ∣ X ) } = argmax { − m + 1 2 ∑ log ( 1 + ( y i − x i T β ) 2 m ) } = argmin { ∑ log ( 1 + ( y i − x i T β ) 2 m ) } . {\displaystyle (40)\quad {\textit {argmax}}\left\{{\mathcal {L}}(\theta ;Y_{\mid X})\right\}={\textit {argmax}}\left\{-{\frac {m+1}{2}}\sum \log \left(1+{\frac {(y_{i}-x_{i}^{T}\beta )^{2}}{m}}\right)\right\}={\textit {argmin}}\left\{\sum \log \left(1+{\frac {(y_{i}-x_{i}^{T}\beta )^{2}}{m}}\right)\right\}.}
使用以下準則函式
( 41 ) f ( β 1 , β 2 ) = ∑ log ( 1 + ( y i − β 1 − β 2 x i ) 2 m ) {\displaystyle (41)\quad f(\beta _{1},\beta _{2})=\sum \log \left(1+{\frac {(y_{i}-\beta _{1}-\beta _{2}x_{i})^{2}}{m}}\right)}
上述方法可以很容易地應用於計算最大似然估計 ( β ^ 1 , M L , β ^ 2 , M L ) {\displaystyle ({\hat {\beta }}_{1,ML},{\hat {\beta }}_{2,ML})} ,最大化 (41)。
Alt, W. (2002): "Nichtlineare Optimierung", Vieweg: Braunschweig/Wiesbaden
Härdle, W. 和 Simar, L. (2003): "應用多元統計分析", Springer: Berlin Heidelberg
Königsberger, K. (2004): "分析 I", Springer: Berlin Heidelberg
Ruud, P. (2000): "古典計量經濟學理論", 牛津大學出版社: 紐約