答張盛東──關於穩定算法的定義

網友張盛東留言:

如果老師有空希望老師可以討論一下數值線性代數中的穩定性 (stability) 和後向穩定性 (backward stability)。Trefethen 的 Numerical Linear Algebra 有兩句話特別令我費解[1]

A stable algorithm gives nearly the right answer to nearly the right question.

A backward stable algorithm gives exactly the right answer to nearly the right question.

請問老師,如何理解這兩句話?這個定義不是太直觀,希望老師給個具體例子說明一下這兩句話的深意。

 
答曰:

我們從線性代數所要解決的問題說起。譬如,解線性方程 A\mathbf{x}=\mathbf{b},計算 A 的特徵值 \lambda_1,\ldots,\lambda_n,其中 A 是給定的 n\times n 階矩陣,\mathbf{b}n 維向量。這些問題可以看成函數

y=f(x)

或記作映射 f:x\mapsto y,這裡 f 表示問題 (problem),x 表示問題的輸入 (或數據資料),y 表示問題的真解 (true solution)。例如,解線性方程是 f:(A,\mathbf{b})\mapsto\mathbf{x},計算特徵值是 f:A\mapsto(\lambda_1,\ldots,\lambda_n)。在數值線性代數中,一個特定問題經常存在許多不同的數值計算方法 (algorithm,簡稱算法)。理想上,我們希望算法能給出完全正確的數值解。但電腦的精確度是有限的,引入誤差在所難免,穩定性的意義即在告知我們何謂「正確的」答案,它是評估算法的主要指標,其中包括前向 (forward)、後向 (backward) 以及混和 (mixed) 穩定性。

 
針對問題 f,任何一個相應的算法可視為不同的函數 \hat{f},它將輸入 x 映至 \hat{y}。現實上,算法 \hat{f} 所得到的結果 \hat{y} 與真解 y 總是有一定的偏差,兩者的差異 \delta y=\hat{y}-y 稱為前向誤差 (forward error)。這個誤差通常來自算法所引入的捨入 (round-off) 誤差和截斷 (truncated) 誤差。下圖中,單箭頭實線表示問題 f,單箭頭虛線表示算法 \hat{f},雙箭頭實線表示前向誤差 \delta y

穩定算法1

一個準確的 (accurate) 算法所產生的 \hat{f}(x) 應該接近真實答案 f(x),也就是說,對於每一輸入 x,前向誤差 \delta y 應該非常小。所謂的「小」是一個相對的概念,不同的應用場合有不同的定義。在一般情況下,我們要求相對前向誤差 \Vert\delta y\Vert/\Vert y\Vert 必須與浮點算術的單位捨入或截斷誤差 (以 u 表示) 處於同一數量級,記為 \Vert\delta y\Vert/\Vert y\Vert=O(u),如此 \hat{f} 便是一個準確的算法。不過,許多問題並不存在滿足前述條件的準確算法。為了進一步界定算法的準確性,我們轉從問題的輸入下手。給定前向誤差 \delta y=\hat{y}-y,反過來問:「哪些輸入 \hat{x} 使得 f(\hat{x})=\hat{y}?」令 \delta x=\hat{x}-x。針對輸入 x,滿足上述關係的 \delta x 經常不止一個,於是我們定義後向誤差 (backward error) 為具有最小「長度」的 \delta x (即 \Vert\delta x\Vert 有最小值) 使得

f(x+\delta x)=\hat{f}(x)

下圖描繪前向、後向誤差,以及問題 f 與算法 \hat{f} 之間的關係。

穩定算法2

對於任一輸入 x,如果後向誤差 \delta x 滿足 \Vert \delta x\Vert/\Vert x\Vert=O(u),我們便說算法 \hat{f} 具有後向穩定性 (backward stability)。換句話說,後向穩定算法的結果 \hat{y} 是幾乎正確的輸入資料 \hat{x} (因為 \delta x=\hat{x}-x 很小) 所給出的正確解答 (因為 \hat{y}=f(\hat{x}))。這種基於後向誤差的作法稱為後向誤差分析。讀者或許有此疑問:為甚麼後向誤差不定義為 \delta x 使得 \hat{f}(x+\delta x)=f(x)?也就是說,輸入誤差 \delta x 的作用在於消彌算法 \hat{f} 產生的前向誤差。將 x=\hat{x}-\delta x 代入上式,可得 \hat{f}(\hat{x})=f(\hat{x}-\delta x),此即定義於輸入 \hat{x} 的後向誤差條件式,所以兩者確實同義,差別在於 \delta x 的符號相反而已。那麼,為甚麼要採用後向誤差分析而非前向誤差分析呢?主要原因有二[2]。第一,我們認為輸入擾動 (perturbation) 與捨入誤差是等價的概念。因為輸入資料不免包含來自多種原因的不確定擾動,如果後向誤差不大於這些不確定量,我們也就無從挑剔這個算法所得到的解。第二,今日已發展成熟的擾動理論可以從後向誤差估計前向誤差的上限 (稍後詳述)。

 
下面用一個例子來說明後向誤差分析。考慮四則運算問題 +, -, \times, /,令 \oplus, \ominus, \otimes, \oslash 表示電腦所執行的浮點計算函數 (即算法)。令函數 \mathrm{fl}(x) 表示 x 儲存於電腦的浮點數,下列性質成立[3]

