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

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.

$y=f(x)$

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

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

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

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

\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}

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

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

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

$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)$

$\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]$

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

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

$\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)$

$\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)$

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

[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.

### 12 Responses to 答張盛東──關於穩定算法的定義

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.

我不很確定你的問題。原本的反矩陣是指
裡面的係數矩陣的反矩陣嗎？所謂「難以實現」所指為何？

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

• chenlogy 說道：

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

12.3.2 Coefﬁcient Selection
At this point the question of how to determine the tap coefﬁcients seems appropri-
ate to consider. In this section we present a method that can be used to establish
coefﬁcient values based on speciﬁc performance criteria. The method, called a
zero forcing solution (ZFS), sets the tap coefﬁcients to force the equalizer output
to match the values desired at all sample points. [Qureshi, 1985; Sklar, 2001].
Development of the algorithm follows.
想到的問題….此外如果雜訊會不會被反矩陣放大至無法判別訊號…
還有一段…
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$。由此應可證明算法具有穩定性 (但並非後向穩定)。

• 張盛東 說道：

謝謝老師解答。