2
$\begingroup$

This is a follow-up to this answer.

Suppose you have a possibly very ill-conditioned matrix $A$, and you compute its inverse with LU/GE to get $X_{\text{lu}}\approx A^{-1}$.

The Newton-Raphson iteration $$ f: X \mapsto 2X-XAX $$ generally converges to the matrix inverse of $A$.

Does computing applying the Newton-Raphson iteration to $X_{\text{lu}}$, that is, computing $f(X_{\text{lu}})$, $f(f(X_{\text{lu}}))$ make the resulting matrix inverse more accurate, or less accurate?

I was unable to find a reference directly, so I'm asking if someone knows for sure.

The absolute error will be proportional to $\|(f')^{-1}\|$ (e.g., Accuracy and Stability of Numerical Algorithms, Chapter 25), but $f'$ vanishes at $X=A^{-1}$. So presumably this means that applying $f$ could only help when the errors came from the growth factors in the LU factorization step? But doesn't that mean that applying $f$ hurts at other times since it relies on computing a root of a function with a vanishing Jacobian, which would make it less accurate than LU? For example, in the one-dimensional case with a root of multiplicity 2, Newton's method would have error $O(\sqrt{\epsilon})$, but LU/GE always has error $O(\rho \epsilon)$.

I tried it on all $16\times 16$ matrices in MatrixDepot.jl, and it seems like it almost always hurts accuracy, never helps by more than $\frac13$, and often fails. (It's not very clear what matrices to test.)

module MatrixInverseNewtonRaphson

using MatrixDepot

F(A, X) = 2X - X*A*X
R(A, X) = norm(A*X - eye(A), 1)

function go(name, m)
    A = full(matrixdepot(name, m))
    B = inv(A)
    B1 = F(A, B)
    X = 0.5A' / norm(A*A', 2)
    for i in 1:2m; X = F(A, X); end
    (R(A, B1) / max(eps(eltype(A)), R(A, B)), R(A, B), R(A, X), name)
end

function main(m=16)
    e = []
    for name in matrixdepot("all")
        try
            push!(e, go(name, m))
        catch err
            println(err)
        end
    end
    sort(e)
end

end
$\endgroup$
  • $\begingroup$ to me it is not clear whether (or why) the convergence is ensured. In general the Newton method is used with start point $X_0 = \alpha A^{T}$, with $\alpha$ choses such that the the eigenvalues of the residual $R_0$ are in module < 1 ,with $R_n = I -AX_n$. So $\lim_{n \rightarrow \infty} (R_0)^{n} = 0$. (See pdf theorem 2.1 for detail). Using $X_0 = X_{lu}$ is this, i.e. the convergence, guaranteed? (If you scale $X_{lu}$ does the situation improve?) $\endgroup$ – Mauro Vanzetto Oct 30 '16 at 18:45
  • $\begingroup$ @MauroVanzetto $X_{\text{lu}}$ is so close to the true root that I think convergence should be fine for typical matrices. (Although MatrixDepot includes pathological matrices.) Typically convergence conditions for Newton methods would be over-strict to guarantee convergence. $\endgroup$ – Kirill Oct 31 '16 at 17:49

Your Answer

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

Browse other questions tagged or ask your own question.