\mathrm{fl}(x)=x(1+\epsilon),~~\vert\epsilon\vert\le u

由此可推論

x\circledast y=(x \ast y)(1+\epsilon),~~\vert\epsilon\vert\le u

上式中,符號 \ast 代表任何一種四則運算,\epsilon 是相對捨入誤差,因為 \epsilon=(\mathrm{fl}(x)-x)/x。下面我們分析浮點算術的減法運算 f(x_1,x_2)=x_1-x_2,對應的算法是

\hat{f}(x_1, x_2)=\mathrm{fl}(x_1)\ominus\mathrm{fl}(x_2)

利用浮點數性質 \mathrm{fl}(x_1)=x_1(1+\epsilon_1)\mathrm{fl}(x_2)=x_2(1+\epsilon_2),以及浮點計算函數性質 \mathrm{fl}(x_1)\ominus\mathrm{fl}(x_2)=\left[\mathrm{fl}(x_1)-\mathrm{fl}(x_2)\right](1+\epsilon_3),其中 \vert\epsilon_j\vert\le uj=1,2,3,可得

\begin{aligned}  \mathrm{fl}(x_1)\ominus\mathrm{fl}(x_2)&=\left[x_1(1+\epsilon_1)-x_2(1+\epsilon_2)\right](1+\epsilon_3)\\  &=x_1(1+\epsilon_1)(1+\epsilon_3)-x_2(1+\epsilon_2)(1+\epsilon_3)\\  &=x_1(1+\epsilon_4)-x_2(1+\epsilon_5),\end{aligned}

其中 \vert\epsilon_4\vert,\vert\epsilon_5\vert\le 2u+O(u^2)。令 \hat{x}_1=x_1(1+\epsilon_4)\hat{x}_2=x_2(1+\epsilon_5),可得

\hat{f}(x_1,x_2)=\hat{x}_1-\hat{x}_2=f(\hat{x}_1,\hat{x}_2)

其中 \vert\hat{x}_1-x_1\vert/\vert x_1\vert=O(u)\vert\hat{x}_2-x_2\vert/\vert x_2\vert=O(u)。所以,浮點減法運算具有後向穩定性。

 
然而,涉及簡單四則運算的問題未必都具有後向穩定性。令 \mathbf{x},\mathbf{y}\in\mathbb{R}^n,考慮秩-1 (rank-one) 矩陣計算問題 A=\mathbf{x}\mathbf{y}^T,即 a_{ij}=x_iy_ji,j=1,\ldots,n。對應的算法如下:

\hat{a}_{ij}=\mathrm{fl}(x_i)\otimes\mathrm{fl}(x_j),~~i,j=1,\ldots,n

因為誤差具有隨機性,\hat{A}=[\hat{a}_{ij}] 極可能不為秩-1矩陣,這時候不論如何選擇 \delta\mathbf{x}\delta\mathbf{y},下式都不會成立:

\hat{A}=(\mathbf{x}+\delta\mathbf{x})(\mathbf{y}+\delta\mathbf{y})^T

這說明秩-1矩陣算法不具後向穩定性。推廣至一般情形,對於問題 f 和算法 \hat{f},如果不存在極小的後向誤差 \delta x 使得 f(x+\delta x)=\hat{f}(x),退而求其次,我們可以放寬限制為

f(x+\delta x)=\hat{f}(x)+\Delta y,~~~\Vert\delta x\Vert/\Vert x\Vert=O(u),~~\Vert\Delta y\Vert/\Vert y\Vert=O(u)

對於每一輸入 x,若存在 \delta x\Delta y 滿足上述關係,則稱算法 \hat{f} 具有混和穩定性 (或簡稱穩定性)。換句話說,穩定算法的結果 \hat{y} 是幾乎正確的輸入資料 \hat{x} (因為 \delta x=\hat{x}-x 很小) 所給出的幾乎正確解答 (因為 \Delta y=f(\hat{x})-\hat{y} 也很小)。下面是混和誤差分析的示意圖。注意,\delta y\Delta y 是兩個不同的度量。

