Mercurial > matrix-functions
view mftoolbox/polar_newton.m @ 4:0aae72032c25 draft
Added files described in changeset d24a00dabdc2.
author | Antonio Pino Robles <data.script93@gmail.com> |
---|---|
date | Fri, 15 May 2015 20:30:21 +0200 |
parents | 8f23314345f4 |
children |
line wrap: on
line source
function [U,H,k] = polar_newton(A) %POLAR_NEWTON Polar decomposition by scaled Newton iteration. % [U,H,k] = POLAR_NEWTON(A), where the matrix A is square and % nonsingular, computes a unitary U and a Hermitian positive % definite H such that A = U*H. k is the number of iterations. % A Newton iteration with acceleration parameters is used. accel_tol = 1e-2; % Precise value not important. need_accel = 1; tol = mft_tolerance(A); X = A; reldiff = inf; maxit = 16; for k = 1:maxit Xold = X; Xinv = inv(X); if need_accel g = ( norm(Xinv,1)*norm(Xinv,inf) / (norm(X,1)*norm(X,inf)) )^(1/4); else g = 1; end X = 0.5*(g*X + Xinv'/g); reldiff_old = reldiff; diff_F = norm(X-Xold,'fro'); reldiff = diff_F/norm(X,'fro'); if need_accel && (reldiff < accel_tol), need_accel = false; end cged = (diff_F <= sqrt(tol)) || (reldiff > reldiff_old/2 && ~need_accel); if cged, break, end if k == maxit, error('Not converged after %2.0f iterations', maxit), end end U = X; if nargout >= 2 H = U'*A; H = (H + H')/2; % Force Hermitian by taking nearest Hermitian matrix. end