On 5/20/06, Joerg Arndt <arndt@jjj.de> wrote:
... Write B := B( 1 + (1-AB) ), (second order) rewrite as B(1+Y) (where Y is small), n-th order iteration is the B(1+Y+Y^2+...+Y^(n-1))
At least with the scalar iteration, never expand as you lose the small terms Y and control of convergence.
The iterations is a special case (a=1) of the division-less iteration for the inverse a-th root. (section 15.3 "Iterations for the inverse a-th root", pp.408ff of the fxtbook gives the general form).
Using the third order iteration B(1+Y+Y^2) will likely give a speedup. ...
Um... But of course, fourth-order would be no faster. Nice that the binomial theorem still works here!
search for "Nicholas J. Higham" http://www.ma.man.ac.uk/~higham/
Some papers of interest are
\jjbibitem{apprmatlog}{Sheung Hun Cheng, Nicholas J.\ Higham, Charles S.\ Kenney, Alan J.\ Laub: \jjbibtitle{Approximating the logarithm of a matrix to specified accuracy}, SIAM Journal on Matrix Analysis and Applications, vol.22, no.4, pp.1112-1125, 2001.} % http://www.ee.ucla.edu/faculty/bios/laub.htm ...
Thanks for the references --- lots of interesting stuff there! Shame that they all posit that matrix inversion is on tap, though. By the way, I hadn't come across the "matrix sign" function before. The iteration given for sign(A) is B := A; B := (B + B^{-1})/2 which raises a couple of questions in my mind. Firstly, is it not possible to involve A in the iteration, and so avoid the computation drifting as a result of rounding error or ill-conditioning? Secondly, as it stands the iteration is inapplicable to rectangular matrices [or singular matrices, besides failing for iimaginary eigenvalues]. Suppose instead we modify it to the more amenable B := A; B := (B + B^{+T})/2 (average B with its transposed pseudo-inverse). For real nonsingular square A, the result appears to converge to an orthogonal matrix B such that A = B C with C symmetric: the "polar decomposition" of A. Presumably this must be well-known; but oddly it doesn't seem to be mentioned in in Higham et al ...
To use it for matrix inversion in general, we'd need to construct an initial approximation B for which d < 1. This unfortunately looks nontrivial. [In the original real number context, it was worked around using the fact that binary representation gives an easy approx. base 2 log; and for reduced a ~ 1 we can just choose b = 1.]
It occurred to me afterwards that what I'm really asking for here is an initial B such that A B is diagonally dominant --- exactly the same as required to make Gauss-Seidel applicable! Now I wonder which box I packed "Numerical Recipies" in ... Thanks again, Fred Lunnon