2
$\begingroup$

From what I've read, my understanding of iterative refinement is as follows (largely based on Björck's Numerical Least Squares text).

We want to solve $Ax=b$ to machine / working precision. As a first step, we compute a numerically approximate solution $x_0 :=$ x_0 using a standard decomposition / solve of $Ax=b$ (i.e., using LU or LU with partial pivoting or QR). Next we "refine" our approximate solution with some number of steps (which I'll call $K$ steps) and seek greater precision in x_k by iteratively computing residuals with high precision inner-products (let's call this extended precision $\epsilon_2$) and then feeding the high precision residual back through an update / solver scheme to compute corrections (albeit in a worse / standard precision, which I'll call $\epsilon_1$).

The intuition for this scheme comes from the observation that $$ r_0 = b-Ax_0 = b-Ax_0 - (b-Ax^*) = A(x^*-x_0) $$ which implies that if we had exact arithmetic, adding $\delta_0 = A^{-1}r_0$ would refine our original $x_0$ solution entirely. In the numerical precision case, if we improve the precision in the residual calculation, then this should "flow through" to refining $x_0$. In pseudocode this is:

for k = 1:K
    r_k = b - A * x_k     # compute k-th residual in precision eps_2
    del_k = A \ r_k       # solve for del_k in precision eps_1
    x_k = x_k + del_k     # update approximate solution
end

Björck discusses the convergence:

  1. If $A$ is very ill-conditioned, $\hat{x}$ has no correct significant digits, and the iterative scheme won't help
  2. If $A$ is not too poorly conditioned and has some accuracy, then the relative error from each iteration decreases by some constant until we reach a solution correct to machine / working precision

A couple of questions:

  1. If we start from a random $\hat{x}$ with no correct significant digits but $A$ is not too poorly conditioned, I think the above procedure should refine the solution, so why not introduce it as a method for a general starting point $\bar{x}$ that is not even an approximate solution (or is this obvious)?
  2. Using $A' = A + \lambda I$ in the intermediate step, does it make sense to solve $A' \delta_k = r_k$ with $\lambda \to 0$ as $k$ increases? That is, in the case $x_0$ comes from another technique (possibly pre-conditioned Conjugate Gradient), can we use the (regularized) decomposition of the preconditioner in our intermediate residual updates? Boyd's paper indicates "yes", but I don't understand the intuition. Further, if this works, can we find some clever way to progressively reduce the regularized decomposition so that it approaches the original problem (that is, can we update the decomposition efficiently for $\lambda \to 0$)? Motivated by p24 (p22 of pdf) of Boyd's discussion wrt monotone operators.
$\endgroup$
1
  • 1
    $\begingroup$ I fixed the link and a few minor matters. Could you elaborate a bit on your thoughts in connection with question 2? $\endgroup$ Jul 12, 2017 at 16:19

0

You must log in to answer this question.