######################################################################## ## ## Copyright (C) 2008-2021 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or . ## ## This file is part of Octave. ## ## Octave is free software: you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation, either version 3 of the License, or ## (at your option) any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ## GNU General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, see ## . ## ######################################################################## ## -*- texinfo -*- ## @deftypefn {} {@var{s} =} logm (@var{A}) ## @deftypefnx {} {@var{s} =} logm (@var{A}, @var{opt_iters}) ## @deftypefnx {} {[@var{s}, @var{iters}] =} logm (@dots{}) ## Compute the matrix logarithm of the square matrix @var{A}. ## ## The implementation utilizes a Pad@'e approximant and the identity ## ## @example ## logm (@var{A}) = 2^k * logm (@var{A}^(1 / 2^k)) ## @end example ## ## The optional input @var{opt_iters} is the maximum number of square roots ## to compute and defaults to 100. ## ## The optional output @var{iters} is the number of square roots actually ## computed. ## @seealso{expm, sqrtm} ## @end deftypefn ## Reference: N. J. Higham, Functions of Matrices: Theory and Computation ## (SIAM, 2008.) ## ## Author: N. J. Higham ## Author: Richard T. Guy function [s, iters] = logm (A, opt_iters = 100) if (nargin == 0) print_usage (); endif if (! issquare (A)) error ("logm: A must be a square matrix"); endif if (isscalar (A)) s = log (A); return; elseif (strfind (typeinfo (A), "diagonal matrix")) s = diag (log (diag (A))); return; endif [u, s] = schur (A); if (isreal (A)) [u, s] = rsf2csf (u, s); endif eigv = diag (s); n = rows(A); tol = n * eps (max (abs (eigv))); real_neg_eigv = (real (eigv) < -tol) & (imag (eigv) <= tol); if (any (real_neg_eigv)) warning ("Octave:logm:non-principal", "logm: principal matrix logarithm is not defined for matrices with negative eigenvalues; computing non-principal logarithm"); endif real_eig = ! any (real_neg_eigv); if (max (abs (triu (s,1))(:)) < tol) ## Will run for Hermitian matrices as Schur decomposition is diagonal. ## This way is faster and more accurate but only works on a diagonal matrix. logeigv=log (eigv); logeigv(isinf (logeigv))=-log (realmax ()); s = u * diag (logeigv) * u'; iters = 0; else k = 0; ## Algorithm 11.9 in "Function of matrices", by N. Higham theta = [0, 0, 1.61e-2, 5.38e-2, 1.13e-1, 1.86e-1, 2.6429608311114350e-1]; p = 0; m = 7; while (k < opt_iters) tau = norm (s - eye (n),1); if (tau <= theta (7)) p += 1; j(1) = find (tau <= theta, 1); j(2) = find (tau / 2 <= theta, 1); if (j(1) - j(2) <= 1 || p == 2) m = j(1); break; endif endif k += 1; s = sqrtm (s); endwhile if (k >= opt_iters) warning ("logm: maximum number of square roots exceeded; results may still be accurate"); endif s -= eye (n); if (m > 1) s = logm_pade_pf (s, m); endif s = 2^k * u * s * u'; if (nargout == 2) iters = k; endif endif ## Remove small complex values (O(eps)) which may have entered calculation if (real_eig && isreal (A)) s = real (s); endif endfunction ################## ANCILLARY FUNCTIONS ################################ ###### Taken from the mfttoolbox (GPL 3) by D. Higham. ###### Reference: ###### D. Higham, Functions of Matrices: Theory and Computation ###### (SIAM, 2008.). ####################################################################### ##LOGM_PADE_PF Evaluate Pade approximant to matrix log by partial fractions. ## Y = LOGM_PADE_PF(A,M) evaluates the [M/M] Pade approximation to ## LOG(EYE(SIZE(A))+A) using a partial fraction expansion. function s = logm_pade_pf (A, m) [nodes, wts] = gauss_legendre (m); ## Convert from [-1,1] to [0,1]. nodes = (nodes+1)/2; wts /= 2; n = length (A); s = zeros (n); for j = 1:m s += wts(j)*(A/(eye (n) + nodes(j)*A)); endfor endfunction ###################################################################### ## GAUSS_LEGENDRE Nodes and weights for Gauss-Legendre quadrature. ## [X,W] = GAUSS_LEGENDRE(N) computes the nodes X and weights W ## for N-point Gauss-Legendre quadrature. ## Reference: ## G. H. Golub and J. H. Welsch, Calculation of Gauss quadrature ## rules, Math. Comp., 23(106):221-230, 1969. function [x, w] = gauss_legendre (n) i = 1:n-1; v = i./sqrt ((2*i).^2-1); [V, D] = eig (diag (v, -1) + diag (v, 1)); x = diag (D); w = 2*(V(1,:)'.^2); endfunction %!assert (norm (logm ([1 -1;0 1]) - [0 -1; 0 0]) < 1e-5) %!test %! warning ("off", "Octave:logm:non-principal", "local"); %! assert (norm (expm (logm ([-1 2 ; 4 -1])) - [-1 2 ; 4 -1]) < 1e-5); %!assert (logm ([1 -1 -1;0 1 -1; 0 0 1]), [0 -1 -1.5; 0 0 -1; 0 0 0], 1e-5) %!assert (logm (10), log (10)) %!assert (full (logm (eye (3))), logm (full (eye (3)))) %!assert (full (logm (10*eye (3))), logm (full (10*eye (3))), 8*eps) %!assert (logm (expm ([0 1i; -1i 0])), [0 1i; -1i 0], 10 * eps) %!test <*60738> %! A = [0.2510, 1.2808, -1.2252; ... %! 0.2015, 1.0766, 0.5630; ... %! -1.9769, -1.0922, -0.5831]; %! warning ("off", "Octave:logm:non-principal", "local"); %! assert (expm (logm (A)), A, 40*eps); %!assert (expm (logm (ones (3))), ones (3), 100 * eps); ## Test input validation %!error logm () %!error logm ([1 0;0 1; 2 2])