Mercurial > matrix-functions

comparison mftoolbox/funm_condest_fro.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 c = funm_condest_fro(A,fun,fun_frechet,its,flag,varargin)
2 %FUNM_CONDEST_FRO Estimate of Frobenius norm condition number of matrix function.
3 % C = FUNM_CONDEST_FRO(A,FUN,FUN_FRECHET,ITS,FLAG) produces an estimate of
4 % the Frobenius norm relative condition number of function FUN at
5 % the matrix A. FUN and FUN_FRECHET are function handles:
6 % - FUN(A) evaluates the function at the matrix A.
7 % - If FLAG == 0 (default)
8 % FUN_FRECHET(B,E) evaluates the Frechet derivative at B
9 % in the direction E;
10 % if FLAG == 1
11 % - FUN_FRECHET('notransp',E) evaluates the
12 % Frechet derivative at A in the direction E.
13 % - FUN_FRECHET('transp',E) evaluates the
14 % Frechet derivative at A' in the direction E.
15 % If FUN_FRECHET is empty then the Frechet derivative is approximated
16 % by finite differences. More reliable results are obtained when
17 % FUN_FRECHET is supplied.
18 % The power method is used, with a random starting matrix,
19 % so the approximation can be different each time.
20 % ITS iterations (default 6) are used.
21 % C = FUNM_COND_EST_FRO(A,FUN,FUN_FRECHET,ITS,FLAG,P1,P2,...)
22 % passes extra inputs P1,P2,... to FUN and FUN_FRECHET.
23 % Note: this function makes an assumption on the adjoint of the
24 % Frechet derivative that, for f having a power series expansion,
25 % is equivalent to the series having real coefficients.
26
27 if nargin < 5 || isempty(flag), flag = 0; end
28 if nargin < 4 || isempty(its), its = 6; end
29 if nargin < 3 || isempty(fun_frechet), fte_diff = 1; else fte_diff = 0; end
30
31 funA = feval(fun,A,varargin{:});
32 d = sqrt( eps*norm(funA,'fro') );
33 Z = randn(size(A));
34 Znorm = 1;
35
36 factor = norm(A,'fro')/norm(funA,'fro');
37
38 for i=1:its
39
40 Z = Z/norm(Z,'fro');
41 if fte_diff
42 W = (feval(fun,A+d*Z,varargin{:}) - funA)/d;
43 else
44 if flag
45 W = feval(fun_frechet,'notransp',Z,varargin{:});
46 else
47 W = feval(fun_frechet,A,Z,varargin{:});
48 end
49 end
50
51 W = W/norm(W,'fro');
52 if fte_diff
53 Z = (feval(fun,A'+d*W,varargin{:}) - funA')/d;
54 else
55 if flag
56 Z = feval(fun_frechet,'transp',W,varargin{:});
57 else
58 Z = feval(fun_frechet,A',W,varargin{:});
59 end
60 end
61
62 Znorm = norm(Z,'fro');
63
64 end
65
66 c = Znorm*factor;