InducedMatrixNorm
InducedMatrixNorm | |
Computes a lower bound of the induced p→q norm of a matrix | |
Other toolboxes required | none |
---|---|
Related functions | InducedSchattenNorm |
Function category | Norms |
Usable within CVX? | no |
InducedMatrixNorm is a function that computes a randomized lower bound of the induced p→q norm of a matrix, defined as follows: \[\|B\|_{p\rightarrow q} := \max\big\{\|B\mathbf{x}\|_q : \|\mathbf{x}\|_p = 1 \big\},\] where \[\|\mathbf{x}\|_{p} := \left(\sum_i|x_i|^p\right)^{1/p}\] is the vector p-norm.
When p = q = 2, this is the usual operator norm, returned by MATLAB's built-in norm function. Similarly, when p = q = 1 or p = q = Inf, this is the maximum absolute column sum or maximum absolute row sum of the matrix, respectively, and for the matrix X it can be computed via the built-in MATLAB function norm(X,1) or norm(X,Inf). However, it most other cases this norm is hard to compute, and this function provides a randomized lower bound of it.
The lower bound is found via the algorithm described here (link to be filled in), which starts with a random vector and performs a local optimization based on that starting vector.
Contents
Syntax
- NRM = InducedMatrixNorm(X,P)
- NRM = InducedMatrixNorm(X,P,Q)
- NRM = InducedMatrixNorm(X,P,Q,TOL)
- NRM = InducedMatrixNorm(X,P,Q,TOL,V0)
- [NRM,V] = InducedMatrixNorm(X,P,Q,TOL,V0)
Argument descriptions
Input arguments
- X: A matrix to have its induced (P→Q)-norm computed.
- P: A real number ≥ 1, or Inf.
- Q (optional, default equals P): A real number ≥ 1, or Inf.
- TOL (optional, default equals sqrt(eps)): Numerical tolerance used throughout the script.
- V0 (optional, default is randomly-generated): A vector to start the numerical search from.
Output arguments
- NRM: A lower bound on the norm of X.
- V (optional): A vector with norm(V,P) = 1 such that norm(X*V,Q) = NRM (i.e., a vector that attains the local maximum that was found).
Examples
Induced norms of the identity matrix
The n-by-n identity matrix has induced p→q-norm equal to $\max\{n^{1/q - 1/p}, 1\}$, which this function finds exactly:
>> X = eye(5); >> InducedMatrixNorm(X,3,3) ans = 1 >> InducedMatrixNorm(X,3,5) ans = 1 >> InducedMatrixNorm(X,5,3) ans = 1.2394 >> 5^(1/3 - 1/5) ans = 1.2394
Source code
Click on "expand" to the right to view the MATLAB source code for this function.
%% INDUCEDMATRIXNORM Computes a lower bound of the induced p->q norm of a matrix
% This function has two required arguments:
% X: a matrix
% P: a real number >= 1 or Inf
%
% NRM = InducedMatrixNorm(X,P) is a lower bound of the induced P-norm of
% the matrix X. This estimate of the norm is computed via a randomized
% algorithm, and thus running this function multiple times may produce
% different lower bounds.
%
% This function has three optional input arguments:
% Q: a real number >= 1 or Inf (by default, Q = P)
% TOL: numerical tolerance used to determine when the algorithm stops
% running (default sqrt(eps))
% V0: a vector that acts as a starting point for the randomized
% algorithm (default is randomly-generated)
%
% This function has one optional output argument:
% V: the best right-multiplication vector that was found (i.e., the
% vector that maximizes norm(X*V,Q) subject to norm(V,P) = 1).
%
% [NRM,V] = InducedMatrixNorm(X,P,Q,TOL,V0) is a lower bound of the
% induced P->Q-norm of the matrix X. This estimate of the norm is
% computed via a randomized algorithm, and thus running this function
% multiple times (with different V0) may produce different lower bounds.
% Smaller values of TOL give better numerical precision, but increase the
% running time of the algorithm.
%
% URL: http://www.qetlab.com/InducedMatrixNorm
% requires: opt_args.m
%
% author: Nathaniel Johnston (nathaniel@njohnston.ca)
% package: QETLAB
% last updated: January 8, 2016
function [nrm,v] = InducedMatrixNorm(X,p,varargin)
[n,m] = size(X); % size of the matrix
% Set optional argument defaults: q=p, tol=10^-8, v0=-1 (randomly-generated v0)
[q,tol,v0] = opt_args({ p, sqrt(eps), -1 },varargin{:});
% Quickly compute in some special cases.
if(p == 1 && q == 1)
[nrm,ind] = max(sum(abs(X),1)); % norm is max abs column sum
v = zeros(m,1);
v(ind) = 1;
return
elseif((p == 2 || strcmpi(p,'fro') == 1) && (q == 2 || strcmpi(q,'fro') == 1))
[~,nrm,v] = svds(X,1); % norm is largest singular value
return
elseif(p == Inf && q == Inf)
nrm = max(sum(abs(X),2)); % norm is max abs row sum
v = ones(m,1);
return
end
% In all other cases, we iterate to compute the induced matrix norm.
% If the user specified a starting guess v0, parse it; otherwise randomly
% generate one.
randv0 = 1;
if(max(size(v0)) > 1)
v0 = v0(:); % make sure it's a column vector
if(length(v0) ~= m)
warning('InducedMatrixNorm:DimensionMismatch','The initial vector v0 must have length equal to the number of columns of X. Using a randomly-generated intial vector instead.');
else
randv0 = 0;
end
end
if randv0 % generate a random starting vector v0, if appropriate
v = randn(m,1);
if(~isreal(X)) % only add imaginary part to v if X is not real (just to make output prettier)
v = v + 1i*randn(m,1);
end
else
v = v0;
end
v = v/norm(v,p); % normalize the starting vector
% Preparation is done; now do the actual iteration.
it_err = 2*tol+1;
nrm = norm(X*v,q);
while it_err > tol
% First, find the best left vector w, keeping the right vector v fixed.
w = X*v;
if(q == Inf)
[~,ind] = max(abs(w));
w = zeros(m,1);
w(ind) = 1;
else
wabs = abs(w); % split w into its phases and magnitudes
wph = w./wabs;
wph(isnan(wph)) = 1; % take care of division by 0 in previous line
wabs = wabs/max(wabs); % pre-process in this way first for numerical reasons
wabs = wabs.^(q-1); % this is the equality condition from Holder's inequality
w = wph.*wabs/norm(wabs,q/(q-1));
end
% Next, find the best right vector v, keeping the left vector w fixed.
v = w'*X;
if(p == 1)
[~,ind] = max(abs(v));
v = zeros(n,1);
v(ind) = 1;
else
vabs = abs(v); % split v into its phases and magnitudes
vph = v./vabs;
vph(isnan(vph)) = 1; % take care of division by 0 in previous line
vabs = vabs'/max(vabs); % pre-process in this way first for numerical reasons
vabs = vabs.^(1/(p-1)); % this is the equality condition from Holder's inequality
v = vph'.*vabs/norm(vabs,p);
end
% Check to see if we made any progress; if so, keep iterating.
new_nrm = norm(X*v,q);
it_err = abs(new_nrm - nrm);
nrm = new_nrm;
end