octave_packages/m/linear-algebra/logm.m

   1 ## Copyright (C) 2008-2012 N.J. Higham
   2 ## Copyright (C) 2010 Richard T. Guy <guyrt7@wfu.edu>
   3 ## Copyright (C) 2010 Marco Caliari <marco.caliari@univr.it>
   4 ##
   5 ## This file is part of Octave.
   6 ##
   7 ## Octave is free software; you can redistribute it and/or modify it
   8 ## under the terms of the GNU General Public License as published by
   9 ## the Free Software Foundation; either version 3 of the License, or (at
  10 ## your option) any later version.
  11 ##
  12 ## Octave is distributed in the hope that it will be useful, but
  13 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  14 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15 ## General Public License for more details.
  16 ##
  17 ## You should have received a copy of the GNU General Public License
  18 ## along with Octave; see the file COPYING.  If not, see
  19 ## <http://www.gnu.org/licenses/>.
  20
  21 ## -*- texinfo -*-
  22 ## @deftypefn  {Function File} {@var{s} =} logm (@var{A})
  23 ## @deftypefnx {Function File} {@var{s} =} logm (@var{A}, @var{opt_iters})
  24 ## @deftypefnx {Function File} {[@var{s}, @var{iters}] =} logm (@dots{})
  25 ## Compute the matrix logarithm of the square matrix @var{A}.  The
  26 ## implementation utilizes a Pad@'e approximant and the identity
  27 ##
  28 ## @example
  29 ## logm (@var{A}) = 2^k * logm (@var{A}^(1 / 2^k))
  30 ## @end example
  31 ##
  32 ## The optional argument @var{opt_iters} is the maximum number of square roots
  33 ## to compute and defaults to 100.  The optional output @var{iters} is the
  34 ## number of square roots actually computed.
  35 ## @seealso{expm, sqrtm}
  36 ## @end deftypefn
  37
  38 ## Reference: N. J. Higham, Functions of Matrices: Theory and Computation
  39 ##            (SIAM, 2008.)
  40 ##
  41
  42 function [s, iters] = logm (A, opt_iters = 100)
  43
  44   if (nargin == 0 || nargin > 2)
  45     print_usage ();
  46   endif
  47
  48   if (! issquare (A))
  49     error ("logm: A must be a square matrix");
  50   endif
  51
  52   if (isscalar (A))
  53     s = log (A);
  54     return;
  55   elseif (strfind (typeinfo (A), "diagonal matrix"))
  56     s = diag (log (diag (A)));
  57     return;
  58   endif
  59
  60   [u, s] = schur (A);
  61
  62   if (isreal (A))
  63     [u, s] = rsf2csf (u, s);
  64   endif
  65
  66   eigv = diag (s);
  67   if (any (eigv < 0))
  68     warning ("Octave:logm:non-principal",
  69              "logm: principal matrix logarithm is not defined for matrices with negative eigenvalues; computing non-principal logarithm");
  70   endif
  71
  72   real_eig = all (eigv >= 0);
  73
  74   k = 0;
  75   ## Algorithm 11.9 in "Function of matrices", by N. Higham
  76   theta = [0, 0, 1.61e-2, 5.38e-2, 1.13e-1, 1.86e-1, 2.6429608311114350e-1];
  77   p = 0;
  78   m = 7;
  79   while (k < opt_iters)
  80     tau = norm (s - eye (size (s)),1);
  81     if (tau <= theta (7))
  82       p = p + 1;
  83       j(1) = find (tau <= theta, 1);
  84       j(2) = find (tau / 2 <= theta, 1);
  85       if (j(1) - j(2) <= 1 || p == 2)
  86         m = j(1);
  87         break
  88       endif
  89     endif
  90     k = k + 1;
  91     s = sqrtm (s);
  92   endwhile
  93
  94   if (k >= opt_iters)
  95     warning ("logm: maximum number of square roots exceeded; results may still be accurate");
  96   endif
  97
  98   s = s - eye (size (s));
  99
 100   if (m > 1)
 101     s = logm_pade_pf (s, m);
 102   endif
 103
 104   s = 2^k * u * s * u';
 105
 106   ## Remove small complex values (O(eps)) which may have entered calculation
 107   if (real_eig && isreal(A))
 108     s = real (s);
 109   endif
 110
 111   if (nargout == 2)
 112     iters = k;
 113   endif
 114
 115 endfunction
 116
 117 ################## ANCILLARY FUNCTIONS ################################
 118 ######  Taken from the mfttoolbox (GPL 3) by D. Higham.
 119 ######  Reference:
 120 ######      D. Higham, Functions of Matrices: Theory and Computation
 121 ######      (SIAM, 2008.).
 122 #######################################################################
 123
 124 ##LOGM_PADE_PF   Evaluate Pade approximant to matrix log by partial fractions.
 125 ##   Y = LOGM_PADE_PF(A,M) evaluates the [M/M] Pade approximation to
 126 ##   LOG(EYE(SIZE(A))+A) using a partial fraction expansion.
 127
 128 function s = logm_pade_pf (A, m)
 129   [nodes, wts] = gauss_legendre (m);
 130   ## Convert from [-1,1] to [0,1].
 131   nodes = (nodes+1)/2;
 132   wts = wts/2;
 133
 134   n = length (A);
 135   s = zeros (n);
 136   for j = 1:m
 137     s += wts(j)*(A/(eye (n) + nodes(j)*A));
 138   endfor
 139 endfunction
 140
 141 ######################################################################
 142 ## GAUSS_LEGENDRE  Nodes and weights for Gauss-Legendre quadrature.
 143 ##   [X,W] = GAUSS_LEGENDRE(N) computes the nodes X and weights W
 144 ##   for N-point Gauss-Legendre quadrature.
 145
 146 ## Reference:
 147 ## G. H. Golub and J. H. Welsch, Calculation of Gauss quadrature
 148 ## rules, Math. Comp., 23(106):221-230, 1969.
 149
 150 function [x, w] = gauss_legendre (n)
 151   i = 1:n-1;
 152   v = i./sqrt ((2*i).^2-1);
 153   [V, D] = eig (diag (v, -1) + diag (v, 1));
 154   x = diag (D);
 155   w = 2*(V(1,:)'.^2);
 156 endfunction
 157
 158
 159 %!assert(norm(logm([1 -1;0 1]) - [0 -1; 0 0]) < 1e-5);
 160 %!assert(norm(expm(logm([-1 2 ; 4 -1])) - [-1 2 ; 4 -1]) < 1e-5);
 161 %!assert(logm([1 -1 -1;0 1 -1; 0 0 1]), [0 -1 -1.5; 0 0 -1; 0 0 0], 1e-5);
 162 %!assert (logm (expm ([0 1i; -1i 0])), [0 1i; -1i 0], 10 * eps)
 163
 164 %% Test input validation
 165 %!error logm ();
 166 %!error logm (1, 2, 3);
 167 %!error <logm: A must be a square matrix> logm([1 0;0 1; 2 2]);
 168
 169 %!assert (logm (10), log (10))
 170 %!assert (full (logm (eye (3))), logm (full (eye (3))))
 171 %!assert (full (logm (10*eye (3))), logm (full (10*eye (3))), 8*eps)