maka89 / matsq

Iterative method to find square root of a matrix

MatSq

Iterative method to find square root of a matrix

• More numerically stable than babylonian method.

• Can make use of initial guess, as opposed to Denman beavers.

• Only matrix multiplications and suc. No inverses in the iteration.

  Uses Newton's method. Step is calculated using Conjugate Gradient, truncated to one or a few iterations.

Theory

Want to solve [N,N] matrix-equation:

$$ F(X) = XX-Y = 0$$

Newtons method:

$$ F'(X) (X_{k+1}-X_k) = - f(X_k)$$

$F'(X)$ is a sparse [N,N,N,N] tensor. But we have a nice expression for the jacobian-vector product (Or jacobian-matrix product in this case?).

$$ F'(X_k)H = XH+HX. $$

This is all we need if we choose to use a matrix-free method like Conjugate Gradient to get the step in Newton's method. This inner iteration can be truncated after only a few iterations.

See N.Higham - Newton's Method for the Matrix Square Root