octave_packages/general-1.3.1/majle.m

   1 ## Copyright (c) 2010 Andrew V. Knyazev <andrew.knyazev@ucdenver.edu>
   2 ## Copyright (c) 2010 Merico .E. Argentati <Merico.Argentati@ucdenver.edu>
   3 ## All rights reserved.
   4 ##
   5 ## Redistribution and use in source and binary forms, with or without
   6 ## modification, are permitted provided that the following conditions are met:
   7 ##
   8 ##     1 Redistributions of source code must retain the above copyright notice,
   9 ##       this list of conditions and the following disclaimer.
  10 ##     2 Redistributions in binary form must reproduce the above copyright
  11 ##       notice, this list of conditions and the following disclaimer in the
  12 ##       documentation and/or other materials provided with the distribution.
  13 ##     3 Neither the name of the author nor the names of its contributors may be
  14 ##       used to endorse or promote products derived from this software without
  15 ##       specific prior written permission.
  16 ##
  17 ## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
  18 ## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  19 ## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  20 ## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
  21 ## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  22 ## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
  23 ## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  24 ## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
  25 ## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  26 ## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  27
  28 %MAJLE  (Weak) Majorization check
  29 %    S = MAJLE(X,Y) checks if the real part of X is (weakly) majorized by
  30 %    the real part of Y, where X and Y must be numeric (full or sparse)
  31 %    arrays. It returns S=0, if there is no weak majorization of X by Y,
  32 %    S=1, if there is a weak majorization of X by Y, or S=2, if there is a
  33 %    strong majorization of X by Y. The shapes of X and Y are ignored.
  34 %    NUMEL(X) and NUMEL(Y) may be different, in which case one of them is
  35 %    appended with zeros to match the sizes with the other and, in case of
  36 %    any negative components, a special warning is issued.
  37 %
  38 %    S = MAJLE(X,Y,MAJLETOL) allows in addition to specify the tolerance in
  39 %    all inequalities [S,Z] = MAJLE(X,Y,MAJLETOL) also outputs a row vector
  40 %    Z, which appears in the definition of the (weak) majorization. In the
  41 %    traditional case, where the real vectors X and Y are of the same size,
  42 %    Z = CUMSUM(SORT(Y,'descend')-SORT(X,'descend')). Here, X is weakly
  43 %    majorized by Y, if MIN(Z)>0, and strongly majorized if MIN(Z)=0, see
  44 %    http://en.wikipedia.org/wiki/Majorization
  45 %
  46 %    The value of MAJLETOL depends on how X and Y have been computed, i.e.,
  47 %    on what the level of error in X or Y is. A good minimal starting point
  48 %    should be MAJLETOL=eps*MAX(NUMEL(X),NUMEL(Y)). The default is 0.
  49 %
  50 %    % Examples:
  51 %    x = [2 2 2]; y = [1 2 3]; s = majle(x,y)
  52 %    % returns the value 2.
  53 %    x = [2 2 2]; y = [1 2 4]; s = majle(x,y)
  54 %    % returns the value 1.
  55 %    x = [2 2 2]; y = [1 2 2]; s = majle(x,y)
  56 %    % returns the value 0.
  57 %    x = [2 2 2]; y = [1 2 2]; [s,z] = majle(x,y)
  58 %    % also returns the vector z = [ 0 0 -1].
  59 %    x = [2 2 2]; y = [1 2 2]; s = majle(x,y,1)
  60 %    % returns the value 2.
  61 %    x = [2 2]; y = [1 2 2]; s = majle(x,y)
  62 %    % returns the value 1 and warns on tailing with zeros
  63 %    x = [2 2]; y = [-1 2 2]; s = majle(x,y)
  64 %    % returns the value 0 and gives two warnings on tailing with zeros
  65 %    x = [2 -inf]; y = [4 inf]; [s,z] = majle(x,y)
  66 %    % returns s = 1 and z = [Inf   Inf].
  67 %    x = [2 inf]; y = [4 inf]; [s,z] = majle(x,y)
  68 %    % returns  s = 1 and z = [NaN NaN] and a warning on NaNs in z.
  69 %    x=speye(2); y=sparse([0 2; -1 1]); s = majle(x,y)
  70 %    % returns the value 2.
  71 %    x = [2 2; 2 2]; y = [1 3 4]; [s,z] = majle(x,y) %and
  72 %    x = [2 2; 2 2]+i; y = [1 3 4]-2*i; [s,z] = majle(x,y)
  73 %    % both return s = 2 and z = [2 3 2 0].
  74 %    x = [1 1 1 1 0]; y = [1 1 1 1 1 0 0]'; s = majle(x,y)
  75 %    % returns the value 1 and warns on tailing with zeros
  76 %
  77 %    % One can use this function to check numerically the validity of the
  78 %    Schur-Horn,Lidskii-Mirsky-Wielandt, and Gelfand-Naimark theorems:
  79 %    clear all; n=100; majleTol=n*n*eps;
  80 %    A = randn(n,n); A = A'+A; eA = -sort(-eig(A)); dA = diag(A);
  81 %    majle(dA,eA,majleTol) % returns the value 2
  82 %    % which is the Schur-Horn theorem; and
  83 %    B=randn(n,n); B=B'+B; eB=-sort(-eig(B));
  84 %    eAmB=-sort(-eig(A-B));
  85 %    majle(eA-eB,eAmB,majleTol) % returns the value 2
  86 %    % which is the Lidskii-Mirsky-Wielandt theorem; finally
  87 %    A = randn(n,n); sA = -sort(-svd(A));
  88 %    B = randn(n,n); sB = -sort(-svd(B));
  89 %    sAB = -sort(-svd(A*B));
  90 %    majle(log2(sAB)-log2(sA), log2(sB), majleTol) % retuns the value 2
  91 %    majle(log2(sAB)-log2(sB), log2(sA), majleTol) % retuns the value 2
  92 %    % which are the log versions of the Gelfand-Naimark theorems
  93
  94 %   Tested in MATLAB 7.9.0.529 (R2009b) and Octave 3.2.3.
  95 function [s,z]=majle(x,y,majleTol)
  96
  97     if (nargin < 2)
  98         error('MAJORIZATION:majle:NotEnoughInputs',...
  99             'Not enough input arguments.');
 100     end
 101     if (nargin > 3)
 102         error('MAJORIZATION:majle:TooManyInputs',...
 103             'Too many input arguments.');
 104     end
 105     if (nargout > 2)
 106         error('MAJORIZATION:majle:TooManyOutputs',...
 107             'Too many output arguments.');
 108     end
 109
 110     % Assign default values to unspecified parameters
 111     if (nargin == 2)
 112         majleTol = 0;
 113     end
 114
 115     % transform into real (row) vectors
 116     x=real(x); xc=reshape(x,1,numel(x)); clear x;
 117     y=real(y); yc=reshape(y,1,numel(y)); clear y;
 118
 119     % sort both vectors in descending order
 120     xc=-sort(-xc); yc=-sort(-yc);
 121
 122     % tail with zeros the shorter vector to make vectors of the same length
 123     if size(xc,2)~=size(yc,2)
 124         checkForNegative = (xc(end) < -majleTol) || (yc(end) < -majleTol);
 125         warning('MAJORIZATION:majle:ResizeVectors', ...
 126             'The input vectors have different sizes. Tailing with zeros.');
 127         yc=[yc zeros(size(xc,2)-size(yc,2),1)'];
 128         xc=[xc zeros(size(yc,2)-size(xc,2),1)'];
 129         % but warn if negative
 130         if checkForNegative
 131             warning('MAJORIZATION:majle:ResizeNegativeVectors', ...
 132                 sprintf('%s%s\n%s\n%s', ...
 133                 'At least one of the input vectors ',...
 134                 'has negative components.',...
 135                 '         Tailing with zeros is most likely senseless.',...
 136                 '         Make sure that you know what you are doing.'));
 137             % sort again both vectors in descending order
 138             xc=-sort(-xc); yc=-sort(-yc);
 139         end
 140     end
 141     z=cumsum(yc-xc);
 142
 143     %check for NaNs in z
 144     if any(isnan(z))
 145         warning('MAJORIZATION:majle:NaNsInComparisons', ...
 146             sprintf('%s%s\n%s\n%s', ...
 147             'At least one of the input vectors ',...
 148             'includes -Inf, Inf, or NaN components.',...
 149             '         Some comparisons could not be made. ',...
 150             '         Make sure that you know what you are doing.'));
 151     end
 152
 153     if min(z) < -majleTol
 154         s=0;  % no majorization
 155     elseif abs(z(end)) <= majleTol
 156         s=2;   % strong majorization
 157     else
 158         s=1; % weak majorization
 159     end
 160 endfunction