Mercurial > matrix-functions

comparison matrixcomp/pnorm.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 [est, x, k] = pnorm(A, p, tol, prnt)
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 input 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, NORMEST1.
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.
20 %
21 % Reference:
22 % N. J. Higham, Estimating the matrix p-norm, Numer. Math.,
23 % 62 (1992), pp. 539-555.
24 % N. J. Higham, Accuracy and Stability of Numerical Algorithms,
25 % Second edition, Society for Industrial and Applied Mathematics,
26 % Philadelphia, PA, 2002; sec. 15.2.
27
28 if nargin < 2, error('Must specify norm via second parameter.'), end
29 [m,n] = size(A);
30 if min(m,n) == 1, est = norm(A,p); return, end
31
32 if nargin < 4, prnt = 0; end
33 if nargin < 3 | isempty(tol), tol = 1e-4; end
34
35 % Stage I. Use Algorithm OSE to get starting vector x for power method.
36 % Form y = B*x, at each stage choosing x(k) = c and scaling previous
37 % x(k+1:n) by s, where norm([c s],p)=1.
38
39 sm = 9; % Number of samples.
40 y = zeros(m,1); x = zeros(n,1);
41
42 for k=1:n
43
44 if k == 1
45 c = 1; s = 0;
46 else
47 W = [A(:,k) y];
48
49 if p == 2 % Special case. Solve exactly for 2-norm.
50 [U,S,V] = svd(full(W));
51 c = V(1,1); s = V(2,1);
52
53 else
54
55 fopt = 0;
56 for th=linspace(0,pi,sm)
57 c1 = cos(th); s1 = sin(th);
58 nrm = norm([c1 s1],p);
59 c1 = c1/nrm; s1 = s1/nrm; % [c1 s1] has unit p-norm.
60 f = norm( W*[c1 s1]', p );
61 if f > fopt
62 fopt = f;
63 c = c1; s = s1;
64 end
65 end
66
67 end
68 end
69
70 x(k) = c;
71 y = x(k)*A(:,k) + s*y;
72 if k > 1, x(1:k-1) = s*x(1:k-1); end
73
74 end
75
76 est = norm(y,p);
77 if prnt, fprintf('Alg OSE: %9.4e\n', est), end
78
79 % Stage II. Apply Algorithm PM (the power method).
80
81 q = dual(p);
82 k = 1;
83
84 while 1
85
86 y = A*x;
87 est_old = est;
88 est = norm(y,p);
89
90 z = A' * dual(y,p);
91
92 if prnt
93 fprintf('%2.0f: norm(y) = %9.4e, norm(z) = %9.4e', ...
94 k, norm(y,p), norm(z,q))
95 fprintf(' rel_incr(est) = %9.4e\n', (est-est_old)/est)
96 end
97
98 if ( norm(z,q) <= z'*x | abs(est-est_old)/est <= tol ) & k > 1
99 return
100 end
101
102 x = dual(z,q);
103 k = k + 1;
104
105 end