diff -r d95a6e491a6c ChangeLog --- a/ChangeLog Wed Jun 09 21:44:21 2010 +0200 +++ b/ChangeLog Wed Jun 09 22:08:49 2010 -0400 @@ -1,3 +1,7 @@ +2010-06-09 Richard Guy + + * new logm.m file + 2010-05-18 Jaroslav Hajek * NEWS: Update. diff -r d95a6e491a6c PROJECTS --- a/PROJECTS Wed Jun 09 21:44:21 2010 +0200 +++ b/PROJECTS Wed Jun 09 22:08:49 2010 -0400 @@ -20,7 +20,7 @@ Numerical: --------- - * Improve logm, and sqrtm. + * Improve and sqrtm. * Improve complex mapper functions. See W. Kahan, ``Branch Cuts for Complex Elementary Functions, or Much Ado About Nothing's Sign diff -r d95a6e491a6c scripts/linear-algebra/logm.m --- a/scripts/linear-algebra/logm.m Wed Jun 09 21:44:21 2010 +0200 +++ b/scripts/linear-algebra/logm.m Wed Jun 09 22:08:49 2010 -0400 @@ -1,35 +1,126 @@ -## Copyright (C) 2003, 2005, 2006, 2007 John W. Eaton +## Copyright (C) 1993, 1995, 1996, 1997, 1999, 2000, 2005, 2006, 2007 +## John W. Eaton ## ## 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. +## +## The Pade approximation was adapted from the mftoobox by D. Higham. +## +## Reference: D. Higham, Functions of Matrices: Theory and Computation +## (SIAM, 2008.) +## ## @end deftypefn + +## Author: Richard T. Guy with significant code modified from mftoolbox (D. Higham) +## Created: April 29, 2010 -function B = logm (A) +function [S,iters] = logm(A,opt_iters) + +if (nargin == 0) + print_usage(); + return; +elseif (nargin < 2) + opt_iters = 100; +elseif (nargin > 2) + print_usage(); + return; +endif - if (nargin != 1) - print_usage (); - endif +if (! issquare(A)) + error("logm: argument must be a square matrix.\n"); + return; +endif - [V, D] = eig (A); - B = V * diag (log (diag (D))) * inv (V); +[U,S] = schur(A); + +if (any(diag(S) < 0)) + warning("Matrix has nonegative eigenvalues. Principal matrix logarithm is not defined. Computing non-principal logarithm."); +endif + + +k = 0; +## Magic number from Higham. +while (k < opt_iters && norm(S - eye(size(S)),1) > 2.6429608311114350e-001) + k = k+1; + S = sqrtm(S); # TODO: There are faster ways to do this. +endwhile + +if (k >= opt_iters) + warning("Maximum number of square roots exceeded. Results may still be accurate."); +endif + +S = logm_pade_pf(S-eye(size(S)),16); +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.). +####################################################################### + +function S = logm_pade_pf(A,m) + +##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. + +[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 + +###################################################################### + +function [x,w] = gauss_legendre(n) + +## 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. + +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]);s