3
$\begingroup$

Suppose that $A\mathbf{x}=\mathbf{b}$ is a diagonally dominant linear system. We can use iterative methods that produce a sequence of approximations $\mathbf{x}^1,\mathbf{x}^2,\dots$ that converge to $\mathbf{x}$.

Doing this for a small system it seems intuitive to stop when the components $x^i_j$ and $x^{i+1}_j$ differ by a tolerance $\varepsilon$ giving solutions correct to $\varepsilon$.

While playing around with larger systems it was clear that this could stop well before convergence and is a very poor stopping rule (and doesn't give accurate answers as claimed).

Researching a little further I see some slightly more sophisticated stopping rules (in the below the norms might be max norms or 2-norms):

  • When $\|\mathbf{x}^{i+1}-\mathbf{x}^i\|< \varepsilon.$
  • When $\|\mathbf{x}^{i+1}-\mathbf{x}^i\|< \varepsilon\|\mathbf{b}\|.$
  • When $\displaystyle \max_j\left|\frac{x_j^{i+1}-x_j^i}{x_j^{i+1}}\right|<\varepsilon$.
  • When the residual $\|\mathbf{r}\|=\|\mathbf{b}-A\mathbf{x}^i\|<\varepsilon$.
  • When the residual $\|\mathbf{r}\|=\|\mathbf{b}-A\mathbf{x}^i\|<\varepsilon\|\mathbf{b}\|$.

I am interested in hearing the pros and cons of each. I am hoping to use one that is relatively easy to implement on VBA with perhaps 20 unknowns.

$\endgroup$
  • 2
    $\begingroup$ I have a perhaps stupid question: why do you want to solve a system with 20 unknowns iteratively? $\endgroup$ – Algebraic Pavel Sep 11 '17 at 22:22
  • 1
    $\begingroup$ There are no stupid questions. Approximating the equilibrium temperature distribution on a plate with fixed boundary conditions (i.e. numerical methods of solving Laplace's Equation). $\endgroup$ – JP McCarthy Sep 12 '17 at 7:50
  • $\begingroup$ Could I ask another question? Why implement in VBA? Is this integrating with some kind of embedded hardware with limited options for programming interface? :) $\endgroup$ – Erick Wong Sep 13 '17 at 20:03
  • $\begingroup$ @ErickWong this is service teaching for an engineering department. They want us to use VBA. $\endgroup$ – JP McCarthy Sep 14 '17 at 11:12
  • $\begingroup$ Small remark. Many algorithms (like CG) already compute the residual. So, residual based convergence criteria contain almost no additional computations. $\endgroup$ – P. Siehr Sep 15 '17 at 9:01
3
+100
$\begingroup$

Let's look at forward error bounds.

From normwise analysis we have $$\frac{\|x-x^i\|}{\|x\|} \leq \frac{2\kappa(A)\epsilon_n}{1-\kappa(A)\epsilon_n},\ \ \epsilon_n = \frac{\|b - Ax^i\|}{\|A\|\|x^i\| + \|b\|} \tag{1}$$ where $\|\cdot\|$ is any vector norm and the corresponding subordinate matrix norm, and $\kappa(A) = \|A\|\|A^{-1}\|$ is the matrix condition number with respect to the norm $\|\cdot\|$. Number $\epsilon_n$ is the normwise backward error obtained at $i$-th iteration. This is the only number, that can be controlled in this inequality.

From componentwise analysis we have: $$\frac{\|x-x^i\|_\infty}{\|x\|_\infty} \leq \frac{2\eta(A)\epsilon_c}{1-\eta(A)\epsilon_c},\ \ \epsilon_c = \max_j\frac{|b - Ax^i|_j}{(|A||x^i| + |b|)_j} \tag{2}$$ where $\eta(A) = \||A^{-1}||A|\|_\infty$ is the Skeel condition number and $\epsilon_c$ is the componentwise backward error obtained at $i$-th iteration.

From (1) and (2) we can see, that none of provided stopping rules are particullary good. The best one is the last one since $$\epsilon_n = \frac{\|b - Ax^i\|}{\|A\|\|x^i\| + \|b\|} \leq \frac{\|b - Ax^i\|}{\|b\|} \tag{3}$$ However condition (3) is too conservative, and as consequence the algorithm will perform too many iterations than required to obtain reasonably accurate solution.

Notice, that you can build a good stopping rule directly from (1) as $$\|b - Ax^i\|_\infty \leq (\|A\|_\infty\|x^i\|_\infty + \|b\|_\infty)\times\epsilon \tag{4}$$ since $\|A\|_\infty$ can be easily estimated. Then the tolerance $\epsilon$ is just the normwise backward error.

Since generally componentwise error bounds are much better, than normwise bounds, a stopping rule based on (2) would be better than (4). Unfortunately this would require calculating $|A||x^i|$, which is not always possible.

$\endgroup$
  • $\begingroup$ So basically the best stopping rules here are the last one I present and also the ones presented here in (4)? $\endgroup$ – JP McCarthy Sep 19 '17 at 7:56
  • 1
    $\begingroup$ Yes. When your last rule show convergence, then it is guaranteed, that desired accuracy is reached. But it is possible, that this rule will show no convergence for accurate enough solution (even for the most accurate solution, that can be obtained by any method). Rule (4) is more precise and robust. $\endgroup$ – Pawel Kowal Sep 19 '17 at 15:10
  • $\begingroup$ What are $|A|$ and $|x_i|$ in your componentwise error analysis? It is not clear from context what those are. $\endgroup$ – shuhalo May 1 at 4:55
  • 1
    $\begingroup$ For a matrix $A$ with elements $a_{ij}$, $|A|$ is a matrix with elements $|a_{ij}|$. $\endgroup$ – Pawel Kowal May 1 at 12:11
  • $\begingroup$ So the problem is that $|A| \cdot |x|$ might be too expensive to compute, right? This looks like it can be implemented straightforward though. $\endgroup$ – shuhalo Sep 18 at 12:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.