diff -r eca2af8a92bc ChangeLog --- a/ChangeLog Thu Apr 29 08:00:42 2010 -0400 +++ b/ChangeLog Thu Apr 29 21:27:46 2010 -0400 @@ -1,3 +1,8 @@ +2010-04-29 Richard Guy + + * new logm.m file + * cleared logm.m, sqrtm.m from PROJECTS + 2010-04-26 Rik * configure.ac: fix bug with shell HERE document introduced in diff -r eca2af8a92bc PROJECTS --- a/PROJECTS Thu Apr 29 08:00:42 2010 -0400 +++ b/PROJECTS Thu Apr 29 21:27:46 2010 -0400 @@ -19,8 +19,6 @@ --------- Numerical: --------- - - * Improve logm, and sqrtm. * Improve complex mapper functions. See W. Kahan, ``Branch Cuts for Complex Elementary Functions, or Much Ado About Nothing's Sign diff -r eca2af8a92bc scripts/linear-algebra/logm.m --- a/scripts/linear-algebra/logm.m Thu Apr 29 08:00:42 2010 -0400 +++ b/scripts/linear-algebra/logm.m Thu Apr 29 21:27:46 2010 -0400 @@ -1,4 +1,5 @@ -## 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. ## @@ -17,19 +18,117 @@ ## . ## -*- 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 -function B = logm (A) +## Author: Richard T. Guy with significant code modified from mftoolbox (D. Higham) +## Created: April 29, 2010 - if (nargin != 1) - print_usage (); - endif +function [S,iters] = logm(A,opt_iters) - [V, D] = eig (A); - B = V * diag (log (diag (D))) * inv (V); +if nargin == 0 + print_usage(); + return; +elseif nargin < 2 + opt_iters = 100; +elseif nargin > 2 + print_usage(); + return; +end + +if !issquare(A) + error("logm: argument must be a square matrix.\n"); + return; +end + +[U,S] = schur(A); + +if(any(diag(S) < 0)) + warning("Matrix has nonegative eigenvalues. Principal matrix logarithm is not defined. Computing non-principal logarithm."); +end + + +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. +end + +if k >= opt_iters + warning("Maximum number of square roots exceeded. Results may still be accurate."); +end + +S = logm_pade_pf(S-eye(size(S)),16); +S = 2^k * U * S * U'; + +if nargout == 2 + iters = k; +end 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)); +end + +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]);