# HG changeset patch # User Jaroslav Hajek # Date 1280301761 -7200 # Node ID cace99cb01abf85c7ada302f215895fe89c5a1bb # Parent 1e6664326d327c6cf97d483fdea6271233eae6a0 rewrite logm (M. Caliari, R.T. Guy) diff -r 1e6664326d32 -r cace99cb01ab scripts/ChangeLog --- a/scripts/ChangeLog Tue Jul 27 10:46:11 2010 -0700 +++ b/scripts/ChangeLog Wed Jul 28 09:22:41 2010 +0200 @@ -1,3 +1,7 @@ +2010-07-28 Jaroslav Hajek + + * linear-algebra/logm.m: Rewrite. Thanks to M. Caliari and R. T. Guy. + 2010-07-26 Rik * deprecated/complement.m, deprecated/intwarning.m, general/arrayfun.m, general/circshift.m, general/colon.m, general/common_size.m, diff -r 1e6664326d32 -r cace99cb01ab scripts/linear-algebra/logm.m --- a/scripts/linear-algebra/logm.m Tue Jul 27 10:46:11 2010 -0700 +++ b/scripts/linear-algebra/logm.m Wed Jul 28 09:22:41 2010 +0200 @@ -1,35 +1,133 @@ -## Copyright (C) 2003, 2005, 2006, 2007 John W. Eaton +## Copyright (C) 2010 Richard T. Guy +## Copyright (C) 2010 Marco Caliari +## Copyright (C) 2008 N.J. Higham ## ## 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 {Function File} {} logm (@var{a}) -## Compute the matrix logarithm of the square matrix @var{a}. Note that -## this is currently implemented in terms of an eigenvalue expansion and -## needs to be improved to be more robust. +## @deftypefn {Function File} {[@var{s}, @var{iters}] =} logm (@var{a}, @var{opt_iters}) +## Compute the matrix logarithm of the square matrix @var{a}. Utilizes a Pade +## approximant and the identity +## +## @code{ logm(@var{a}) = 2^k * logm(@var{a}^(1 / 2^k)) }. +## +## Optional argument @var{opt_iters} is the number of square roots computed +## and defaults to 100. Optional output @var{iters} is the number of square +## roots actually computed. +## ## @end deftypefn -function B = logm (A) +## Reference: N. J. Higham, Functions of Matrices: Theory and Computation +## (SIAM, 2008.) +## - if (nargin != 1) +function [s, iters] = logm (a, opt_iters) + + if (nargin == 0) + print_usage (); + elseif (nargin < 2) + opt_iters = 100; + elseif (nargin > 2) print_usage (); endif - [V, D] = eig (A); - B = V * diag (log (diag (D))) * inv (V); + if (! issquare (a)) + error ("logm: argument must be a square matrix."); + endif + + [u, s] = schur (a); + + if (isreal (a)) + [u, s] = rsf2csf (u, s); + endif + + if (any (diag (s) < 0)) + warning ("Octave:logm:non-principal", + ["logm: Matrix has nonegative eigenvalues. Principal matrix logarithm is not defined.", ... + "Computing non-principal logarithm."]); + endif + + 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 (size (s)),1); + if (tau <= theta (7)) + p = 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 = 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 = logm_pade_pf (s - eye (size (s)), m); + + s = 2^k * u * s * u'; + + if (nargout == 2) + iters = k; + 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 = wts/2; + + n = length(a); + s = zeros(n); + for j=1:m + s = 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); +%!assert(norm(expm(logm([-1 2 ; 4 -1])) - [-1 2 ; 4 -1]) < 1e-5); +%! +%!error logm (); +%!error logm (1, 2, 3); +%!error logm([1 0;0 1; 2 2]);