穩定算法3

 
接著我們討論如何利用後向誤差來估計前向誤差。考慮單變量問題 y=f(x),令近似解為 \tilde{y}=f(x+\delta x) (若 \Delta y=0,則 \tilde{y}=\hat{y},見上圖)。若 f 是二次連續可導函數,則

\displaystyle  \tilde{y}-y=f(x+\delta x)-f(x)=f'(x)\delta x+\frac{f''(x+\alpha\delta x)}{2}(\delta x)^2,~~\alpha\in[0,1]

上式等號兩邊同時除以 y,可得

\displaystyle  \frac{\tilde{y}-y}{y}=\left(\frac{xf'(x)}{f(x)}\right)\frac{\delta x}{x}+O\left((\delta x)^2\right)

如不計誤差平方項,相對前向誤差 \vert\tilde{y}-y\vert/\vert y\vert 和相對後向誤差 \vert\delta x\vert/\vert x\vert 的比值即是

\displaystyle  \kappa(x)=\left|\frac{xf'(x)}{f(x)}\right|

稱為 f 的相對條件數 (condition number)。將 f'\delta f/\delta x 取代,推廣至多變量 xf,我們定義相對條件數為

\displaystyle  \kappa(x)=\lim_{\eta\to 0}\max_{\Vert\delta x\Vert\le\eta}\left(\frac{\Vert\delta f\Vert}{\Vert f(x)\Vert}\bigg/\frac{\Vert\delta x\Vert}{\Vert x\Vert}\right)

其中 \Vert\cdot\Vert 是任何合法的向量或矩陣範數 (norm)。若 \delta x\delta f 都是無窮小,則 \kappa(x) 可簡化為 (見“條件數”)

\displaystyle  \kappa(x)=\max_{\delta x}\left(\frac{\Vert\delta f\Vert}{\Vert f(x)\Vert}\bigg/\frac{\Vert\delta x\Vert}{\Vert x\Vert}\right)

上式的意義是 \kappa(x) 等於 f 的相對誤差 (即前向誤差) 與 x 的相對誤差 (即後向誤差) 的最大比值,故

\displaystyle  \frac{\Vert\delta f\Vert}{\Vert f(x)\Vert}\le \kappa(x)\frac{\Vert\delta x\Vert}{\Vert x\Vert}

若相對條件數很大,輸入 x 的微小擾動便可能造成 f 的巨大改變,我們稱這樣的問題是病態的 (ill-conditioned);反之,則稱為良置的 (well-conditioned)。病態問題不能歸罪於算法引入的誤差,根本原因在於問題的「天生」特性 (見“病態系統”)。病態問題不存在可行的數學解決方案,數值計算方法僅適用於良置問題。

 
最後補充一個概念:前向穩定性。如果一個算法所產生的前向誤差大小近似其他的後向穩定算法的結果,我們說此算法具有前向穩定性。但前向穩定算法本身未必具有後向穩定性,也就是說,後向穩定性蘊含前向穩定性,但相反陳述並不成立。

 
引用來源:
[1] Lloyd N. Trefethen and David Bau III, Numerical Linear Algebra, 1997, pp 104.
[2] Nicholas J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd edition, 2002, pp 6-8.
[3] Gene H. Golub and Charles F. Van Loan, Matrix Computation, 2nd edition, 1989, pp 61-62.

This entry was posted in 特別主題, 答讀者問 and tagged , , , , . Bookmark the permalink.

12 則回應給 答張盛東──關於穩定算法的定義

  1. chenlogy 說:

    感覺是控制組研究所等級的問題了

    • ccjou 說:

      以上這些數值分析概念都是英國數學家 James H. Wilkinson 提出來的,目前是計算機上各種數值計算最常用的誤差分析手段 (維基百科如是說)。列入研究所等級問題一點也不為過。

      • chenlogy 說:

        我忘記要打上想請教的問題了….
        假如要設計一個zero forcing equalizer 發現原本的反矩陣不存在或難以實現,怎麼辦?
        如果採用MMSE的方法要怎麼樣不用DSP,只靠一般元件實現?

        • ccjou 說:

          上面部落格的副標題:I seek not to know only answers, but to understand the questions.

          我不很確定你的問題。原本的反矩陣是指
          http://zh.wikipedia.org/wiki/%E7%AD%89%E5%8C%96%E5%99%A8
          裡面的係數矩陣的反矩陣嗎?所謂「難以實現」所指為何?

          我可以寫篇RLS或與最佳估計相關的文章,或許這能回答你的問題。

          • chenlogy 說:

            是的,其實是從一本書上看到的
            [Stephen_H._Hall,_Howard_L._Heck]
            ADVANCED SIGNAL INTEGRITY FOR HIGH-SPEED DIGITAL DESIGNS

            12.3.2 Coefficient Selection
            At this point the question of how to determine the tap coefficients seems appropri-
            ate to consider. In this section we present a method that can be used to establish
            coefficient values based on specific performance criteria. The method, called a
            zero forcing solution (ZFS), sets the tap coefficients to force the equalizer output
            to match the values desired at all sample points. [Qureshi, 1985; Sklar, 2001].
            Development of the algorithm follows.
            想到的問題….此外如果雜訊會不會被反矩陣放大至無法判別訊號…
            還有一段…
            The adaptive ZFS approach has the advantage that it is
            very easy to implement but has the drawback that it does not comprehend ISI
            that occurs outside the length of the equalizer.
            還有,要怎麼樣用基本元件實現?
            我覺得我可能不適合唸研究所………

            • ccjou 說:

              你們學校有關於Wiener filters的課程嗎?Channel equalization(即ZFS)是它的一個典型應用。Wiener filters的基本構想是最小化估計誤差的MSE,估計誤差定義為目標響應與實際輸出的偏差。如果我沒有會錯意,你說的問題不會發生。

              除了DSP之外,我不知道牽涉矩陣運算的硬體實現是否還有其他實作方法。

  2. 張盛東 說:

    最近比較忙,今天才有機會細讀老師的大作。萬分感謝老師解惑!

  3. 張盛東 說:

    老師,請教一個問題。如何分析一個沒有輸入的算法的穩定性呢?如果連輸入都沒有,這個算法是否不可能是backward stable? 例如我在Numerical Linear Algebra 的第15章exercise看到兩道問題:
    15.1(e) Data: None. Solution: e, computed by summing \sum_{k=0}^{\infty}1/k! from left to right using \bigoplus and \bigotimes , stopping when a summand is reached of magnitude < \epsilon_{machine}.
    (f) Data: none. Solution: e, computed by the same algorithm as above except with the series summed from right to left。

    • ccjou 說:

      e 可以看成在 x=1 計算 e^x。當 x>0,以泰勒展開計算 e^x 是穩定的。

      我不太懂這個問題隱藏的含意,第一:(e) 和 (f) 的計算順序相反,這暗示兩個結果不同,為何運算順序造成差異?第二:\sum_{k=0}^n1/k! 的計算如何僅使用\bigoplus\bigotimes達成?

      • 張盛東 說:

        關於老師的兩個問題我也很困惑。第一,k!在電腦中是以整數表示而不是浮點數,不應該會引入誤差,除非k非常大而不得不採用浮點數。第二,每一個summand都需要一次除法獲得。我查了很多大學的網站都沒有關於這兩題的答案或者相關說明,不知道這兩道問題是否有誤。

        • ccjou 說:

          s_n=1/0!+1/1!+\cdots+1/n!,則
          \hat{s}_0=\mathrm{fl}(1)\oslash\mathrm{fl}(0!)=(1/1)(1+\epsilon_1)
          \hat{s}_1=(\hat{s}_0+(1/1)(1+\epsilon_2))(1+\epsilon_3)=((1/1)(1+\epsilon_1)+(1/1)(1+\epsilon_2))(1+\epsilon_3)=(1/1)(1+\epsilon_1)(1+\epsilon_3)+(1/1)(1+\epsilon_2)(1+\epsilon_3)
          其中 \epsilon_i=O(u)。以此類推,
          \hat{s}_n=(1/0!)(1+\epsilon)^{n+1}+(1/1!)(1+\epsilon)^{n}+(1/2!)(1+\epsilon)^{n-1}+\cdots+(1/n!)(1+\epsilon)^2
          停止時 \vert (1/n!)(1+\epsilon)\vert<\epsilon_{machine}。暫時先這樣。
          (如果採用相反方向計算,則有
          \tilde{s}_n=(1/n!)(1+\epsilon)^{n+1}+(1/(n-1)!)(1+\epsilon)^{n}+(1/(n-2)!)(1+\epsilon)^{n-1}+\cdots+(1/0!)(1+\epsilon)^2)
          再令 (1+\epsilon)^k\approx 1+k\epsilon,則後向誤差為 (1/0!)(n+1)\epsilon,(1/1!)(n)\epsilon,\ldots,(1/n!)(2)\epsilon。由此應可證明算法具有穩定性 (但並非後向穩定)。

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com Logo

你正使用 WordPress.com 帳號留言。 登出 / 變更 )

Twitter picture

你正使用 Twitter 帳號留言。 登出 / 變更 )

Facebook照片

你正使用 Facebook 帳號留言。 登出 / 變更 )

Google+ photo

你正使用 Google+ 帳號留言。 登出 / 變更 )

連結到 %s