Mercurial > matrix-functions

comparison toolbox/pnorm.m @ 2:c124219d7bfa draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Re-add the 1995 toolbox after noticing the statement in the ~higham/mctoolbox/ webpage.
author Antonio Pino Robles <data.script93@gmail.com>
date Thu, 07 May 2015 18:36:24 +0200
parents 8f23314345f4
children
comparison
equal deleted inserted replaced
1:e471a92d17be 2:c124219d7bfa
1 function [est, x, k] = pnorm(A, p, tol, noprint)
2 %PNORM Estimate of matrix p-norm (1 <= p <= inf).
3 % [EST, x, k] = PNORM(A, p, TOL) estimates the Holder p-norm of a
4 % matrix A, using the p-norm power method with a specially
5 % chosen starting vector.
6 % TOL is a relative convergence tolerance (default 1E-4).
7 % Returned are the norm estimate EST (which is a lower bound for the
8 % exact p-norm), the corresponding approximate maximizing vector x,
9 % and the number of power method iterations k.
10 % A nonzero fourth argument causes trace output to the screen.
11 % If A is a vector, this routine simply returns NORM(A, p).
12 %
13 % See also NORM, NORMEST.
14
15 % Note: The estimate is exact for p = 1, but is not always exact for
16 % p = 2 or p = inf. Code could be added to treat p = 2 and p = inf
17 % separately.
18 %
19 % Calls DUAL and SEQA.
20 %
21 % Reference:
22 % N.J. Higham, Estimating the matrix p-norm,
23 % Numer. Math., 62 (1992), pp. 539-555.
24
25 if nargin < 2, error('Must specify norm via second parameter.'), end
26 [m,n] = size(A);
27 if min(m,n) == 1, est = norm(A,p); return, end
28
29 if nargin < 4, noprint = 0; end
30 if nargin < 3, tol = 1e-4; end
31
32 % Stage I. Use Algorithm OSE to get starting vector x for power method.
33 % Form y = B*x, at each stage choosing x(k) = c and scaling previous
34 % x(k+1:n) by s, where norm([c s],p)=1.
35
36 sm = 9; % Number of samples.
37 y = zeros(m,1); x = zeros(n,1);
38
39 for k=1:n
40
41 if k == 1
42 c = 1; s = 0;
43 else
44 W = [A(:,k) y];
45
46 if p == 2 % Special case. Solve exactly for 2-norm.
47 [U,S,V] = svd(full(W));
48 c = V(1,1); s = V(2,1);
49
50 else
51
52 fopt = 0;
53 for th=seqa(0,pi,sm)
54 c1 = cos(th); s1 = sin(th);
55 nrm = norm([c1 s1],p);
56 c1 = c1/nrm; s1 = s1/nrm; % [c1 s1] has unit p-norm.
57 f = norm( W*[c1 s1]', p );
58 if f > fopt
59 fopt = f;
60 c = c1; s = s1;
61 end
62 end
63
64 end
65 end
66
67 x(k) = c;
68 y = x(k)*A(:,k) + s*y;
69 if k > 1, x(1:k-1) = s*x(1:k-1); end
70
71 end
72
73 est = norm(y,p);
74 if noprint, fprintf('Alg OSE: %9.4e\n', est), end
75
76 % Stage II. Apply Algorithm PM (the power method).
77
78 q = dual(p);
79 k = 1;
80
81 while 1
82
83 y = A*x;
84 est_old = est;
85 est = norm(y,p);
86
87 z = A' * dual(y,p);
88
89 if noprint
90 fprintf('%2.0f: norm(y) = %9.4e, norm(z) = %9.4e', ...
91 k, norm(y,p), norm(z,q))
92 fprintf(' rel_incr(est) = %9.4e\n', (est-est_old)/est)
93 end
94
95 if ( norm(z,q) <= z'*x | abs(est-est_old)/est <= tol ) & k > 1
96 return
97 end
98
99 x = dual(z,q);
100 k = k + 1;
101
102 end