Mercurial > matrix-functions
view toolbox/pnorm.m @ 3:d24a00dabdc2 draft
Added first test and Makefile after figuring out how to!!! Also removed
fontconfig rubbish created by a bug in my GNU/linux system.
Note: review E.B. Davies's paper "Approximate Diagonalization" for tricky
diagonalization problems.
author | Antonio Pino Robles <data.script93@gmail.com> |
---|---|
date | Fri, 15 May 2015 20:24:25 +0200 |
parents | 8f23314345f4 |
children |
line wrap: on
line source
function [est, x, k] = pnorm(A, p, tol, noprint) %PNORM Estimate of matrix p-norm (1 <= p <= inf). % [EST, x, k] = PNORM(A, p, TOL) estimates the Holder p-norm of a % matrix A, using the p-norm power method with a specially % chosen starting vector. % TOL is a relative convergence tolerance (default 1E-4). % Returned are the norm estimate EST (which is a lower bound for the % exact p-norm), the corresponding approximate maximizing vector x, % and the number of power method iterations k. % A nonzero fourth argument causes trace output to the screen. % If A is a vector, this routine simply returns NORM(A, p). % % See also NORM, NORMEST. % Note: The estimate is exact for p = 1, but is not always exact for % p = 2 or p = inf. Code could be added to treat p = 2 and p = inf % separately. % % Calls DUAL and SEQA. % % Reference: % N.J. Higham, Estimating the matrix p-norm, % Numer. Math., 62 (1992), pp. 539-555. if nargin < 2, error('Must specify norm via second parameter.'), end [m,n] = size(A); if min(m,n) == 1, est = norm(A,p); return, end if nargin < 4, noprint = 0; end if nargin < 3, tol = 1e-4; end % Stage I. Use Algorithm OSE to get starting vector x for power method. % Form y = B*x, at each stage choosing x(k) = c and scaling previous % x(k+1:n) by s, where norm([c s],p)=1. sm = 9; % Number of samples. y = zeros(m,1); x = zeros(n,1); for k=1:n if k == 1 c = 1; s = 0; else W = [A(:,k) y]; if p == 2 % Special case. Solve exactly for 2-norm. [U,S,V] = svd(full(W)); c = V(1,1); s = V(2,1); else fopt = 0; for th=seqa(0,pi,sm) c1 = cos(th); s1 = sin(th); nrm = norm([c1 s1],p); c1 = c1/nrm; s1 = s1/nrm; % [c1 s1] has unit p-norm. f = norm( W*[c1 s1]', p ); if f > fopt fopt = f; c = c1; s = s1; end end end end x(k) = c; y = x(k)*A(:,k) + s*y; if k > 1, x(1:k-1) = s*x(1:k-1); end end est = norm(y,p); if noprint, fprintf('Alg OSE: %9.4e\n', est), end % Stage II. Apply Algorithm PM (the power method). q = dual(p); k = 1; while 1 y = A*x; est_old = est; est = norm(y,p); z = A' * dual(y,p); if noprint fprintf('%2.0f: norm(y) = %9.4e, norm(z) = %9.4e', ... k, norm(y,p), norm(z,q)) fprintf(' rel_incr(est) = %9.4e\n', (est-est_old)/est) end if ( norm(z,q) <= z'*x | abs(est-est_old)/est <= tol ) & k > 1 return end x = dual(z,q); k = k + 1; end