InducedMatrixNorm

From QETLAB
Jump to: navigation, search
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, which starts with a random vector and performs a local optimization based on that starting vector.

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 (PQ)-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.

  1. %%  INDUCEDMATRIXNORM    Computes a lower bound of the induced p->q norm of a matrix
  2. %   This function has two required arguments:
  3. %     X: a matrix
  4. %     P: a real number >= 1 or Inf
  5. %
  6. %   NRM = InducedMatrixNorm(X,P) is a lower bound of the induced P-norm of
  7. %   the matrix X. This estimate of the norm is computed via a randomized
  8. %   algorithm, and thus running this function multiple times may produce
  9. %   different lower bounds.
  10. %
  11. %   This function has three optional input arguments:
  12. %     Q: a real number >= 1 or Inf (by default, Q = P)
  13. %     TOL: numerical tolerance used to determine when the algorithm stops
  14. %          running (default sqrt(eps))
  15. %     V0: a vector that acts as a starting point for the randomized
  16. %         algorithm (default is randomly-generated)
  17. %
  18. %   This function has one optional output argument:
  19. %     V: the best right-multiplication vector that was found (i.e., the
  20. %        vector that maximizes norm(X*V,Q) subject to norm(V,P) = 1).
  21. %   
  22. %   [NRM,V] = InducedMatrixNorm(X,P,Q,TOL,V0) is a lower bound of the
  23. %   induced P->Q-norm of the matrix X. This estimate of the norm is
  24. %   computed via a randomized algorithm, and thus running this function
  25. %   multiple times (with different V0) may produce different lower bounds.
  26. %   Smaller values of TOL give better numerical precision, but increase the
  27. %   running time of the algorithm.
  28. %
  29. %   URL: http://www.qetlab.com/InducedMatrixNorm
  30.  
  31. %   requires: opt_args.m
  32. %             
  33. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)
  34. %   package: QETLAB
  35. %   last updated: January 8, 2016
  36.  
  37. function [nrm,v] = InducedMatrixNorm(X,p,varargin)
  38.  
  39. [n,m] = size(X); % size of the matrix
  40.  
  41. % Set optional argument defaults: q=p, tol=10^-8, v0=-1 (randomly-generated v0)
  42. [q,tol,v0] = opt_args({ p, sqrt(eps), -1 },varargin{:});
  43.  
  44. % Quickly compute in some special cases.
  45. if(p == 1 && q == 1)
  46.     [nrm,ind] = max(sum(abs(X),1)); % norm is max abs column sum
  47.     v = zeros(m,1);
  48.     v(ind) = 1;
  49.     return
  50. elseif((p == 2 || strcmpi(p,'fro') == 1) && (q == 2 || strcmpi(q,'fro') == 1))
  51.     [~,nrm,v] = svds(X,1); % norm is largest singular value
  52.     return
  53. elseif(p == Inf && q == Inf)
  54.     nrm = max(sum(abs(X),2)); % norm is max abs row sum
  55.     v = ones(m,1);
  56.     return
  57. end
  58.  
  59. % In all other cases, we iterate to compute the induced matrix norm.
  60.  
  61. % If the user specified a starting guess v0, parse it; otherwise randomly
  62. % generate one.
  63. randv0 = 1;
  64. if(max(size(v0)) > 1)
  65.     v0 = v0(:); % make sure it's a column vector
  66.     if(length(v0) ~= m)
  67.         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.');
  68.     else
  69.         randv0 = 0;
  70.     end
  71. end
  72. if randv0 % generate a random starting vector v0, if appropriate
  73.     v = randn(m,1);
  74.     if(~isreal(X)) % only add imaginary part to v if X is not real (just to make output prettier)
  75.         v = v + 1i*randn(m,1);
  76.     end
  77. else
  78.     v = v0;
  79. end
  80. v = v/norm(v,p); % normalize the starting vector
  81.  
  82. % Preparation is done; now do the actual iteration.
  83. it_err = 2*tol+1;
  84. nrm = norm(X*v,q);
  85.  
  86. while it_err > tol
  87.     % First, find the best left vector w, keeping the right vector v fixed.
  88.     w = X*v;
  89.     if(q == Inf)
  90.         [~,ind] = max(abs(w));
  91.         w = zeros(m,1);
  92.         w(ind) = 1;
  93.     else
  94.         wabs = abs(w); % split w into its phases and magnitudes
  95.         wph = w./wabs;
  96.         wph(isnan(wph)) = 1; % take care of division by 0 in previous line
  97.  
  98.         wabs = wabs/max(wabs); % pre-process in this way first for numerical reasons
  99.         wabs = wabs.^(q-1); % this is the equality condition from Holder's inequality
  100.         w = wph.*wabs/norm(wabs,q/(q-1));
  101.     end
  102.  
  103.     % Next, find the best right vector v, keeping the left vector w fixed.
  104.     v = w'*X;
  105.     if(p == 1)
  106.         [~,ind] = max(abs(v));
  107.         v = zeros(n,1);
  108.         v(ind) = 1;
  109.     else
  110.         vabs = abs(v); % split v into its phases and magnitudes
  111.         vph = v./vabs;
  112.         vph(isnan(vph)) = 1; % take care of division by 0 in previous line  
  113.  
  114.         vabs = vabs'/max(vabs); % pre-process in this way first for numerical reasons
  115.         vabs = vabs.^(1/(p-1)); % this is the equality condition from Holder's inequality
  116.         v = vph'.*vabs/norm(vabs,p);
  117.     end
  118.  
  119.     % Check to see if we made any progress; if so, keep iterating.
  120.     new_nrm = norm(X*v,q);
  121.     it_err = abs(new_nrm - nrm);
  122.     nrm = new_nrm;
  123. end