Mercurial > matrix-functions

comparison mftoolbox/signm_newton.m @ 0:8f23314345f4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Create local repository for matrix toolboxes. Step #0 done.
author Antonio Pino Robles <data.script93@gmail.com>
date Wed, 06 May 2015 14:56:53 +0200
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:8f23314345f4
1 function [X,k] = signm_newton(A,scal)
2 %SIGNM_NEWTON Matrix sign function by Newton iteration.
3 % [S,k] = SIGNM_NEWTON(A,SCAL) computes the matrix sign
4 % function S of A using the scaled Newton iteration, with
5 % scale factors specified by SCAL:
6 % SCAL = 0: no scaling.
7 % SCAL = 1: determinantal scaling (default).
8 % SCAL = 2: spectral scaling.
9 % SCAL = 3: norm scaling.
10 % k is the number of iterations.
11
12 if nargin < 2, scal = 1; end
13
14 tol = mft_tolerance(A);
15 accel_tol = 1e-2; % Precise value not important.
16 need_accel = 1;
17 n = length(A);
18 maxit = 16;
19
20 X = A;
21 reldiff = inf;
22
23 for k = 1:maxit
24
25 Xold = X;
26 Xinv = inv(X);
27 if need_accel && scal > 0
28 % In practice should estimate spectral radius and 2-norms;
29 % here they are computed exactly.
30 switch scal
31 case 1
32 g = abs(det(X))^(-1/n);
33 case 2
34 s1 = max(abs(eig(Xinv)));
35 s2 = max(abs(eig(X)));
36 g = sqrt(s1/s2);
37 case 3
38 g = sqrt( norm(Xinv) / (norm(X)) );
39 end
40 X = g*X; Xinv = Xinv/g;
41 end
42 X = 0.5*(X + Xinv);
43 diff_F = norm(X-Xold,'fro');
44 reldiff_old = reldiff;
45 reldiff = diff_F/norm(X,'fro');
46
47 if need_accel && (reldiff < accel_tol), need_accel = false; end
48 cged = (diff_F <= sqrt( tol*norm(X)/norm(Xinv) ) || ...
49 reldiff > reldiff_old/2 && ~need_accel);
50 if cged, return, end
51
52 end
53 error('Not converged after %2.0f iterations', maxit)